Multi-Frequency Resonance Behaviour of a Si FractalNEMS Resonator.

Novel Si-based nanosize mechanical resonator has been top-down fabricated. The shape of the resonating body has been numerically derived and consists of seven star-polygons that form a fractal structure. The actual resonator is defined by focused ion-beam implantation on a SOI wafer where its 18 vertices are clamped to nanopillars. The structure is suspended over a 10 μm trench and has width of 12 μm. Its thickness of 0.040 μm is defined by the fabrication process and prescribes Young’s modulus of 76 GPa which is significantly lower than the value of the bulk material. The resonator is excited by the bottom Si-layer and the interferometric characterisation confirms broadband frequency response with quality factors of over 800 for several peaks between 2 MHz and 16 MHz. COMSOL FEM software has been used to vary material properties and residual stress in order to fit the eigenfrequencies of the model with the resonance peaks detected experimentally. Further use of the model shows how the symmetry of the device affects the frequency spectrum. Also, by using the FEM model, the possibility for an electrical read out of the device was tested. The experimental measurements and simulations proved that the device can resonate at many different excitation frequencies allowing multiple operational bands. The size, and the power needed for actuation are comparable with the ones of single beam resonator while the fractal structure allows much larger area for functionalisation.


Introduction
Top down nano-electromechanical structures used as NEMS resonators are extensively investigated due to their application in sensing devices. They have small mass, high quality factor and resonance frequencies. Also nanoresonators ensure precision and ultra-fast response when used as mass sensors or AFM probes [1][2][3]. The material dimensions of the nanostructures are extremely small and they exhibit material properties of dissimilar values compared to their bulk analogues [4]. Moreover, the importance of crystalline orientation can not be neglected [5]. Thus, detectable changes of mechanical and electromagnetic properties are also employed for detection [3,6].
Due to its simple geometric shape and approximate 1D-dimension, nanobeam is a preferable choice for a resonating structure. Nanobeams became functionalised for wide variety of mass detectors of extremely small objects of interest in Physics, Chemistry and Biology [1,[7][8][9][10]. Moreover, there is sufficiently good prediction of nanobeams' resonant frequency and modeshapes by using beam models or finite element method (FEM) softwares such as ANSYS or COMSOL [11,12]. That is why nanobeam resonators are the most developed so far and most of the materials suitable for nanoresonator fabrication had already been tested as a nanobeam resonator [13,14]. A drawback of beam geometry is that coupled mode vibrations easily occur. Then, the individual harmonics that produce the overall oscillation are hard to be differentiated [15,16]. The problem can be resolved by an array of nanobeams where each beam is responsible for a particular frequency [17,18]. This way, each nanobeam is responsible for a particular frequency of interest and coupled modes can be avoided. Due to relatively large circuit elements, most often, these solutions does not compromise the overall size of the devices. However, for the sake of miniaturisation, the size of the multifrequency devices must be further optimised. Twenty years ago, radio frequency devices get miniaturised by the fractal shaped antennas. That opened a scientific field that is still not fully explored. By designing fractal shaped antennas engineers are able to have the functionality of an array of dipole antennas while occupying much smaller volume [19,20]. Also, during the last ten years, fractal designs has been utilised in many novel micro-and nano-devices. More specifically, fractal micromixers have been developed to rapidly stir nano-and micro-liter volume solutions [21,22]. In photonics, fractal shapes are widely investigated where surface plasmon resonance is observed [23,24]. In the field of MEMS RF capacitors, fractal geometries are studied for suppressing residual stress, miniaturisation, excitation by single layer capacitance and optimisation of pull-in instability effects [25][26][27].
Regarding NEMS, to the knowledge of the authors, fractal shape has not been yet utilised as an electrostatically driven mechanical resonator. Hence, the current work is a pioneer investigation of the functionality and behaviour of oscillating mass that has self-similar geometry of nanoscale dimension. Our device consists of seven identical ( 6 2 ) star polygons generated by an Iterated Function System (IFS) that forms Sierpinsky flake fractal of second iteration [28]. In our work we show that as in the case of large scale fractal structures [29] the fractal geometry allows multiple resonant frequencies. The quality factor of the resonant peaks is comparable with the one measured in nanobeams fabricated by the same technique. This result can be associated with several features that are present in our device. Firstly, the truss geometry and superior span-to-volume ratio of fractals that allows the resonating body to span relatively wide area for the given mass while still being a rigid structure. A different truss device with similar properties was also exploited in the work of Heritier [2] where two nanowires have been connected by nanostruds forming a nanoladder for AFM probe application. They report sensitivity of a nanowire detector combined with rigidity of a cantilever beam. Secondly, as our fractal structure is clamped at all its vertices and the fabrication technology introduces high residual stress, we expect to some degree the effects of dissipation dilution and higher Q due to clamp-tapering [30,31]. In the recent studies of 1D and hexagonal phonon crystals (PnC) [32][33][34] these enhancing Q effects are combined with soft clamping which further reveals the importance of the lattice geometry for the amount of dissipated mechanical energy through the anchors.
In our work we report fabrication process that allows nanoscale functional Silicon structures of complicated fractal shapes. Our experimental setup consist of electrostatic excitation in vacuum and an interferometric system that provides optical readout. By variation of the AC-voltage, using lock-in amplifier we systematically explore different resonance phenomena. We experimentally detected multiple resonance peaks in the range of [2 MHz,24 MHz] where linear and nonlinear frequency response behaviour was examined. We report quality factors and low excitation energy that are comparable with nanobeams fabricated by the same technique. The FEM simulations were done by COMSOL allowing for stationary study to define the residual stress. Then, the resonance modes are explored by using eigenvalue together with frequency domain studies. We have used COMSOL to show how the symmetry and stress of the resonator affect the resonant modes. Moreover, we use COMSOL to show that the resonant peaks can be detected not only by optical but also, by electrical measurement setup.
The main advantage of the structure is its self-similar geometry. Our fractal resonator anchored with nanopillars has the ability to resonate at many different frequencies while keeping quality factors high enough for detection purposes. The specific shape also ensures very compact design and high span-to-volume ratio. We want to emphasise our versatile fabrication process that allows functional Si-based devices with complex design ready for investigation. Therefore, we have shown that using this top-down technique one can experiment with sophisticate geometries.
The work is structured in five main sections and supporting material is given. In Section 2 , the details of the fabrication process and design characteristics are given. In Section 3 , the measurement setup is described followed by the experimental results. In Section 4 , the simulation results are presented and alternative approaches are discussed. In Section 5 , complementary and future work is discussed as well as side results that are included in the supplementary material.

Device and Fabrication
The present work is a proof of concept of multifrequency mechanical resonator at nanoscale whose top and side view SEM images are shown in Figure 1a,b. To this end we use the top-down fabrication process based on direct focused ion beam (FIB) implantation which was previously employed for nanobeams [35]. The resonating body of the device consists of seven star-polygons that form a fractal structure. Its all 18 vertices are clamped to nanopillars that are 2 µm long and connect the resonator with the insulating SiO 2 layer. The gap between the hanged structure and the insulator is produced when the 2 µm thick top Si-layer has been etched. The middle SiO 2 -layer has thickness of 2 µm too. The structure is suspended over a 10 µm trench and has width of 12 µm. The nanobeams that form the device are 0.12 µm wide and 0.04 µm thick. The thickness of 0.04 µm is defined by the fabrication process and prescribes Young's modulus of 76 GPa which is significantly lower than the value of the bulk material. For the fabrication of the devices we diced 1 × 1 cm 2 chips from a silicon on insulator wafer with 2 µm thick device layer oriented in the 100 crystallographic direction. Before the fabrication, each dice have been cleaned in an organic solvent (acetone) to remove potential organic contamination and rinsed in isopropanol. The samples have been directly patterned by focused ion beam (FIB) implantation of Ga+ at 30 keV, beam current of 10 pA and dose of 1 × 10 16 per cm 2 in a Cross-Beam system (1560xB from Zeiss, Oberkochen, Germany). The ion beam has been controlled with the nanolithography kit Elphy Quantum (from Raith, Dortmund, Germany) to enable the definition of different geometries [11]. It is remarkable to mention that this resist free lithographic step permits to modify the crystallinity of the irradiated silicon creating an etching mask that can be functional from the electrical point of view with the appropriate treatment [35]. According to SRIM (The Stopping and Range of Ions in Matter) simulations, the implanted and damaged volume in silicon goes up to 40nm in depth when ions goes perpendicular to the surface. This has been concluded from the TEM measurements after the ion implantation [35] and the SEM images after the silicon etching [36].
Working at the above mentioned conditions the thickness of the devices is fixed at 40 nm. Supporting pillars have been defined by ion milling and lateral ion implantation to sustain the fractal resonator, in a similar way as in [37] where a suspended lateral electrode was developed for the gatebility of suspended single charged transistors. Tetramethylammonium hydroxide at 80 degree Celsius and 25% concentration has been prepared to selectively etch the non-implanted volume of silicon. This approach permits to fabricate free suspended mechanical structures in only two steps [38]. The fractal pattern has been generated by an equation based home-written Python code which is able to produce wide variety of polygonal fractals while ensuring that the substructures of the resulted fractal shape do not overlap; see [28]. We decided for a structure that consist of ( 6 2 )-polygons where the self-similarity is up to second iteration [39] in the sense of IFS. Finally, in order to improve the electrical conductivity, we put through the devices under an annealing process at high temperature (up to 1000 • C) with boron in a nitrogen reach atmotmosphere to promote the recrystallization and doping (p-type) of the device [35,38]. This final step serves two important purposes for the electrical excitation of the devices to be achieved. It promotes the crystallisation of the amorphous material and it introduces electrical carriers through their doping.

Experimental Set-Up and Characterisation
The devices have been measured in an in-house heterodyne interferometric set-up which allows the optical characterisation of the resonant displacement of NEMS devices in air and in vacuum environment. This system is composed by a vacuum chamber to place the sample, a 633 nm He-Ne laser, an acousto-optic modulator with 40 MHz of frequency shift, an avalanche photodetector of 1 GHz bandwidth which is connected to a 50 MHz Zurich Instrument lock-in amplifier that performs the read-out of the interference signal produced at the photodetector. Our interferometer is placed over an optical table with active vibration isolation. The estimated maximum resolution for the out-of-plane amplitudes measurement is at about 4 nm.
The layout of the device includes two pads where the electrical current can be feed-in while the bottom Si-layer of the chip has been grounded, [40]. Using these pads, IV-curve measurements of nano-crystalline doped silicon (nc-Si) and amorphous silicon (a-Si) structures were conducted. First, an a-Si structure was measured showing an ohmic behaviour with resistance of 33.33 MΩ. Then, an nc-Si structure with improved conductivity was also measured exhibiting an ohmic behaviour with much lower resistance of 420 kΩ. Taking the geometry of the structure into account the resistivity of the annealed structure was computed by COMSOL to be 0.0024 Ω/m. Consequently, the a-Si structure was not able to be excited due to the poor electrical conductivity while the electrically improved nc-Si structures were tested to be functional. The performed measurements use one of those pads to apply alternating current V AC producing potential difference with the ground in order to excite the structure electrostatically leading to transverse mechanical vibration. The improvement of the excitation due to doping in the Si-nanodevices shows that the annealing step reduces the resistivity of the material [41] allowing higher voltage to be distributed further away from the excitation pad of the device. Thus, the electrostatic force between the suspended structure and the ground layer becomes large enough for the resonator to function.
By using the described set-up, the experimental measurements confirm multiple resonance peaks from 2 MHz up to 24 MHz; see Figure 2. The laser spot of about 12 µm just covers the whole structure and we have aligned its centre with the centre of the structure in order to equally cover the membrane. The usual measurement condition we have used between 2 MHz and 17 MHz was peak-to-peak voltage V pp = 18 V and offset voltage V o f f = 0. However, the magnitude of the response is very sensible to the experimental conditions and ranges above 18 MHz need additional measurements with V o f f = 3 V for the targeted resonant peaks to become clear among the others. The resulted frequency response spectrum has multiple nonlinear resonance peaks where the peaks up to 15 MHz are presented in detail. When measurements at the specific peaks were conducted the softening (negative amplitude dependent frequency shift) and hardening (positive amplitude dependent frequency shift) were detected, hence they were swiped-up and -down for different voltages so the hysteresis jumps to become evident. We have also detected the motion of the structure in air by applying up to 40 V DC bias voltage, see Appendix A.2.

MHz
(a.u.) Figure 2. Experimentally measured frequency response spectrum produced by using electrostatic actuation and optical interferometric detection.

The Peaks at 2.37 MHz and 5.1 MHz
In the frequency range from 2.3 to 5.2 MHz we measured linear and nonlinear responses. The peak at 2.37 MHz, see Figure 3a, is an example for nonlinear hardening, where the hysteresis loop becomes larger when V o f f is applied. For lower excitation rates, linear behaviour was detected and Q factor of 805 was calculated. All measured Q-factors can be seen in Appendix A.3. Interestingly, the peak at 5.1 MHz shown in Figure 3b decreases when V o f f is applied. The peak at 5.1 MHz preserves its linear behaviour for all excitation voltages we tried while spanning the allowed range of our equipment. The Q factor for V pp = 20 V, V o f f = 0 V was estimated to be Q ≈ 912 computed at −3 dB [3]. It should be mentioned that other peaks at higher frequencies also decrease its maximum amplitude when the offset voltage V o f f is applied; see Appendix A.1.  The peak at 6.1 MHz in Figure 4a is an example for nonlinear softening combined with hardening, where there are two hysteresis loops that form an inclined plateau. For lower V pp the peak does not have the hysteresis loop from the right-hand side and looks like nonlinear softening. It has been used as a control peak during the measurements. According to the simulations, this peak corresponds to the odd-(1,1) ellipse membrane mode; see Section 4.1.

MHz
The double peak at 7.1 MHz shown in Figure 4b appears from a single peak when V pp is increased up to V pp = 18 V. Also, when V o f f increases it further splits the peak making the hysteresis loops more evident. Few other peaks have similar behaviour, where increased V o f f divides a plateau-like peak. Barely visible hysteresis loops were observed, hence they are not plotted for clarity. If V o f f is applied, another peak at 9.6 MHz merges with the ones shown here. For V pp = 10 there are two peaks (see the blue curve), one just below 9.8 MHz and another one at 9.83 MHz. When V pp increases to 18V, the double inclined plateau and a hardening hysteresis loop become visible; see the cyan and red curves in Figure 5a.

MHz
A resonance peak appears at 10.1 MHz (Figure 5b) which is twice the frequency of the peak shown in Figure 3b. When the peak-to-peak voltage is increased a hardening reveals; see the cyan curve of Figure 5b.
There is a peak at 11.2 MHz but it was impossible to be detected in a great detail. Thus, we did not include a dedicated figure for the corresponding range of the spectrum. It does not have significant hysteresis loops but sometimes appears as a very noisy triple peak.
3.4. The Peaks at 12.3 MHz and 13.05 MHz.
The peak at 12.3 MHz shown in Figure 6a appears at frequency twice bigger than the frequency of the control peak; see Figure 4a. Moreover, the softening type nonlinearity is preserved and clearly visible for higher V pp . However, no significant hysteresis loop was found, hence the swipe-down is skipped in the figure.
In Figure 6b the peak at about 13.05 MHz is shown. It exhibits interesting phenomena where the increase of V pp shifts the linear peak to higher frequencies. We were not able to experience nonlinear hardening or softening in the range of voltages we operate. This way, the peak's resonance frequency increases with the applied voltage without experiencing nonlinearity. The quality factor for V pp = 20 V (yellow peak) was estimated to be Q ≈ 961. The peak at 13.45 MHz (see Figure 7a) exhibits notable softening type nonlinearity when V pp is increased. The detected hysteresis loop spans over a range of 50 kHz. Figure 7b shows us how at 14.3 MHz when V o f f is increased, the double peak and the corresponding hysteresis loops become more notable; see the cyan and the red peaks in 7b. This double peak appears at twice the frequency of the double peak shown in Figure 4b. Again, if only V pp is increased for V o f f = 0 V, a single peak evolves to an inclined plateau and then if further increased it divides into two separate peaks.
We have initially swiped up and down every peak under the conditions of V o f f = 0, V pp = 10 and V o f f = 0, V pp = 18. Then, we applied many different V pp voltages up to V pp = 20 where we reached the limit of our instrument. In that fashion, we recorded the interesting nonlinear behaviour presented in this section and in Appendix A.1. However, for higher frequencies above 18 MHz we needed to use V o f f in order to obtain better response. This way, we have seen that V o f f influences the peaks dissimilarly from V pp . For example, double peaks widen when V o f f is increased and some peaks get lower when V o f f is increased. Hence, we did more measurements of the peaks below 18 MHz for different V o f f conditions. These experimental findings are interesting and their theoretical explanation would be beyond the scope of the present study. One approach would be if nonlinear approximations using Duffing-type equations are applied to the nonlinear peaks which might help to determine the Young modulus; see the Discussion Section 5. In the following section we find this important material property by using our FEM model.

FEM Simulations
We use COMSOL finite element method (FEM) software in order to define the fractal resonator and to investigate its eigenfrequencies and material properties. Different studies were set so we can distinguish how the geometry of the device and the residual stress affect the frequency spectrum of the device. In the following Section 4.1 we compare the resonant frequencies of two hexagonal-based fractal resonators. One, where sections AD = BE = 12 µm but CF = 10 µm with another one where the sections AD = BE = CF = 12 µm, see Figure 8. Then, in Section 4.2 an initial stress is introduced to the device. As it was previously shown with a nanobeam [11], the stress can alter the order at which the different resonant modes appear while swiping the frequency spectrum. Finally, in Appendix B.2 Frequency Domain study was set to show that resonant peaks affect the resistance of the device, hence they can be detected by electrical measurements.
The extremely thin resonator we are modelling needs careful investigation of how the reduced size influence the material properties. An important size dependent property that we take into account is the Young's modulus of polysilicon which for the bulk material has to be at about 160 GPa. According to [4] a Si-film of 50 nm thickness has Young's modulus of about 70 Gpa, which is significantly smaller compared to the one of bulk Si. In order to define the value of Young's modulus we have modelled a double clamped nanobeam of length 4.16 µm, width 540 nm and thickness 40 nm that was fabricated by the same top-down technique explained above in Section 2. The purpose was to fit the resulted eigenvalues with the measurements shown in [11]. As it is shown from the comparison between the simulated and experimental natural frequencies of Si-nanobeam in Section 4.2, the first four resonant frequencies have very similar values. This was possible by reducing the Young's modulus down to 76 GPa and introducing stress to the structure in order to resemble the effect of post fabrication residual stress. Thus we presume the elastic properties of polysilicon were calibrated well for our fabrication technique.
At this point we have changed the nanobeam geometry and introduced the clamped fractal structure geometry in the COMSOL toolkit. This was achieved by a vertical extrusion of previously computed 2D vector shape of the structure's horizontal section. Then, by the Solid Mechanics Interface we defined the material properties Young's modulus, Poisson ratio and density, and the boundary conditions. In our model we have the resonator that has fixed boundary conditions at its vertices and two block domains that correspond to the gap and the insulation SiO 2 layer. The insulation layer has fixed boundary conditions as well. Also, we define the Electrostatic Interface where the excitation pad and the ground(the lower boundary of the SiO 2 layer) are specified. The next step is to mesh the domains. As the height of our device is relatively small in comparison with its length and width, we first use triangle mesh for the upper boundary of the structure which has been swiped down through the domains. When the mesh was evaluated, we need to define the different solvers. For the eigenfrequency study of the relaxed structure shown in the following subsection we use only the eigenfrequency solver. For the investigation of the natural frequencies of the compressed structure that computes the results of Section 4.2, a previous study that defines the bending must be run. The bending is computed by a Stationary study where the boundary conditions that define the side electrodes are displaced towards each other. Furthermore, DC voltage is applied and then gradually decreased to 0V so the structure is firstly attracted towards the bottom Si-layer and then left at its bent equilibrium. The resulted bent shape has its specific stresses that resemble the postfabrication stress, hence this numerical output is plugged again in the Eigenfrequency study where we compute the eigenfrequencies of the compressed structure. For both, Eigenfrequency and Stationary studies we use the standard solver settings of COMSOL. The eigenvale solver uses MUMPS direct solver while the stationary study uses combination of MUMPS and segregated solver where the mechanics and electrostatics are decoupled

Symmetry Axes and Mode Shapes
The simulations shown in Figure 8 represent the mode shapes of the fabricated device from Figure 1. Due to the length difference between the section CF and AD = BE the resulted geometry is symmetric only about two axes. Namely, − → Ox(FC) and − → Oy where O is the center of the fractal. However, if AD = BE = CF then these three sections will lie on three different axes of symmetry. The effect of such symmetry breaking can be seen in Table 1 where the resulted resonant frequencies are given. In both cases the first five mode shapes correspond to elliptic membrane mode shapes (see Table 1). About membrane's modes notation, see [42]. The difference is that the odd-and even-symmetrical variations of (1,1)-and (2,1)-mode split and appear at different frequencies if AD = BE = CF. Thus, if AD = BE = CF the frequency spectrum has less double resonance peaks than in the case of AD = BE = CF.

Mode Alternation Due to Compression Stress
The fabrication process of the device consists of three steps; see Section 2. Previous study shows, that this fabrication process results in a compression stress due to the Gallium implantation that causes nanobeams to buckle [11]. By using COMSOL Eigenvalue study the compression stress was introduced and the first four natural frequencies of the nanobeam were successfully simulated. It was shown that the mode with three nodes appears before the two nodes mode. In other words, what is known as second mode shape appears at lower frequency than the first mode shape, see Table 2, columns 1, 2 and 3. The needed effect of compression was achieved numerically by displacing the side electrodes toward each other. By using the same technique, an eigenvalue study of compressed fractal structure was defined. In COMSOL, the distance between the resonator's pads was decreased by 4 nm. That way we defined buckling with 82 nm displacement at the center of the structure. However, for the buckling to be achieved numerically some voltage must be applied. Then, the voltage was decreased to zero so the model could reach the stable equilibrium at 0 V where the structure has buckled shape. Few different simulations were run with various voltages to verify that the achieved buckled shape is independent of the applied force. Also, as the fabrication technique is the same as in the case of nanobeam, the values of Young modulus = 76 GPa, Poisson ratio = 0.3 and density = 2328 kg/m 3 are also kept the same as in the successful nanobeam model [40]. Finally, an eigenvalue study was performed where the temperature was set to 296 K. In Table 2 we compare the frequency measurements of the first five modes with the frequencies resulted from the numerical study and the first three of them fit very well. When the natural frequencies of the device before and after compression get compared (see Tables 1 and 2), we detected that the standard order of appearance of the elliptical membrane modes is altered due to the stress of the device. Without compression we have even-(0,1), odd-(1,1), even-(1,1), even-(2,1) and odd-(2,1), and after compression we have even-(1,1), even-(0,1), odd-(1,1), odd-(2,1) and even-(2,1) modes.

Discussion
Our resonator employs a nontrivial shape aiming to show that self-similar structures of very small size can be used as mechanical resonators. The functionality get proved and multiple resonant peaks were detected. We comment on the nonlinearities detected at most of the resonances but we do not try to model them due to uncertainties that are beyond the scope of the present work. Nevertheless, let us discuss the evident nonlinear behaviour of our system. Experimentally, we can categorise the different types of nonlinearities that are evident. Nonlinear effects of softening and hardening, and superharmonics that get excited at twice the resonant frequency are evident. Double peaks that get divided by increased bias voltage are also evident. Our measurement set-up does not have enough resolution to detect the displacement of the individual substructures of the device. Thus, we were unable to get clear answer what caused these nonlinear effects. One option is the frequency of the elliptic-like membrane modes to shift with the amplitude. However, if there is a prominent localisation phenomena [43,44], then some resonances of the substructures can get excited and coupled modes will explain the nonlinear peaks.
Theoretically, by utilising nonlinear Duffing-type equations for modelling nanobeams and 2D-membranes, one can explain spring hardening and softening phenomena with the initial stress and the excitation voltage, see [45,46] and the references therein. Thus, by using the coefficients in front of the nonlinear terms, one can fit a nonlinear resonant curve for a particular mode and obtain some material properties such as Young's modulus [46]. Moreover, doping and crystal orientation can influence the nonlinear response of capacitive doped Si MEMS resonators [45]. Also, superharmonic resonance is a nonlinearity induced by the second power of the harmonic term of the capacitive forcing [45]. And finally, double peaks could be due to symmetry breaking.
The fractal structure of our interest does not allow proper vibration modelling by utilising well-known membrane equations due to its specific geometry that brings the following arguments. As seen in Figure A12 of the Appendix, if our device is stretched along x-axis this would lead to elongation along y-axis too. Thus, our device can be associated to the materials with structural hierarchy (truss structures) which have negative Poisson ratio [47]. At resonance, this property will lead to strain and stress components along y-direction that are of opposite sign if compared with the ones produced by non-fractal 2D-membrane. Another factor that can affect the nonlinear response is the doping and our structure was doped once for the shape to be defined, and once again during the annealing step. As doping affects the resistivity [41] it affects the voltage distribution which variations along the structure lead to electrostatic nonlinearities. Also, doping changes the material's mechanical properties due to induced defects in the crystal lattice, hence, can lead to nonlinear response. Indeed, doping can cause both spring hardening and softening phenomena; as shown in [45]. Taking the mentioned above into account we do not focus on modelling the nonlinear behaviour. Instead, we used the opportunity to define the geometry of the device, its material properties and physical interactions into a multiphysics FEM model. In [12,29,34,43] FEM is successfully employed for the analysis of the properties of a vibrating fractal structures and phononic crystals. Thus, we use COMSOL FEM software to vary the stress in x-direction so the resulted eigenmodes fit with the experimental ones. Moreover, in the Appendix B.2 the piezoresistivity of the device is simulated by our FEM model for the case of non-compressed fractal structure. This way, we show that many applications of our device that would rely on electrical resonant peak detection can be achieved.
The two most important application examples are mass and force sensing which are very convenient for NEMS devices. The broadband spectrum of the device allows different types of multi-frequency experiments to be conducted [48]. Also, the multiple double peaks controlled by V o f f (see Figures 4b, 5a and 7b) can be employed as bandpass filters. Arrays of nanobeams are employed for bio-applications on molecular level where high sensitivity is needed but single nanobeam can not provide enough area for a practical device [6]. Taking this into account, the fractal geometry has a beneficial property of high span-to-volume ratio. This ensures relatively large area for sensitive biological element to be placed without compromising with additional mass. If localisation effects are able to be employed, each substructure can be functionalised differently resulting in a multi-analyte bio-sensor.

Conclusions
The present work is dedicated to the fabrication, measurement and analysis of a NEMS device which uses a fractal nanostructure as a suspended resonating body that is capacitively excited. The fabrication of the device is defined by using Ga FIB-implantation, wet etching and annealing in Ba environment. It has been characterised by an interferometric set-up which allows optical characterization in air and vacuum environment. Then, we optically proved the functionality of all our suspended fractal structures as NEMS mechanical resonators. We have detected broadband spectrum where we specified the quality factors of most of the peaks. We also associated the first five peaks with the resonant modes found by our COMSOL FEM model. Then, by utilising our model, we found how the symmetry of the device can affect the resulted spectrum and how the residual stress can affect the distribution of the mechanical modes. Moreover, with COMSOL we found that due to piezoresistivity, resonant peaks can be detected by measuring the resistance of the device built between the pads. Acknowledgments: I would like to acknowledge Martin Riverola for conducting the resistance measurements of the device.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Extended Measurements
Appendix A. 1

. Resonance Peaks between 16 MHz and 24 MHz
The peak at 16.6 MHz is another example of a response that diminishes when V o f f increases, see Figure A1, panel (a). In panel (b) an initial left-side tilt of the peak can be recognised, however no softening was detected within the range of voltages we can operate.
The peak at 17.4 MHz was not detected with sufficient detail. When V o f f increases, the peak moves towards lower frequencies but hysteresis was not detected, see Figure A2. At 18.6 MHz a peak was detected that exhibits linear behaviour when V o f f is increased, see Figure A3.

Appendix A.2. Measurements in Air and Spectra of Other Devices with Defects
Few resonant peaks were detected in Air. The resulted frequency spectrum is shown in Figure A8. For the lower frequencies, DC-voltage variation from 20 V to 40 V was possible. For the higher ones, the peaks were visible only at 40 V. The smaller Q factor is notable, hence fewer peaks were detected. The examined spectrum between 8 MHz and 9 MHz is not plotted because no resonant peaks were found there.
Aside of the device that has been extensively tested for this study, another two fractal structures were fabricated. Unfortunately, both have some defects, however they are still functional. One of them has contamination deposited over the suspended structure and the other has one of its side pillars collapsed. The resulted spectra of these two structures can be seen in Figure A9.

Appendix A.3. Q Factors
In Table A1 we show the estimated Q-factors of most of the detected resonant peaks. As we were able to detect linear frequency response at lower excitation voltages the signal was postprocessed by MATLAB smoothdata function and the resulted linear peak bandwidths were analysed. In Figures A10  and A11 we show the linear response of the resonances from Sections 3.1 and 3.2 together with the smoothed peaks processed by the MATLAB smoothdata function.   Figure A11. Matlab approximation by 'smoothdata' function with Gaussian filter followed by Q-factor derivation that result in 518 and 544 for the 6.03 MHz and 7.121 MHz peaks respectively. Figure A12. FEM simulation of the displacement of the structure when it is stretched along x-axis.
The colours denote the y-component of the displacement where red is positive and blue negative displacement in micrometers, see the legend. Note that the center of the structure is at x = 0, y = 0.

Appendix B.2. Piezoresistance Simulations
As our experimental set up uses laser-doppler vibrometer to detect the resonant peaks, we have no experimental evidence if the fractal structure can be measured electrically. As this detection possibility is important for the practical NEMS applications we use our COMSOL multiphysics model to justify it. To this end, we use Frequency Domain study with coupled Structural Mechanics, Electrostatic forces and Piezoresistive effect. The structural dynamics coupled with electrostatic forces simulate the displacement when excited electrostatically at different resonance frequencies, namely the ones shown in Table 1 column 3. Also, the Structural Dynamics is coupled with the Piezoresistive effect that computes how modal shapes alter the resistance of the structure. In Figure A13 we show the resistance change at few resonant frequencies by the expression R − ∆R where R = 420,219.19 Ω is the resistance of the structure at rest and ∆ R is the stress dependent resistance. Thus, by only taking the Silicon's piezoresistive coupling coefficient into account, we have confirmed that, the detection of the resonant frequency of our device is possible by measuring the alternation of the resistance or the current while varying the excitation frequency. The relation between the electric field, E, and the current, J, within a piezoresistor is: E = ρJ + ∆ρJ. where ρ is the resistivity and ∆ρ is the induced change in the resistivity, and both are rank 2 tensors. The change in resistance is related to the stress, σ: ∆ρ = Πσ where Π is the piezoresistance tensor (SI units: Ω m/Pa) [49]. Silicon has cubic symmetry, and as a result the Π matrix can be described in the following manner: MHz MHz Figure A13. The first three resonant modes of the device that has no initial compression (see Table 1, column 3) and the corresponding resistance peaks. The resistance changes are represented by the expression R − ∆R where R = 420,219.19 Ω is the resistance of the structure at rest and ∆ R is the stress dependent resistance.