Ultra-coherent nanomechanical resonators based on inverse design

Engineered micro- and nanomechanical resonators with ultra-low dissipation constitute a promising platform for various quantum technologies and foundational research. Traditionally, the improvement of the resonator’s performance through nanomechanical structural engineering has been driven by human intuition and insight. Such an approach is inefficient and leaves aside a plethora of unexplored mechanical designs that potentially achieve better performance. Here, we use a computer-aided inverse design approach known as topology optimization to structurally design mechanical resonators with optimized performance of the fundamental mechanical mode. Using the outcomes of this approach, we fabricate and characterize ultra-coherent nanomechanical resonators with, to the best of our knowledge, record-high Q ⋅ f products for their fundamental mode (where Q is the quality factor and f is the frequency). The proposed approach - which can also be used to improve phononic crystals and coupled-mode resonators - opens up a new paradigm for designing ultra-coherent micro- and nanomechanical resonators, enabling e.g. novel experiments in fundamental physics and extreme sensing.

Topology optimization is a computational morphogenesis procedure widely applied in engineering to determine the best possible structural design and material distributions within a prescribed design domain to maximize a set of performance targets [21]. Examples include the maximization of the structural stiffness of an object under certain design and manufacturing constraints to determine the optimal design of a full-scale aeroplane wing [22] or a girder of a suspension bridge [23], and the maximization of light concentration to develop the optimal design of nanophotonic resonators [24].
The basic strategy of topology optimization is to define a design domain in which material can be distributed. Material is being added to or removed from this domain, and founded on a physical model for the system, a gradient-based computational method is used to optimize the figure-of-merit. Through iterations, material is gradually redistributed towards the optimal design for which the figure of merit is either maximized or minimized, depending on the problem to be solved.
We use topology optimization to optimally design a nanomechanical resonator towards maximizing its Q·f (Qf) product [25][26][27]. Previously, improving the resonator's performance has been done through a combination of human intuition and trial-and-error based on experience and approximative analytical expression for the * denho@fysik.dtu.dk † ulrik.andersen@fysik.dtu.dk different dissipation mechanisms of the resonator. Such an intuition-based approach has recently led to impressive progress in increasing the Qf product of mechanical resonators by using a combination of dissipation dilution [6], soft-clamping [10], thin-clamping [28] and strain engineering [11]. Despite these recent successes, the approach inevitably leaves out many, possibly counter-intuitive, designs that might exhibit superior behavior. Topology optimization counteracts this problem as it directly develops the optimized structure under given initial design constraints and loss models with no geometrical preassumptions. Aiming at maximizing the Qf product of the fundamental mode of a nanomechanical resonator suitable for opto-mechanical experiments, we consider the preconstrained structure illustrated in Fig. 1a. It comprises an area of 700 × 700 µm 2 with a single pad of size 100 × 100 µm 2 (that allows for the interaction with light via radiation pressure force) and a narrow frame of 5 µm to ease fabrication. The remaining space is free to evolve through topology optimization. Furthermore, we assume that the resonator is made of pre-stressed silicon nitride with a thickness of 50 nm. The pre-stressed resonator is numerically discretized using finite (quadrilateral shell) element method, see details in the method section.
Two damping mechanisms associated with intrinsic losses and phonon tunneling losses have been included in the model (Fig. 1b). The intrinsic losses are modeled as bending losses in the form of hysteretic damping, i.e. using a lossy Young's modulus. The phonon tunneling loss (PTL) associated with radiation of the phonons into the substrate was modeled by coupling the boundary outof-plane displacements to continuously distributed lossy springs (illustrated in Fig. 1a) while keeping all the other degrees of freedom fixed. We calibrate the resonator models against previous measurements and from this consider five different cases associated with different ratios between intrinsic and photon tunneling loss, denoted D1-D5. For D1 and D5 the system is purely limited by intrinsic loss and by phonon tunneling loss, respectively, while for D2-D4 the ratio is gradually changed. The exact ratios can be found in the method section. To illustrate the iterative procedure of the topology optimization, in removing buckling-prone features and smoothing sharp features to prevent high tensile stresses at the boundaries.
We note, interestingly, that the optimised design of D1 is similar to the membrane design suggested and experimentally tested in Ref. [29] but using a completely different approach. They arrive at this geometry based on considerations on soft-clamping using a hierarchical design concept [12].
The post-processed designs were patterned on highstress (σ 0 ≤ 1.2 GPa) silicon nitride with a thickness of 12-50 nm grown by low pressure chemical vapor deposition on a silicon wafer. We release the resonators by backetching the silicon substrate in a window of 1.4×1.4 mm 2 (see methods). The fabricated structures are shown in Fig. 3a. To measure the mechanical frequency and quality of the fundamental mode, ring-down measurements were carried out in vacuum (pressure less than 10 −7 mbar) at room temperature using high-sensitivity fiberbased homodyne detection (see methods). An example of a ring-down measurement of a fundamental mode of frequency 240 kHz exhibiting an amplitude ring-down time of about 160s is illustrated in Fig. 2a. This corresponds to a Q factor of 1.18 ± 0.01 × 10 8 and a Qf product of 2.83×10 13 . We also present an example of a thermal noise spectrum including some higher-order modes in Fig. 2b.
We performed ring-down measurements of the fundamental mode of 967 devices that include all the topologically optimized resonators, D1-D5, as well as the conventional non-optimized trampoline resonator [8,9] which is used as reference structure. A collection of some of our measurements on frequency, quality factor and Qf product is presented in Fig. 3. It is clear from these measurements that the topologically optimized resonators are superior to the reference trampolines and that they are all deeply into the regime where the resonator is able to undergo coherent oscillations (corresponding to Qf > 6 × 10 12 Hz) as required for quantum coherent experiments [1].
To understand the limitations of their performance, we fit our best attained result to a theory for the intrinsic and phonon tunneling losses. As the intrinsic loss )UHTXHQF\>N+]@ Overview of the measured frequencies and quality factors across all designs and thicknesses together with selected microscope images. Solid lines correspond to theory fitted to the measured frequencies and best attained quality factors. Dotted and dashed lines are associated with the phonon tunneling and intrinsic loss contributions, respectively. For some designs, the theory curves for the phonon tunnelling loss is not visible on the plots. b. Quality factor plotted against frequency for the five best samples for each design. The shaded area marks the parameter regime in which the resonator may undergo quantum coherent oscillations at room temperature. The inset is the nominal trampoline design fabricated in this work and used as a reference [8].
∆W is mainly dominated by the clamping losses near the boundaries, we directly estimate these losses from the expression ∆W = πφ 12 where h, E, ν are the thickness, Young's modulus and Poisson's ratio of the resonator material [10]. u(x, y) is the mode shape and the loss angle is modelled as φ = 1/(hβ) where β is related to the intrinsic damping at the surface. The mode profiles of all designs are simulated using the COMSOL Multiphysics package with the results shown in Fig. 1c. Phonon tunneling losses are simulated by coupling the resonator to the substrate using a spring with the complex spring constant defined as where S, m , n s and Q s are the substrate's area, mass per unit area, modal density and intrinsic quality factor (see methods). These two loss contributions (intrinisic and phonon tunnelling losses) are then adjusted to match the best experimental data using the loss factors β and Q s as fitting parameters. We find β = (2.93 ± 0.19) × 10 11 m −1 and Q s = (1.27 ± 0.31) × 10 5 , and the resulting theory curves for all designs are shown in Fig. 3a where dotted (dashed) lines correspond to phonon tunneling (intrinsic) losses while the total contribution is represented by solid lines. It is clear that the best performing resonators of all five designs are mainly dominated by intrinsic losses. However, for some resonators we observe markedly lower performance which we attribute to a near-resonant coupling to the substrate modes, consequently leading to significantly higher phonon tunnelling losses which eventually become the dominating loss factor. This randomized coupling to the substrate modes can be circumvented by inserting a damping shield encapsulating the resonator [30]. We highlight the source of intrinsic losses by plotting the bending loss distribution of design D1 in Fig. 4a. First we note that there is a significant amount of bending loss near the boundaries (as highlighted by the inset) and near the intersection between the circular frame and the tethers. This dilution of the bending loss into two areas (resulting from the strong mode confinement) is likely the origin of the quality enhancement, and is similar to the effect observed in resonators based on fractal structures [29]. The observed bending loss at the central pad is due to its low stress leading to a locally reduced stiffness and consequently, sharper bending. In Fig. 4b we Bending loss distribution of D1 on a logarithmic scale. The inset highlights the high bending losses at the boundary. b. Static von Mises distribution. The bars indicate the direction of the first principal stress component. c.
Numerically predicted intrinsic quality factor Qint partitioned into boundary (Q bound ) and distributed (Q dist ) bending losses with illustrate the stress distribution from which we observe a large stress component on the circular frame perpendicular to the tether. The wavelength predicted by the stress and frequency is around 2 mm which is larger than the dimensions of the resonator. Therefore, it cannot exist on the circular frame resulting in mode confinement and dilution of losses. Finally, we compared the amount of boundary bending losses (localized along the outer boundary) to the estimated amount of distributed bending losses (far away from the boundary) as shown in Fig. 4c. It is clear that the resonator is limited by the boundary losses. Micromechanical oscillators with a Qf product of more than 10 13 for the fundamental mode will have a number of intriguing applications in quantum optomechanics and precision sensing. One of the main requirements in quantum optomechanics, e.g. for cooling the oscillator to the quantum ground state, interrogating macroscopic quantum superpositions and entangling different systems, is that the decoherence time exceeds the mechanical oscillation period. This translates into the requirement that Qf > k B T / = 6 × 10 12 Hz at room temperature (where k B is Boltzman's constant, is the reduced Planck's constant and the temperature is T = 300 K) [1, 2]. While most of the resonators fulfill this requirement, our best performing device yields around 4 coherent oscillations which is the largest number ever reported for the fundamental mode of a membrane at room temperature. Our devices will also exhibit exceptional performance in force sensing measurement as for example used in magnetic resonance force microscopy of electron and nuclear spins [31]. In such measurements the sensitivity is limited by the thermal noise ( 4mk B T 2πf Q where m is the mass) which we find to be at 10 aN/ √ Hz for the best devices which is significantly beyond what is attainable with currently available room temperature force microscopes.
The topology optimization method, that we have here employed to maximize the Qf product of the fundamental mode of a membrane, is applicable to many other similar morphogenesis problems in engineering of highperformance micro-and nanomechanical resonators. It can for example be applied to finding the optimal structure for maximizing the dissipation dilution effect -and thus the Qf product -in phononic crystal resonators [10,11] where Q factors of around one billion (at low temperature) have already been achieved without topology optimization [32]. Another interesting avenue for new studies using our methodology is to optimize other application-specific parameters instead of the Qf product. An example is the optimization of the co-operativity parameter associated with the coupling of a specifically functionalized mechanical oscillator to spins [33], light [34] or charges [19] with the aim of significantly enhancing quantum transduction or sensing. Finally, it is also possible to optimize structures with addtional constraints, either structural constrains enabling certain applications or parameter constraints, e.g. fixing the mass to a large value with the aim of maximizing the coupling to gravity as required for example for interrogating the quantum nature of gravity. Our methodology thus has the potential to revolutionize the way nano-and micro-mechanical systems are being designed enabling radically new applications and fundamental explorations.

A. Topology optimization implementation
We employed a density-based topology optimization approach [21] to design ultrahigh coherent resonators. The basic methodology and the detailed optimization formulation are described in the following.
Prestressed membrane resonators are simulated using finite element methods with the 4-node MITC (Mixed Interpolation of Tensorial Components) quadrilateral shell element [35]. The mechanical dynamic problem is solved in two steps: 1) Establish static equilibrium of a prestressed membrane resonator under prescribed stress; 2) Identify resonating modes using linear eigenvalue analysis. The FE equations are stated in discrete form as, Here F 0 is the equivalent force vector resulting from a prestress σ 0 , K 0 represents the linear stiffness matrix and K σ (U 0 ) represents the initial stress stiffness matrix that depends on the displacement U 0 of the prestress problem in Eq. (3). C and M denote damping and mass matrixes, ω j and φ j are the angular frequency and modal profile of the j-th resonating mode and i = √ −1 is the imaginary unit.
The damping matrix, C, covers intrinsic and phonon tunneling losses. The intrinsic losses are considered via a relaxation mechanism described by a complex-valued Young's modulusẼ = (1 + iη s ) E. The phonon tunneling losses are modeled using damped springs distributed along the boundary with a total stiffness ofk b = (1 + iη b ) k b and k b = 8.315 × 10 7 kN/m 2 . The detailed calculation formulations of quantities in Eqs (3) and (4) can be found in [27]. The quality factor and frequency of the j-th resonating mode are calculated by In the density-based topology optimization approach, an element-wise design variable, x e ∈ [0, 1], is introduced to indicate the material occupation in element e. To avoid checkerboard pattern and mesh dependence [21] and enhance design discreteness, the design variables are first filtered using a density filter [36] and then smoothly projected using a hyperbolic tangent threshold function [37], given as Here,x e is the filtered design variable, y k are the center coordinates of element k. v k and x k are the corresponding volume and design variable of element k, respectively. N e is the neighborhood of element e within a certain filter radius specified by N e = {k| x k − y e ≤ r}, and w e (y k ) = r − y k − y e .x e is the projected design variable of element, e. When β 1 is large,x e ≈ 1 if x e > η representing Si 3 N 4 , andx e ≈ 0 ifx e < η indicating void. The projection suppresses gray element density regions induced by the density filter when β 1 is sufficiently large and ensures black-white designs when the optimization converges. Moreover, it mimics the manufacturing process and the manufacturing errors can be taken into accounts in the optimization by choosing different thresholds, η, as discussed later. The Young's modulus of element e is directly related to the projected design variable using the Rational Approximation of Material Properties (RAMP) [38] and the mass density is linearly interpolated as Spurious modes caused by inappropriate stiffness-tomass ratios in low-density regions are suppressed by setting E 0 = 10 −6 E and ρ 0 = 10 −7 ρ to represent void in this study. Wrinkling-like instabilities in low-density regions are alleviated using a displacement interpolation with detailed formulations presented in [27].
To enhance the design robustness with respect to manufacturing errors and impose a minimal length scale in the nominal design, a three-case robust formulation is employed [37]. Three design realizations are generated to mimic an eroded, normal and dilated manufacturing processes. The optimization problem for designing ultrahigh coherent resonators is formulated to maximize the Qf product of the fundamental mode for the worst case of the three design realizations, subjected to frequency constraints and a volume fraction constraint, given as The three design realizations are generated using η ∈ {0.55, 0.5, 0.45} with a filter radius of r = 15 µm. This corresponds to a minimal feature size of 6.7 µm in both solid and void regions of the nominal design. The prescribed frequency lower bound and volume fraction upper bound are f * = 240 kHz and v * = 0.5. Gradients of the objective and constraint functions are calculated using the adjoint sensitivity analysis and the chain rules [27,36,37]. The design variables are iteratively updated using the deterministic mathematical programming approach, Method of Moving Asymptotes (MMA) [39] based on the gradients of the objective and constraints. β 1 is updated until the convergence criterion is satisfied by β

B. Fabrication
We deposit stoichiometric silicon nitride onto a 100 mm single-crystal silicon wafer of 500 µm thickness using low-pressure chemical vapour deposition. This is followed by spincoating photoresist onto the wafer and transfer of the different resonator designs using UV-lithography. The photoresist is developed and the silicon nitride is etched in these regions by means of reactive ion etching. Residual photoresist is removed using oxide plasma, and finally, the trampolines are released in potassium hydroxide at 80 • C followed by cleaning in hydrochloric acid and sulfiric acid mixed with ammonium persulfate.

C. Characterization
Measurements of the frequency and quality factor of the resonators are performed using optical interferometry driven by a laser with a wavelength of 1550 nm. The laser beam is reflected off the vibrating membrane (located inside a vacuum chamber at low pressure < 10 −7 mBar), and the resulting phase shift is detected with highsensitivity using a phase-locked homodyne detector and recorded with a spectrum analyzer. Excitation of the mechanical oscillator is done by modulating the intensity of the laser at the resonance frequency. Once excited, the modulation is switched off and the amplitude decay is subsequently measured.

D. Phonon tunneling loss model
The implemented phonon tunneling loss model used in evaluating the designs post topology optimization, is derived by treating the out-of-plane forces induced by the resonator onto the substrate as a point source on a large but finite silicon substrate. This assumption is valid as long as the characteristic size of the membrane as well as the substrate thickness are much smaller than the wavelength of the radiating wave in the substrate. This is the case for the treated designs when considering their fundamental mode. Furthermore, at these low frequencies, Kirchoff-Love plate theory can be used to model the vibrations induced into the substrate.
In the point-source assumption the amplitude of vibrations u s at an excited point (x, y) on the finite substrate is related to the excitation force F by where S is the area of the substrate and m = ρ s h s is the mass per unit area of the substrate with ρ s and h s being the substrate density and thickness, respectively [40]. ψ n (x, y) are the eigen-modes of the substrate with eigen-frequencies ω n and Λ n = ψ 2 n (x, y) dxdy/S. To simplify the model an effective spatial overlap is assumed, ψ 2 n (x, y) → Λ n . Furthermore, a new variable is introduced describing the spectral distance ∆ω n between the excitation frequency and substrate modes defined as ω n = ω + ∆ω. When only the closest substrate is considered, an effective spring constant describing the coupling between a resonator and the substrate can be derived as When a large substrate is used (like a 100 mm silicon wafer) the spectral distance to the closest mode will be difficult to estimate, due to the high density of modes. Instead, a stochastic approach is used where ∆ω ≡ XY 2 . Here, X is a uniformly distributed variable between 0 and (IIHFWLYHWHQVLOHSUHVWUHVV>*3D@ , where n s is the modal density of the substrate. E s and ν s are the substrate's Young's modulus and Poisson's ratio. This leads to a distribution of possible stiffness for the spring. Due to the assumption that the wavelength of the excited wave in the substrate is much larger than the resonator dimensions, we can treat the outer boundary of the resonator as a single rigid frame with only 1 degree of freedom normal to the plane of the resonator. This frame is then coupled to a reference via a spring defined by k PTL . This enables easy implementation of the PTL model into finite element models. For simulations the expectation value of k PTL is used (see Eq (2)). Note that this implementation is only valid for the fundamental mode (and possibly a few higher-order modes) of the resonator.

E. Stress-thickness dependency
To model the frequency dependency of the silicon nitride thickness in Fig. 3a, the tensile pre-stress dependency of the thickness is needed. The stress-thickness dependency is believed to be caused by the oxidization layer which introduces a compressive stress contribution onto the silicon nitride film dependent on its thickness [41]. Assuming that the oxidized layer is much smaller than the total film thickness, we model the effect by the expression σ(h) = σ 0 − β σ /h where σ 0 is the asymptotic pre-stress parameter and β σ is a coefficient that determines how fast the pre-stress changes with thickness. We fit these two parameters against data attained from the measurement of tensile pre-stress from 2573 samples of different thicknesses as shown in Fig. 5. The tensile stress was derived by measuring the resonance frequency and comparing to predicted values from finite element simulations noting the f ∝ √ σ dependency. This approach has some inherent uncertainties related to fabrication and the assumptions of the material parameters of silicon nitride. We find σ 0 = 1.235 ± 0.002 GPa and β σ = 4.52 ± 0.06 Pa · m.