Spiderweb nanomechanical resonators via Bayesian optimization: inspired by nature and guided by machine learning

From ultra-sensitive detectors of fundamental forces to quantum networks and sensors, mechanical resonators are enabling next-generation technologies to operate in room temperature environments. Currently, silicon nitride nanoresonators stand as a leading microchip platform in these advances by allowing for mechanical resonators whose motion is remarkably isolated from ambient thermal noise. However, to date, human intuition has remained the driving force behind design processes. Here, inspired by nature and guided by machine learning, a spiderweb nanomechanical resonator is developed that exhibits vibration modes which are isolated from ambient thermal environments via a novel"torsional soft-clamping"mechanism discovered by the data-driven optimization algorithm. This bio-inspired resonator is then fabricated; experimentally confirming a new paradigm in mechanics with quality factors above 1 billion in room temperature environments. In contrast to other state-of-the-art resonators, this milestone is achieved with a compact design which does not require sub-micron lithographic features or complex phononic bandgaps, making it significantly easier and cheaper to manufacture at large scales. Here we demonstrate the ability of machine learning to work in tandem with human intuition to augment creative possibilities and uncover new strategies in computing and nanotechnology.

Major advances in nanotechnology have allowed mechanical resonators to improve dramatically over the last decades. One of the most sought after characteristics for a mechanical resonator is noise isolation from thermal environments, namely at room temperature conditions where thermomechanical noise can dominate. The degree of mechanical isolation is characterized by a resonator's mechanical quality factor, Q m . Typically Q m is defined as the ratio of energy stored in a resonator over the energy dissipated over one cycle of oscillation. Inversely, mechanical quality factors can indicate the dissipation of mechanical noise into a resonator from ambient environments. For mechanical sensors, a resonator's isolation from ambient thermal noise can greatly enhance their ability to detect ultra-small forces, pressures, positions, masses, velocities, and accelerations. For quantum technologies, mechanical quality factor dictates the average number of coherent oscillations a nanomechanical resonator (in the quantum regime) can undergo before one phonon of thermal noise enters the resonator and causes decoherence of its quantum properties [1]. From microchip sensing to quantum networks, cryogenics are conventionally required to counteract thermal noise but enabling these burgeoning technologies to operate in ambient temperatures would have a significant impact on their widespread use.
In room temperature environments, on-chip mechanical resonators with state-of-the-art quality factors have mostly consisted of high-aspect-ratio suspended nanostructures fabricated from tensile thin-films. Silicon nitride (Si 3 N 4 ) films have been the material of choice for their high intrinsic stress, yield strength, temperature stability, chemical inertness and prevalence in nanotechnology. Over the years, researchers have developed improved design principles that manipulate the strain, bending and mode shape in nanomechanical resonators to improve quality factors, which are ultimately limited by bending losses as the resonator oscillates in vacuum. Another important characteristic is the mechanical frequency of a nanomechanical resonator's vibrational mode, f m . For high-precision detectors of fundamental forces like gravity and dark matter [2,3] and quantumlimited commercial sensors [4,5], resonators with high Q m and low f m are a long-standing goal. Minimizing f m /Q m is a key figure-of-merit towards quantumlimited force [6] or acceleration [7] sensitivities (S f , S a ∝ (f m /Q m ) 0.5 ) and for enabling quantum sensing of forces like dark matter [8,9] and gravity [10,11] where low frequency and higher quality factor are advantageous. For phonon-based quantum technologies, a mechanical resonator's vibrational modes are initialized into the quantum regime, where their motion harbors less than one quanta of vibration (phonon) [12,13]. Mechanical resonators in these quantum regimes must have sufficiently high Q m × f m > k B T room /h to suppress the effects of room temperature, T room , thermal noise on their frag-ile quantum properties. While there are only a handful of platforms [14][15][16][17][18][19][20][21] to overcome these stringent requirements on Q m × f m at room temperature, a general goal has been to achieve the highest Q m and lowest f m possible while still maintaining Q m × f m above 6 × 10 12 Hz.
Previously proposed nanomechanical resonators follow strategies largely motivated by one-dimensional analytical models of resonators [23] because they provide easyto-interpret design rules. It is important to note that while silicon nitride has been conventionally used for high-Q resonators, these design principles are valid for nearly any strained thin-film material. These analytical models show that higher strain (i.e. higher stress σ and low Young's modulus E), longer (L) and thinner (t) geometries generally lead to higher quality factors in both nanomechanical membranes and strings. For example, f m /Q m of the double clamped beam's fundamental mode is proportional to √ E/L 3 (L + 1.4t E/σ), when assuming a thin long pre-stressed beam (explicit formulation can be found in the Supporting Information). While increasing the aspect ratio of resonators usually leads to smaller f m /Q m , it also makes them much more challenging to fabricate reliably. Figure 1A illustrates conventional design strategies. When considering the fundamental mode of resonators, mechanical quality factors are typically improved by pre-stressing -a form a strain engineering called dissipation dilution [24][25][26][27] which increases stored energies and lowers dissipation compared to unstressed resonators. It overcomes the fundamental limit of the material's intrinsic damping from its bulk and surface, enabling a higher quality factor by orders of magnitude. Simultaneously, resonator models [23] also explain that for high-aspect ratio prestressed resonators (t << L) the dominating loss that decreases Q m is due to the mode's sharp curvature at the clamped boundary between the oscillating element and the substrate (on which the resonators are fabricated). This observation motivated the use of phononic crystals (i.e. phononic bandgap) which confine a higherorder mode from the clamping regions using a periodic pattern around the mode. Now rather than having large curvature near the edges, the phononic crystals enable a soft-clamping [15,19,21] to reduce the mode's curvature close to the rigid clamp, and thereby eliminating this dissipation mechanism when operating at high order vibration modes as depicted in the schematic in Figure  1A. Phononic crystals enable higher quality factors that approach hundreds of millions, at the cost of operating at higher frequencies and typically requiring higher aspectratio resonators which are more difficult to manufacture reliably.
To illustrate the design space, Figure 1B shows the mechanical quality factor Q m versus vibration frequency f m for a 50 nm thick and 3 mm long Si 3 N 4 beam. The color scale indicates the various out-of-plane vibration modes, from fundamental to higher-order modes with shorter wavelengths. The Supporting Information ("Lessons from string resonators") provides additional in-formation about the bounds shown in the figure. The considered pre-stressed double clamped beam shows the improvement of Q m caused by enhancement of stored energy (intermediate solid line) when compared to the unstressed beam (horizontal solid line at the bottom, Q 0 ) and the improvement caused if perfect soft clamping is achieved around the region shown in Figure 1A (top solid line). These lines help illustrate the design region in the graph's top left corner that remains largely unexplored in current resonators that aim mainly for high quality factor. Previous works have focused on increasing Q m of the fundamental mode with strain engineering [6,20], including design strategies as topology optimization [17,28] or hierarchical designs [16,22]. Here we pursue high quality factors at lower frequencies by following a new approach inspired by nature and guided by machine learning.
Spiderweb designs have unique geometries that make them one of the most well-known and fascinating classes of micro-mechanical structures found in nature. Despite their ubiquitous presence, experts from physics, materials science, and biology are still uncovering the elusive mechanics of spiderwebs that enable them to be remarkably robust vibration sensors [29,30]. Spider silk threads have high toughness and stiffness, reaching yield strengths on the order of a gigapascal -about 5 times higher than steel [31] and about the same as Si 3 N 4 . They are used to create lightweight fibrous web structures which harbor an extraordinary strength-to-weight ratio rarely observed among other structures found in nature or science [32][33][34][35]. Furthermore, in the case of spiders that sense their prey via the webs, these structures are designed to be most sensitive to vibrations emanating from the web and not from surrounding environmental vibrational disturbances [36][37][38][39]. Since their unique sensing capabilities have been relentlessly optimized over millions of years of complex evolutionary competition [40][41][42], spiderwebs stand as a promising starting point for machine learning algorithms to design nanomechanical sensors [43][44][45].
Without making any assumptions about how a spiderweb functions as a vibrational sensor in nature, we propose a web-like structure composed of radial beams, lateral beams, and junctions between them, as shown in Figure 2. Instead of spiderweb threads which are microns thick, we consider highly-stressed Si 3 N 4 that can be as thin as 20 nm while being suspended over several millimeters. The properties of Si 3 N 4 were considered to be E=250 GPa, ν=0.23, ρ=3100 kg/m 3 with an initial released stress of 1.07 GPa based on measurements. The parameterized model shown in Figure 2 includes six design parameters, d, w 1 , w 2 , l 1 , l 2 and N r . The two inner rings (i.e. the rings formed by the lateral beams) were constrained to have at least a distance of 8 µm between them. We also considered even numbers between 4 to 16 for the number of lateral beams per ring, N r . Note that the even number of N r was considered to include not only the symmetric but also the anti-symmetric periodic boundary condition. The beam width at the outer and FIG. 1. A Illustration of the mechanical resonator and target vibration modes for a high quality factor resonator. Unlike fundamental modes [20,22], or higher-order modes [15], the target mode shape for lower-order mode hadn't been discovered. B Quality factor (Qm) versus frequency (fm) for a double clamped 50 nm thick 3 mm long silicon nitride beam. The bottom solid line corresponds to the intrinsic quality factor of Si3N4, the intermediate solid line to the effect of high pre-stress, and the upper solid line to the effect of complete elimination of clamp losses (perfect soft clamping) on the string resonator. The dashed line corresponds to the mechanical decoherence constraint, indicating that the resonator having the quality factor above can complete one full coherent oscillation without a thermal phonon entering the resonator. The figure highlights the target unexplored region of designs with high quality factors and lower-order modes.
inner parts of the radial beams (d o and d i ) were taken to be 4 µm and 1 µm, respectively. The width of the structure at the resonator-substrate interface, w oh , was set to 2 µm (half the maximum beam width) in the finite element model to reflect the inevitable overhang originating from the fabrication process. The fixed boundary condition around the resonator is modeled to reflect the overhang attached to a fixed substrate. Additionally, we gave a 1 µm radius fillet for every corner at the junction of lateral and radial beams. By limiting the model's features such as tether widths and fillets to micron scales (rather than sub-micron), it ensures that these structures can ultimately be defined using photolithography which allow for significantly easier, large-scale fabrication. Finally, the simulation model in the paper considers L, t as 3 mm and 50 nm, respectively. From the simulation, we estimated the mechanical quality factor of the resonator by calculating the dissipation dilution [27,46] of the out-of-plane vibration modes. The quality factors are calculated as, with α and β being defined as, β = u 2 z,xx + u 2 z,yy + 2νu z,xx u z,yy + 2(1 − ν)u 2 z,xy with u z being the out-of-plane displacement during vibration, and σ the stress distribution resulting from static analysis for the initial stress. The comma denotes a partial derivative with respect to that coordinate. Q 0 is the intrinsic quality factor defined as where Q volume is the bulk material loss of Si 3 N 4 , and Q surface is the surface loss that varies linearly with the resonator's thickness. For thin resonators in room temperature, we assume Q 0 ≈ 6900 t/100 nm [47]. Note that α is proportional to the elastic energy in tension, which corresponds to the energy stored in the resonator, and β is proportional to the bending loss. The detailed derivation is provided in the Supporting Information ("Derivation of the quality factor for two-dimensional structures"). With the parameterized model, we aimed to find the highest quality factor considering general types of mode shape below 1 MHz, which is the tenth mode frequency of the same size of the pre-stressed string in Figure 1B. Given that the length L of our spiderweb nanomechanical resonators are limited to 3 mm, this ensures a target frequency in the hundreds of kHz regime.
The choice of optimization algorithm to guide the datadriven design process represents the most crucial part for solving real application problems and depends on the characteristics of the problem and data availability. For example, recently machine learning algorithms have proved their success in material design problems with abundant data [48][49][50]. On the contrary in this case, a new resonator is designed based on a new spiderweb model shown in Figure 2 and therefore, no prior data is available. Trial-and-error experimentation is difficult because conducting a single experiment of a particular design takes several days of fabrication and testing. In addition, fast analytical predictions of the quality factors and vibrational mode frequencies of designs are also not possible due to the complexity of the twodimensional geometry, which also requires to consider various vibration mode shapes. In fact, finite element analyses of these structures are computationally expensive; taking between 10 and 30 minutes using 20 CPU cores of our high performance computing cluster. Despite only considering periodic boundary conditions and simulating only a fraction of the total structure, these long simulation times arise from the fine mesh of elements required by the high aspect ratio structure such that the subtle shape curvatures are captured in small but crucial regions, such as joints. The geometry is meshed with 4 to 6 elements in the beam's width direction using shell elements, and the mesh resolution near the joints is about double that (detail of the finite element model information can be found in Supporting Information). Therefore, a week of computation can only generate data corresponding to less than 1000 design iterations. This is then classified as a data-scarce optimization problem where each new design iteration should be as informative to the design goal as possible. Under these conditions, using data-scarce machine learning to guide the optimization process is particularly effective, as achieved by the Bayesian Optimization method [51][52][53].
Bayesian optimization [54] constructs a machine learning regression model usually from Gaussian processes [55], by predicting model uncertainty and seeking the optimum solution in fewer iterations than competing algorithms [56,57]. Applying the algorithm to new problem domains, which requires new kinds of surrogate models without pre-domain knowledge as in our problem, is especially beneficial [57]. For readers unfamiliar with the topic, the Supporting Information includes a short introduction to the method. In the context of designing the spiderweb nanomechanical resonator, Bayesian optimization is expected to not only explore the design space to find new vibrational modes that induce soft clamping with a compact design, but also use them to reach high quality factors in the low frequency regime for a given resonator size. In this work we used the GPyOpt python implementation of the method [58], and MATLAB for the pre-and post-processing of our spider-web design. The finite element analysis was performed by COMSOL [59]. Figure 3 refers to the optimization history wherein the spiderweb nanomechanical resonator's quality factor is maximized. The process starts with a random search of 40 iterations to train the model (shown in light blue), followed by the Bayesian optimization phase (shown in light red). Figure 3B shows the evolution of the quality factor Q m , while Figure 3D plots the distance from a previous optimized point to the point considered in that iteration (the distance between the normalized input vectors). As seen in Figure 3A, the random search from iterations 1-40 find vibrational modes with the highest quality factors in iteration 26 and 27 which vibrate in the outer lateral ring. When the Bayesian optimization starts at iteration 41 (red markers), it begins to follow the vibrational modes concentrated on the inner lateral beams. Even though iteration 41 starts with a quality factor similar to the one found in the random search (light blue region), the algorithm continues exploiting the optimal design region without going too far from this successful iteration, as can be seen in Figure 3D between iterations 41 and 51, mostly by finely tuning design parameters as can be seen in Figure 3E. In these iterations, the design improves the quality factor by more than 180% compared to the best value obtained from the first 40 random searches and the machine learning model promotes local optimization.
Note that the optimal modes found during the Bayesian optimization process correspond to vibrational modes (iterations 41, 51, 136 in Figure 3A) which harbor only a slight deformation near the clamping points because the major bending elements are in the intermediate ring. Surprisingly, these lateral vibrational modes mimic actual vibrations utilized in spider webs for prey detection [60]. Without encoding any prior knowledge about how spiderwebs function, the machine learning algorithm was able to find how actual spiderwebs work in nature and adapt it to silicon nitride nanostructures. After iteration 51, the algorithm starts exploring the design ; B evolution of the quality factor Qm (the red markers and red line indicate the highest quality factor until that iteration); C the frequency fm; D the distance from a previous optimized point to the point considered in that iteration; and E values of the 6 design parameters at every iteration, providing an idea of how the designs change in the optimization process (4 ≤ Nr ≤ 16, 1 µm ≤ (d, w1, w2) ≤ 4 µm, 0 mm < l1 < l2 < 1.5 mm). The abscissa for B-E is the same and corresponds to the design iterations as the optimization process evolves. The blue region in B-D corresponds to the 40 initial designs that were randomly selected, i.e. before starting the Bayesian optimization process; while the red region corresponds to the Bayesian optimization iterations. space more to search for a better design far from the previous optimal, which can be seen from the high values in Figure 3D. This trade-off between exploration and exploitation is often responsible for the competitive advantage of Bayesian optimization when compared to other algorithms. The green markers in Figure 3 are clear examples of the Bayesian optimization exploring far from previous optima. Note that the lateral beam's vibrational mode did not always have the highest quality factor, as optimum performance arises from the discovery of this new mode in combination with geometric parameters that promote bending and torsion in a particular way, as discussed later. By exploring and exploiting simultaneously to optimize the quality factor of the mechanical resonator, the algorithm reaches a maximum quality factor in iteration 136. Bayesian optimization consistently balances exploitation and exploration, so the solution for a large number of iterations could lead to continuous improvement. The result here considers up to 200 iterations for the optimization, considering a few days for optimization. In the case of having thousands of data sets regardless of the computational cost, the consideration of additional design parameters could be interesting for future work. As seen in Figure 3C, we also found that the highest quality factors Q m occur at lower frequencies, which agrees with the simplified one-dimensional model of the beam resonator. Note that the convergence speed and the optimum result could depend on the initial searching points. The comparison study can be found in the Supporting Information ("Optimization convergence dependency on the initial random points"). Moreover, the optimization process also shows that thinner structures are not necessarily better when using micro-wide tethers -a counter-intuitive finding given that every previ-ous design of nanomechanical resonators had out-of-plane mechanical modes with Q m that benefited from thinner geometries. The detailed discussion can be found in the Supporting Information ("Conversion of energy loss regarding the thickness of the spiderweb nanomechanical resonator").
The optimum spiderweb nanomechanical resonator is predicted to have a quality factor Q m above 1.75 billion at 134.9kHz for a design considering a diagonal size of 3 mm and a thickness of 50 nm. Table I shows the design parameters corresponding to this design. Comparatively, what is striking about this design is that it is able to achieve a quality factor Q m above a billion without requiring any tether widths under a micron. This allows it to be readily defined using large-scale photolithography which further makes manufacturing faster and cheaper. This is extremely beneficial for real applications in that it decreases the resonators' microchip footprint. For thermal management, lower aspect ratio is very beneficial as typically Si 3 N 4 nanoresonators are interfaced with optics for high-precision sensing and quantum application. Although Si 3 N 4 is preferred for its low optical absorption [26,61], even small amounts of optical heating can have deleterious consequences for high precision experiments, and having a smaller aspect ratio allows for enhanced thermal conduction (to the substrate) which scales as t/L. What is important to note is that previous nanomechanical resonators followed a common design paradigm wherein the maximal amplitude of the mode is in the center of the resonator and where the aim is to reduce bending losses from that center to the substrate. Here the algorithm takes an entirely different route by looking at modes that oscillate laterally like a mechanical whispering gallery mode, allowing for tiny distances between amplitude maxima and substrate, providing new insight into the nanomechanical resonator design. Figure 4 provides additional details about the optimum design, where Figure 4A highlights the novel "torsional soft clamping" mechanism found by the datadriven strategy. This design yields an unprecedentedly high quality factor because the resonator vibrates with an out-of-plane deformation that is localized in the inner ring of lateral beams while undergoing torsional deformation of the radial beams. Figure 4B and C also support this observation, where the displacement magnitude (B) and normalized bending loss (C) clearly demonstrate the low displacement and bending losses at the boundary (blue marker) while the joint of the inner lateral ring (yellow marker) undergoes significant deformation. As a consequence of the radial beams' torsional motion, the curvature between the bending lateral beam with the ra- dial beam at the clamping point is not highly concentrated, thus significantly diminishing the clamping losses at these points. Although the torsional motion of the radial beams leads to 70 % of the energy dissipation, it is comparable to the bending losses in the deforming lateral beams. Note that the normalized bending loss in Figure 4C indicates that the bending energy near the joint of the inner lateral beam (yellow marker region) is spread out in the region near the joint, which avoids sharp curvatures that can ultimately limit Q m . Additionally, the blue marker region of Figure 4A shows that the outer lateral beams near the boundary prevent deformation and bending loss from propagating toward the boundary in a subtle way by preventing torsional deformation in the radial beams from propagating to the boundary. Unlike the bending loss density of the string's bending modes [23] (where it is highly concentrated on the clamping region), the torsional bending loss density in the spiderweb resonator does not highly concentrate where the vibration stops. For this reason, simulations without the outer lateral beams also gave a similar quality factor. Nonetheless, the optimized position of the outer ring was used to block the torsion propagation from the inner ring to the chip, enhancing the resonator's isolation from the substrate. The optimized N r is 4, which maximizes the side beam length when all other parameters are the same. This trend shows that the optimized mode aimed to make the vibrating lateral beams as long as possible, as the quality factor of a string resonator [27] increases for a longer beam. Compared with the one-dimensional approach [14,15], which is a subset of our structure by considering N r = 2, the web-like structure in the two-dimensional domain has the potential to achieve higher Q m by exploring novel soft clamping motions compared to the limited number of vibration modes [16] of the one-dimensional structure. The optimized l 2 , which defines the distance between the centre and the inner radial beam, was close to the maximum limit, while the number of radial beams N r was minimized. This trend shows that the optimized mode aimed to make the vibrating lateral beams as long as possible, as the quality factor of a string resonator [27] increases for a longer beam. Also, the w 2 was optimized at 1 µm which represents the minimum limit allowed. Allowing for widths thinner than 1 micron would likely increase the quality factor. While the lateral size, L, was limited to 3 mm, a study of increasing L (discussed in the Supporting Information) also shows that the optimized Q m scales quadratically with the size of the resonator while f m /Q m is significantly reduced. All of these mechanisms work together to achieve a compact design with an ultra-high quality factor at a lower order vibration mode via a novel soft clamping approach that does not require the use of phononic crystals or sub-micron lithographic features. Notwithstanding, up to now this finding remains a computational prediction.
While we have so far described the computational design process which occurred without experimental trial- and-error, the performance of the novel resonator needs to be experimentally validated. To do this we fabricate the optimal resonator and experimentally determine the quality factor of the system from a ringdown measurement in the high vacuum setup shown in Figure 5F to avoid air damping. Ringdown measurements involve using piezoelectric stages to resonantly excite the motion of nanomechanical resonators, stopping the drive, and observing decay of the resonators' motion via interferometric optical readout. The rate of decay of the resonators amplitude gives its rate of energy dissipation and thus its mechanical quality factor. The spiderweb nanomechanical resonator is fabricated on high-stress Si 3 N 4 grown by low-pressure chemical vapor deposition on a silicon wafer ( Figure 5A). The pattern is first written on a resist mask ( Figure 5B), then transferred on the Si 3 N 4 layer with a directional CHF 3 plasma etching ( Figure  5C). After that, the resist mask is removed ( Figure 5D) and the spiderweb nanomechanical resonator is released by a fluorine-based (SF 6 ) dry etching step ( Figure 5E), which does not require any mask or additional cleaning steps, making the fabrication considerably easier, higheryield and higher-quality. Crucially, this fabrication process allows for remarkable agreement between experimental results and idealized simulations that allow us to reliably use the latter as data-points for the machine learning algorithm. A detailed explanation of the fabrica-tion process and the mechanical characterization setup is summarized in the Methods Section. Figure 6A shows a scanning electron microscope image of the fabricated device, where the suspended Si 3 N 4 spiderweb nanomechanical resonator is highlighted in blue, surrounded by dummy Si 3 N 4 islands disconnected from the suspended structure, used to prevent overetching and overexposure (in light gray). The thermomechanical noise spectra as obtained by interferometric optical readout at the center of the spiderweb nanomechanical resonator (orange marker) and the inner lateral beam (pink marker) are plotted in Figure 6C. As shown in the figure, the optimal vibration mode occurring around 133.6 kHz is visible in the spectrum of the inner ring but has no amplitude in the center of the resonator and therefore well confined, contrary to the fundamental mode at 81.8 kHz. This corroborates the presence of the novel soft clamping mode shown in Figure 4. Furthermore the simulated out-of-plane vibration mode frequencies agrees with experiments to around 1% as shown in Figure 6B. The presence of the torsional soft clamping mode is further proved by the ringdown measurement shown in Figure  6D, where the quality factor of the optimized mode at 133.6 kHz is measured to be 1.8 billion, which is in excellent agreement with the computational predictions (1.75 billion) and the highest mechanical quality factor yet measured in this frequency range at room temperature. Compared to the fundamental mode's Q m , it has more than 60 times higher value. In this research, we simulated the Q m considering the dissipation dilution but the acoustic radiation (Q rad ) [16,20] and the loss from gas damping (Q gas ) [25] can also affect the ringdown in the experiment. The excellent match between bending loss simulations and experimental results supports our hypothesis that bending loss is the dominant source of mechanical loss for the spiderweb nanomechanical resonator. While acoustic radiation loss (through the substrate) affects mechanical resonators with motion near the resonator-substrate boundary [16,17,20], it is expected to be negligible in our optimal design because the resonator motion is isolated from the substrate, with the support of the outer lateral beams, as described in Figure 4. At the same time, the gas damping effect can be ignored by performing the measurement under a sufficiently high vacuum of 4.0·10 −9 hPa. The high Q m result is especially striking when considering the short length and larger thickness of the resonator than existing solutions in the literature, making it more practical to fabricate and operate. Figure 6E compares our result with the state-of-the-art nanomechanical resonator's experiment values at room temperature by plotting their Q m , f m and their aspect-ratio t/L via marker area when t and L represent the thickness and the size of each reported resonator, respectively. It supports that our spiderweb resonator has obtained a high Q m not using more chal-lenge fabrication but considering novel vibration mode.
In conclusion, a simulation based data-driven optimization approach was used to design a spiderweb nanomechanical resonator with ultralow dissipation in room temperature environments. Our approach relies on the observation that spiderwebs have evolved over millions of years through evolutionary competition to be remarkable vibration sensors [29,30]. Using silicon nitride as a base material, our machine learning algorithm hitchhikes on this natural optimization, and discovers nanomechanical designs tailored for high precision sensors. While silicon nitride is one of the most widely used thin-films for nanomechanical resonators, the design approach in this work could be extended to other materials such as diamond [68], gallium arsenide [64,69], silicon carbide [70,71], indium gallium phosphide [65,72], fused silica glass [73], silicon [74], phosphorus carbide [75,76], and even superconducting films [77,78]. The enhancement of mechanical quality factor results from the discovery of a soft-clamping mechanism that uses a torsional motion to isolate a nanomechanical mode from ambient thermal noise. This enables high-Q m nanomechanical resonators that have smaller aspect ratios than previous state-of-the-art designs, making them significantly easier, cheaper and faster to manufacture. Our experimental validation demonstrates a new class of mechanical resonators that exhibit mechanical quality factor exceed-FIG. 6. Experimental characterization of the optimal spiderweb nanomechanical resonator. A False colored scanning electron microscope images of the optimal design. B Qm from the simulation and experimental results of the out-of-plane vibration modes with the figure of the fundamental and optimized mode shapes. C The thermomechanical noise spectra measured at the center of the resonator (orange marker) and at the center of the inner lateral beam (pink marker). The y axis are the normalized power with respect to the maximum power. D Ringdown measurement of the optimized 3 mm spiderweb resonator excited in its 133.6 kHz mode with an extracted quality factor higher than 1.8 billion, compared with the quality factor measured for the fundamental mode at 81.8 kHz. The y axis is the normalized power with respect to the power at time=0 for each curve. E Comparison of the presented experiment result with state-of-the-art experiment reports [6, 7, 14-21, 23, 62-67]. Marker area corresponds to their aspect-ratio t/L. The dashed line corresponds to the mechanical decoherence constraint in Figure 1 .
ing a billion in room temperature environments. This is achieved via a torsional soft clamping mechanism that avoids radiation losses without using phononic crystals or sub-micron lithographic features. While other stateof-the-art resonators require tethers which are hundreds of nanometers in width, our resonators (with micronsized features) can be reliably fabricated at large scales with photolithography. While high-Q m resonators typically require ∼ 20 − 30nm thicknesses, we design ours with 50 nm thickness to simplify the fabrication. Undoubtedly, by designing our resonators tethers with submicron tethers and thinner geometries, we could further improve Q m at the cost of making the fabrication less accessible for general use. The low dissipation rates of the resonator, with f m /Q m ≈ 75 µHz, also represent an important step towards high-precision sensing applica-tions and room temperature quantum technologies. This includes quantum-limited force microscopy [2], "cavityfree" cooling scheme [79] and quantum control of motion at room temperature [80]. What is fascinating is that the machine learning algorithm independently hones in on torsional vibration mechanisms which are actually used by spiderwebs in nature without any knowledge of how a spiderweb functions detect prey. Notwithstanding, we recognize that this data-driven exploration guided by machine learning is just a first step towards rational design of the next-generation of nanomechanical resonators. The demonstrated approach for realizing high-Q m resonance modes is not restricted to the specific spiderweblike design studied in this work. The design strategy might be applied to a wide range of geometries and design problems involving low-throughput simulations or experiments (the most common scenario in engineering and science). We expect future developments in machine learning and optimization together with novel fabrication techniques to lead to unprecedented nanotechnology within the next decade.

METHODS
Fabrication process: Our nanomechanical resonators are fabricated from 58 nm thick high-stress (1.07 GPa) Si 3 N 4 deposited by low pressure chemical vapor deposition (LPCVD) on a silicon substrate ( Figure 5A). The disposition is carried out in-house, which allows to obtain Si 3 N 4 films of arbitrary thickness in stoichiometric form (3/4 ratio of silicon to nitrogen) leading to the uniform high tensile stress in the film. The pattern in the resonator is first written in a positive tone resist (AR-P 6200) by electron beam lithography ( Figure 5B) to create a mask. To do so, the resist is spin-coated on top of the Si 3 N 4 , baked at 155°C, exposed and developed in pentylacetate. Note that we used e-beam lithography instead of optical lithography due to the high level of control and flexibility for the device geometry which is highly beneficial during the development of a new design. However the minimum features of the resonators are designed to be 1 µm to ensure a easier, large-scale fabrication with optical lithography. The pattern is then transferred into the silicon nitride thin-film layer using an inductively-coupled plasma (ICP) etching based on CHF3 plasma etch ( Figure 5C). Next, the resist is removed with dimethylformamide followed by two cleaning steps with hot piranha solution to remove organic residues and diluted hydrofluoric acid solution to remove surface oxides ( Figure 5D). Last, the Si 3 N 4 layer is released from the silicon substrate with an ICP etching with SF6 at -120°( Figure 5E) [81], performed at high pressure and low DC bias to etch isotropically the silicon substrate. This last step does not require any mask given the high selectivity of the chosen chemical against silicon nitride, avoiding any additional cleaning steps and removing any limitation posed by capillary force or stiction usually encountered in isotropic wet etchings (such as KOH or TMAH). Moreover, it is not constrained by the crystal planes of the silicon substrate, enabling the fabrication of arbitrary shapes and avoiding multiple exposures. The final thickness of the Si 3 N 4 films is expected to be 50 nm.
Mechanical characterization setup: All the measurements presented were performed using a custom balanced homodyne detection interferometer ( Figure 5F). The mechanical displacement is probed with a fiber coupled infrared laser (1550 nm). The power is divided into two arms: the local oscillator (90 %) used as interference reference and the signal arm (10 %) terminated with a lensed fiber. The signal arm and the device are mounted on two separate 3-axis nanopositioners placed perpen-dicular to each other, in order to align the device to the focal plane of the lensed fiber. In this way the signal which comes out from the lensed fiber is focused on the device and its reflection collected back inside the fiber. A piezoelectric plate is connected to the sample holder to actuate the devices mechanically. To reduce the effect of gas damping on the measurements, the lensed fiber and sample stage are placed inside a vacuum chamber. With the aid of a turbomolecular and a diaphragm pump, the system can reach a pressure of 4.0·10 −9 hPa. The sensitivity of the measured signal to phase oscillations is maximal in the linear region of the interference signal. To this end, the phase of the local oscillator signal is controlled with a fiber stretcher driven by a proportional-integralderivative (PID) controller implemented with a FPGA board (RedPitaya 125-14) in order to stabilize the interferometer's low-frequency fluctuations using the signal measured from the balanced photodetector as an error signal for a feedback loop. Thermomechanical noise spectra were acquired with an electronic spectrum analyzer without mechanical excitation applied to the piezoelectric plate. On the contrary, for the ringdown measurements the device was first actuated close to the mechanical resonance frequency of interest with the piezoelectric plate until it reached and excited steady-state. Second, the mechanical actuation was turned off and the decay in time of the measured displacement signal was measured with an electronic spectrum analyzer by setting a resolution bandwidth larger than 5 Hz. The bandwidth needs to be larger than the expected linewidth of the resonator, but small enough to increase the signal-to-noise ratio.

Acknowledgements
The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement Nos. 785219 and 881603 Graphene Flagship. This work has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme (No. 17FUN05 PhotoQuant). This publication is part of the project, Probing the physics of exotic superconductors with microchip Casimir experiments (740.018.020) of the research programme NWO Start-up which is partly financed by the Dutch Research Council (NWO). A.C and R. N. acknowledge valuable support from the Kavli Nanolab Delft, in particular from C. de Boer, and from the Technical Support Staff at PME 3mE Delft, in particular from Gideon Emmaneel and Patrick van Holst. A.C and R. N. would like to thank Minxing Xu and Martin Lee for stimulating discussions and early assistance with fabrication and experiments. R. N. would also like to thank Simon Groeblacher for initial support. D.S., M.A.B and R.N. would like to acknowledge the TU Delft's 3mE Faculty Cohesion grant that enabled to start this project.

Conflict of Interest
The authors declare no conflict of interest. Generally, room temperature resonators with state-of-the-art quality factors have converged towards shared similarities in their geometry, material properties and mechanical modes of interest. To date, they all exhibit some of the highest aspect-ratios in microchip technology with free-standing structures nanometers thick but laterally suspended over several millimeters. As discussed in the main text, increasing the lateral size (L) of the resonator with respect to its thickness in the direction of motion (t) improves Q m for beam and membrane nanomechanical resonators since generally Q m ∝ L/t. Another common attribute of all high Q m resonators is the use of highly stretched silicon nitride (Si 3 N 4 ) to fabricate these high-aspect-ratio nanostructures. A resonators' Q m is defined as 2πW/∆W where W is the total stored energy of the resonator, and ∆W corresponds to the energy dissipated for each vibration cycle. Fabricating resonators from tensile materials enables extremely low mechanical dissipation (∆W ) while increasing the resonator's stored mechanical energy W as a consequence of geometric strain engineering [27]. The large initial stress before the vibration increases the potential energy dominantly, whereas the energy loss for each cycle is determined by the bending of the resonator's vibration mode. It means that the critical design parameters for the nanomechanical resonator's Q m are the strain (stress) distribution after releasing the initial pre-stress and the vibration mode shape. The analytical formulation for a one-dimensional beam's bending f m and the corresponding Q m follows the equations below [23].

Author Contributions
Q 0 is the intrinsic quality factor defined as Q −1 0 = Q −1 volume + Q −1 surface , where Q volume is the bulk material loss of Si 3 N 4 and Q surface is the surface loss, linear to the resonator's thickness. For thin resonators operating at room temperature, we assume Q 0 ≈ 6900 t/100 nm [47]. n, σ, ρ are bending mode number, initial stress and density, respectively. λ is E/12σt/L where E is the Young's modulus. Note that the dominant factor 2λ in eq. S2 happens because of the sharp curvature change near the boundaries around a narrow length L c = √ 2Lλ [23], which reduces the quality factor dominantly. A similar formulation for the membrane structure could be found in [46]. Figure 1 in the main text shows the Q m versus vibration frequency (f m ) graph for a 50 nm thick 3 mm Si 3 N 4 beam considering 1.07 GPa initial stress, when we vary the mode number. As it can be seen in it, the pre-stressed double clamped beam shows a significant improvement of Q m compared to the unstressed beam's intrinsic quality factor Q 0 . Also, the perfect soft clamping assumption, which excludes the clamping loss part in eq. S2, emphasizes the major loss of the quality factor from the sharp curvature change in the boundary, especially at lower order modes. Finally, all of these Si 3 N 4 resonators have relied on mechanical modes which have a highest amplitude of oscillation in the resonator center and where most of the induced losses come from how the mechanical mode connects to the microchip substrate where bending is highest. Given that the dominant channel of mechanical resonator loss comes from bending losses, until now most of the design strategies have typically consisted of engineering the decay of the resonator's mechanical mode to reduce bending in the Si 3 N 4 resonators near the clamping points via hierarchical [16,22] structures or phononic crystal designs [14,15].

DERIVATION OF THE QUALITY FACTOR FOR TWO-DIMENSIONAL STRUCTURES
The formulation here is the generalized form of the uniform membrane equations from Yu and co-workers [46]. Considering the out-of-plane deformation as the displacement field corresponding to the vibration mode shape, the strain can be derived from the generalized Euler-Bernoulli equation. Note that we calculate the strain assuming plane stress, constant deformation at the center of the plate u z , as well as constant thickness. The oscillation of the plate u z (x, y)e i2πfmt induces strain in the plate due to bending and elongation during the vibration. Then, the corresponding variation of stress is derived from generalized Hooke's law with the complex Young's The mechanical work done per oscillation, which corresponds to the energy loss, can be derived as Inserting eq. S3 and eq. S4 into eq. S5 and by ignoring the elongation term which is a few orders of magnitude lower than the other terms for the small deformation regime, eq. S5 can be summarized as The system's stored energy can be obtained by calculating the maximum kinetic energy or the maximum elastic energy. Between the two, in this research, we considered the maximum elastic energy of the stretched structure due to bending and elongation. The work done due to elongation happens to be a few orders of magnitude higher than the bending energy, so we ignore the bending energy here. Therefore, From this result we obtain 2πU/∆U : Leading to eqs. (1)-(3) in the main text, when Q 0 =E/E 2 considering E 1 E 2 .

BAYESIAN OPTIMIZATION
FIG. S1. An overview of Bayesian optimization with a one-dimensional minimization example. A Gaussian process regression model is fitted to the observed data (blue markers) and its acquisition function is maximized to select the design parameter (red markers) for the next iteration. The dotted black line corresponds to the unknown black-box function and the solid blue line corresponds to the regression result based on Gaussian process regression. The regression model is plotted as a posterior mean with shaded areas representing 95 percent confidence interval.
Handling simulation-based optimization of spider web resonators requires the optimizer to use simulations with high computational cost. Bayesian optimization [54,57] is a global optimizer that is expected to avoid many local solutions associated to different modes of vibration. Using an online machine learning approach, this method performs optimization while updating limited information in the design space. If the response surface of Q m is unknown, we can train the model with scarce information, starting with a few initial observations and add additional feature evaluations sequentially. The goal is to optimize and track optimal design parameters using several evaluations of the finite element model of the spider web resonator. As shown in Figure S1, Bayesian optimization uses Gaussian process regression (solid blue line) to approximate the unknown response of the function (dashed line) at each iteration. With the obtained data and the probability distribution taking into account the variance of the unmeasured region, the optimizer uses an acquisition function that determines the design variables for the next iteration. Here we used Expected Improvement (EI) as the most commonly used acquisition function. It estimates where the most considerable improvement over the current best results will be so that both exploitation for local area search as well as the exploration for global optimization could be performed simultaneously. Compared to other optimization algorithms, additional computational resources are required to determine the next optimization point. However, when considering expensive function calls and small or medium design dimensions as in our study, they are negligible compared to the finite element simulation time. In this study, the Python library GPyOpt [58] was used. A more detailed information about the method can be found at [57]. The initial design of experiments affects the convergence to the optimum solution. Since our problem is dealing with global optimization, the curse of dimensionality is a crucial issue to be addressed. Even though we use Bayesian optimization to reduce the total number of design evaluations, the number of initial points (=40) is small compared to the design space (six design parameters). The result in Figure S2 shows the optimization history considering four different sets of randomly selected initial points. The figure shows that changing the initial points affects the convergence speed, as expected. However, the optimized design obtained for all the four cases exhibits the same vibration mechanism described in the main text with a similar range of quality factors.

CONVERSION OF ENERGY LOSS REGARDING THE THICKNESS OF THE WEB-LIKE RESONATOR
The optimization considering thickness as a design parameter has shown that minimizing the thickness was not giving the best quality factor, contrary to what is generally observed for straight strings (eq. S2), or phononic crystal resonators [15]. Figure S3A shows the optimization iteration history with the same conditions of the main text  Figure 4 in the main text. C The full motion of the vibration mode shape and an illustration of the portion of energy loss of a 10 nm thickness spiderweb nanomechanical resonator. Other design parameters refer to Table 1 in the main text. except considering the thickness as an additional design parameter. The optimized result converged at around 60 nm (local optimum), which differs from the general trend of minimizing the resonator thickness. To study more in details this effect, we performed two optimizations constraining the resonator thickness to 30 nm and 50 nm, respectively. Afterward we swept the thickness of the two obtained designs from 10 nm to 100 nm. Figure S3B shows the trend. As usual, the quality factor increases for thinner resonators when the thickness is larger than 50 nm. However, below 40 nm, it starts to decrease sharply. The reason of this effect can be understood from Figure S3C, which shows the same design in the main text Figure 4, but modifying the thickness to 10 nm for comparison. As shown in the figure, most of the energy loss starts to focus on the lateral vibration beam, especially because the torsional motion of the lateral beam becomes significant. Furthermore, unlike the thicker designs, the joint region of the inner lateral beams starts to have a concentration of bending energy which starts to keep a large portion of energy loss because the sharp curvature change occurs, even though the soft clamping mechanism isolates the resonator. Considering this motion, it is expected to improve the quality factor for thinner resonators when considering more challenging fabrication with a width below 1 µm. Note that the research in the main text designed a 50 nm thickness spiderweb nanomechanical resonator considering our fabrication process.

MECHANICAL QUALITY FACTOR OF THE WEB-LIKE RESONATOR REGARDING THE RESONATOR SIZE
The length of the mechanical resonator has been discussed as an essential factor for the f m and Q m . Following the trend of straight beams (eqs. S1 and S2), Q m and 1/f m is linear to L when L t, and Q m becomes quadratic to L if we assume perfect soft clamping. To measure the length effects on the mechanical quality factor and the frequency of the spiderweb nanomechanical resonator, we optimized the quality factor of each design starting from 100 µm to 10 mm length of the resonator. Other design parameters and constraints were considered the same as in the main text. As shown in Figure S4, increasing the resonator's length increased the quality factor quadratically, while the frequency was proportional to the inverse of the length. This result also supports that our new resonator follows the expected trend of the effect of soft clamping. Since soft clamping becomes more crucial for longer resonators, we expect to achieve a Q m around 10 billion with a 1 cm web-like resonator.

FUNDAMENTAL AND CENTER DEFECT MODE SHAPE
FIG. S5. The optimal web-like soft clamping mode for A Design A and B Design B obtained from simulation. The normalized bending loss density around the edge of Design A and Design B shows that the bending loss is concentrated on the boundary and the joints between the lateral and radial beams.
To analyze the effect of the target mode shapes on the resulting Q m , we performed two additional optimizations. The first (Design A) focuses on optimizing the Q m of the fundamental vibration mode, while the second (Design B) aims to optimize the Q m when the out-of-plane deformation at the center of the web is maximized. The optimal Design A, shown in Figure S5A, exhibits a small d while pushing the outer lateral beam near the boundary with a large w 2 . At the same time, the w 1 is nearly maximized with the inner lateral beam around the center of the resonator (see Table S2). The result shows a similar strategy to the trampoline design presented in [20], using strain engineering. It maximizes the stress of the radial directional beam, which primarily vibrates in the fundamental mode by adding the mass around the center and the boundary. This suggests that the machine learning optimizer is capable of finding similar solutions to physics-based design approaches.
The optimal Design B increases the number of radial beams (see Table S2), showing that we could achieve a higher Q m at higher order modes, unlike uniform string [23] or membrane [46] resonators. The optimized result adds the lateral beams near the wave propagation node and couples the wave propagation in the radial direction to the lateral direction by the lateral beams' bending motion, as it can be seen in Figure S5B. This demonstrates that the optimizer finds soft clamping modes without any pre-information, by alternating the sharp curvature change near the boundary and distributing it to several joint deformations. Note that for both of the designs, the outer lateral beams distribute the bending loss concentrated at the boundaries, similarly to what was recently proposed by fractal [22] and hierarchical [16] soft clamping.  Figure S6 shows an enlarged view of the design (black) and the vibrational modes (red) obtained at different iterations. The first two designs at iteration 26 and 27 are obtained during the random search phase and the mode is localized in the outer later beam. The following designs at iterations 41 and 51 result from the Bayesian optimization phase and show a localized mode in the inner lateral beam. The last two designs at iterations 78 and 97 show a clear example of the Bayesian optimization exploring, resulting in vibration modes different from the previous iterations.

SEM FROM TILTED VIEW
The spiderweb nanomechanical resonator, fabricated in high stress silicon nitride, is suspended over the silicon substrate as shown in Figure S7. The release step of the suspended nanomechanical resonator, highlighted in blue, is performed from the top, isotropically etching the silicon underneath. It follows that the fixed boundary is etched creating an overhang, as in can be seen in the inset on the top right of Figure S7.

PHOTOTHERMAL EFFECT ON Q FACTOR
The interaction between the nanomechanical resonator and the optical field can lead to photothermal effects, which in turn can mask the intrinsic damping rate of the oscillators. Those effects are proportional to the laser power. Therefore when measuring extreme high quality factor, it is important to ensure that the measured decay is not affected by the laser power. To verify it, we performed different ringdown measurements of the mechanical resonance at around 133.6 kHz of the spiderweb nanomechanical resonator shown in the main text in Figure 6A, by varying the laser power incident on the resonator from 500 µW to 5 µW. The result in Figure S8 shows a comparable decay rate for all the measurements, suggesting that photothermal effects are negligible. The plotted power on the y axis is normalized with respect to the power at the the beginning of the ringdown for each curve to provide a clear comparison. Note that the decaying envelope of the oscillation at carrier frequency is linear until the measured power equals the nanomechanical resonator's thermal fluctuations level leading to large fluctuations of the measured signal. Moreover strong fluctuations caused by unwanted temperature drifts and mechanical vibrations of the setup can bring the measured signal outside the linear region of the interference signal for a short interval of time, resulting in occasional spikes as visible in the curve acquired with a laser power equal to 500 µW at around 1800 seconds.  Figure 3A in the main text.

FINITE ELEMENT SIMULATION MODEL OF THE SPIDERWEB DESIGN
As mentioned in the paper, the optimized result in this research was obtained from the finite element analysis performed by COMSOL [59]. We used plate elements for two-dimensional structural mechanic simulation dealing with a pre-stressed eigenfrequency analysis, considering the thickness to the width ratio. The plate theory was based on the Mindlin theory with a linear elastic material model. To speed up the simulation, we simulated the partial part of the spiderweb resonator. Figure S9A is an example of our iteration 27. The figure shows that it is modeling a portion of the resonator, which depends on N r . The blue lines in the rotational directional edges considered continuity and antiperiodicity periodic conditions for the simulation, and the red line shows the fixed boundary condition. For each simulation, we considered initial in-plane forces ([N/m]) as 1.07 GPa · 50 nm, which is the product of the initial stress and the thickness. We performed the static analysis as the first step, including geometric nonlinearity because of the significant deformation happening from the static analysis [20]. Note that during the stationary analysis, it is always required to give the periodic condition for the periodic edges even if the antiperiodicity boundary condition is considered for the eigenvalue analysis. After performing the stationary analysis, we performed eigenfrequency simulation, also including geometric nonlinearity. We have set up the maximum number of eigenfrequencies to be a hundred while setting the smallest real part as 100 Hz and the largest real part as 1 MHz. Concerning meshing, we divided the model with small squares around every joint to adjust a more fine mesh to consider the bending loss energy more precisely, as shown in Figure S9B . To avoid numerical errors during the optimization, we meshed the upper half of the structure and copied the mesh domain to the lower half of the structure. After the finite element simulation, we calculated the Q m based on eq. 1 in the main text. After calculating Q m of all the out-of-plane modes, we selected the maximized Q m as the performance of the design.