Nanomechanical Measurement of the Brownian Force Noise in a Viscous Liquid
Atakan B. Ari, M. Selim Hanay, Mark R. Paul, Kamil L. Ekinci
NNanomechanical Measurement of the Brownian Force Noise in a Viscous Liquid
Atakan B. Ari, M. Selim Hanay, Mark R. Paul, and Kamil L. Ekinci ∗ Department of Mechanical Engineering, Division of Materials Science and Engineering,and the Photonics Center, Boston University, Boston, Massachusetts 02215, USA Department of Mechanical Engineering, Bilkent University, 06800, Ankara, Turkey Department of Mechanical Engineering, Virginia Tech, Blacksburg, Virginia 24061, USA (Dated: December 17, 2020)We study the spectral properties of the thermal force giving rise to the Brownian motion of acontinuous mechanical system — namely, a nanomechanical beam resonator — in a viscous liquid.To this end, we perform two separate sets of experiments. First, we measure the power spectraldensity (PSD) of the position fluctuations of the resonator around its fundamental mode at itscenter. Then, we measure the frequency-dependent linear response of the resonator, again at itscenter, by driving it with a harmonic force that couples well to the fundamental mode. These twomeasurements allow us to determine the PSD of the Brownian force noise acting on the structurein its fundamental mode. The PSD of the force noise extracted from multiple resonators spanninga broad frequency range displays a “colored spectrum”. Using a single-mode theory, we showthat, around the fundamental resonances of the resonators, the PSD of the force noise followsthe dissipation of a blade oscillating in a viscous liquid — by virtue of the fluctuation-dissipationtheorem.
Brownian motion, the random steps taken by a micron-sized particle in a liquid, is a distinct reality of the mi-croscopic world. The Brownian particle is incessantlybombarded by thermally-agitated liquid molecules fromall sides, with the momentum exchange giving rise to arapidly fluctuating Brownian force. One can find an ap-proximation for the Brownian force by integrating outthe many degrees of freedom of the liquid and writinga Langevin equation for the motion of the particle [1].For a single-degree-of-freedom moving along the z axis,the power spectral density (PSD) of the Brownian forcenoise G F ( ω ) is related to the PSD of the particle’s posi-tion fluctuations as G Z ( ω ) = | ˆ χ ( ω ) | G F ( ω ). Here, ˆ χ ( ω )is the complex linear response function or the force sus-ceptibility of the particle and describes how the particleresponds to a harmonic force at angular frequency ω .The simplest description of the dynamics of the Brown-ian particle comes from the assumption of a “white” PSDfor the Brownian force noise that satisfies the fluctuation-dissipation theorem [2]. This approximation, while ne-glecting all effects of inertia and flow-structure inter-action, captures the long-time diffusive behavior of theBrownian particle [2–4].Brownian motion also sets the limits of mechanical res-onators in physics experiments. Mechanical resonatorswith linear dimensions over many orders of magnitude— from meter-scale mirrors [5–7] all the way down toatomic-scale nanostructures [8, 9] – have been used fordetecting charge [10] and mass [11], and for studyingelectromagnetic fields [12] and quantum mechanics [13].A typical continuous mechanical resonator can be de-scribed as a collection of mechanical modes, with eachmode behaving like a particle bound in a harmonic po-tential, i.e., a harmonic oscillator [14]. With the nor- ∗ Electronic mail: [email protected] mal mode approximation, Brownian motion of a contin-uous mechanical resonator can be easily formulated forsmall dissipation [15, 16], when spectral flatness (i.e., awhite PSD) and modal orthogonality can be assumed forthe Brownian force. For a multi-degree-of-freedom sys-tem with large and spatially varying dissipation, however,the theoretical formulation of Brownian motion is non-trivial [6, 17]. Here, the modes are coupled strongly andmotions in different modes become correlated [15, 18]. Apossibility for finding the characteristics of the thermalforce comes from the fluctuation-dissipation theorem, as-suming that one can determine the dissipation in the sys-tem from separate theory, e.g., fluid dynamics [19].While the Brownian force acting on single-degree-of-freedom particles has been directly measured in liquids[3, 4, 20], the few reports on continuous mechanical res-onators remain in the small dissipation limit [21–23]. Inthe presence of large dissipation, experimental challenges,such as dampened signal levels and lack of reliable mo-tion actuation methods, have so far precluded the di-rect measurement of the Brownian force on continuousmechanical systems. Here, we employ optical and elec-tronic measurement techniques to extract the PSD of theBrownian force noise acting on a nanomechanical res-onator in a viscous liquid. The force noise exerted bythe liquid on the resonator has a “colored” PSD and fol-lows the viscous dissipation of the resonator as dictatedby the fluctuation-dissipation theorem [19]. A single-mode approximation obtained from fluid dynamics [24]captures the observed colored PSD at low frequency butdeviates from the experiment with increasing frequencywhere higher mode contributions and details of our ex-ternal driving approach become appreciable.Our experiments are performed on nanomechanical sil-icon nitride doubly-clamped beam resonators under ten-sion. Figure 1(a) shows a false-colored scanning electronmicroscope (SEM) image of a typical beam that has di- a r X i v : . [ c ond - m a t . o t h e r] D ec
500 nm 10 µm (a) ~ xz , W -29 -28 -27 -26 D i s p . N o i s e PS D G W [ m / H z ] Frequency [MHz] D i s p . A m p li t ude W w [ p m ] Frequency [MHz] (c)(b) | c | [ m / N ] Frequency [MHz]
10 20 30 40 50 601234 S p r i ng C on s t. k [ N / m ] Length l [ m m] v F o r c e F [ p N ] Voltage V [mV] m m50 m m40 m m30 m m20 m m15 m m
62 mV50 mV40 mV25 mV31 mV
FIG. 1. (a) False-colored SEM image of a silicon nitride nanomechanical beam ( l × b × h ≈ µ m ×
950 nm ×
93 nm), showingthe structure (green) and the actuators (yellow). (Inset) Close-up image of an electrothermal actuator near one of the clampsof the beam. This is a u-shaped gold film resistor deposited on top of the beam with a thickness of 100 nm and a width of120 nm. (b) PSDs of the displacement noise of beams with different lengths (15 µ m ≤ l ≤ µ m) measured in air around theirfundamental resonant modes and at their centers ( x = l/ k as a function of lengthfor each beam. (c) Driven displacement amplitudes of the 60- µ m-long beam in water at different drive voltages V . Left insetshows the magnitude of the susceptibility | ˆ χ | as a function of frequency, and the right inset shows the responsivity of the forcetransducer. The line is a fit to F = AV with A = 1 . × − N / V . mensions of l × b × h ≈ µ m ×
950 nm ×
93 nm. There isa 2 µ m gap between the beam and the substrate. Thereare two u-shaped metal (gold) electrodes on each end ofthe beam, through which an AC electric current can bepassed (Fig. 1(a) inset). This causes Ohmic heating cy-cles, which in turn generate thermal bending moments.The result is efficient actuation of nanomechanical os-cillations at exactly twice the frequency of the appliedAC current [25, 26]. Both the driven and Brownian mo-tions of the beams are measured using a path-stabilizedMichelson interferometer that can resolve a displacementof ∼ / √ Hz after proper numerical background sub-traction [27]. We focus on the out-of-plane motions ofthe beams at their centers ( x = l/ W ( t )(Fig. 1(a)). For measurements in liquid, the device chipis immersed in a small bath of liquid. Table I lists the di-mensions and experimentally-determined mechanical pa-rameters of all the devices used in this study. In thefollowing analysis, we use a density of ρ s = 2750 kg/m ,Young’s modulus of E = 300 GPa, and a tension force of S = 7 . µ N for all the beams. All experimental detailsand data are available in the SI [28].We first describe how the spring constants k and effec-tive masses m are obtained for the fundamental mode ofthe resonators from thermal noise measurements in air.Figure 1(b) shows the PSD of the displacement noise (po-sition fluctuations) G W (in units of m / Hz) at x = l/ ω π for each resonator around itsfundamental mode resonance frequency. Since the reso-nances are sharply-peaked, G W can be integrated accu-rately over frequency to obtain the mean-squared fluctu-ation amplitude, h W i , for the fundamental mode. Us-ing the classical equipartition theorem, we determine thespring constants of the resonators as k = k B T / h W i , TABLE I. Experimentally-obtained mechanical properties ofthe measured devices.Device l × b × h k m ω / π ( µ m ) (N/m) (pg) (MHz)60 µ m 60 × . × .
093 0.93 10.82 1.4850 µ m 50 × . × .
093 1.24 9.63 1.8040 µ m 40 × . × .
093 1.42 6.53 2.3530 µ m 30 × . × .
093 1.66 3.33 3.5520 µ m 20 × . × .
093 2.09 1.83 5.3815 µ m 15 × . × .
093 3.42 1.26 8.28 where k B is the Boltzmann constant and T is the tem-perature. Thus, k is the spring constant for the fun-damental mode when measured at x = l/
2. The insetof Fig. 1(b) shows k as a function of beam length. Inair, the frequency of the fundamental mode ω and itseffective mass m are assumed to be very close to theirrespective values in vacuum [29]. Thus, with k and ω in hand, m can be found from m = k /ω where m is the effective mass of the beam in the absence of asurrounding liquid. We discuss how the experimentallymeasured values of k , m and ω relate to the theoreticalpredictions for an Euler-Bernoulli beam under tension inthe SI [28].We now turn to the calibration of the forced responsein water. Under a harmonic force F ( t ) = F sin ωt , wecan write the oscillatory displacement of the beam atits center as W ( t ) = W ω sin ( ωt + ϕ ω ), with W ω and ϕ ω being the frequency-dependent displacement amplitudeand phase, respectively. Linear response theory yields W ω = | ˆ χ ( ω ) | F [30]. Assuming that the fundamentalmode response dominates at low frequency ( ω → W dc = F /k Frequency [MHz] G W [f m / H z ] (a) (b)(d) (f)(e) m m m m m m (c) c | m m ˆ | χ | [ m ² / N ² ] ˆ A m p . - S q . S u sc ep t i b ili t y D i s p . N o i s e PS D w f /2 p w f /2 p m m m m G W w f /2 p w f /2 pw f /2 p w f /2 p w f /2 pw f /2 p FIG. 2. Amplitude-squared susceptibility | ˆ χ | (black) and PSD of the displacement fluctuations G W (blue) for each beamin water; beam length is indicated in the lower left of each panel. Both quantities are measured at the center of the beam x = l/
2. Square symbols are experimental measurements. The continuous lines are theoretical predictions using the single-modedescription. The arrows show the approximate positions of the peaks of the second ( ω f / π ) and third mode ( ω f / π ) in fluid(when in range). (see Eqs. (1) and (2) below). Figure 1(c) shows the drivenresponse of the 60- µ m beam at x = l/ V and sweep-ing the frequency of the voltage. The displacement am-plitude at low frequency, W dc , is determined from eachtrace, and F is found as F ≈ k W dc with k from ther-mal noise measurements. From the measured displace-ment amplitudes at different drive amplitudes, one canobtain the force transducer responsivity, shown in the up-per right inset of Fig. 1(c). More importantly, under theassumption that F is constant [28], one can extract theamplitude of the force susceptibility as | ˆ χ ( ω ) | ≈ W ω /F .The left inset of Fig. 1(c) shows that the response ofthe device at different drive voltages (forces) can be col-lapsed onto | ˆ χ ( ω ) | as measured at its center using theproper force calibration.We show measurements of the driven response andthe Brownian fluctuations for each nanomechanical beamin water in Fig. 2. Each double-logarithmic plot corre-sponds to a separate beam, showing the PSD of the dis-placement noise, G W (blue, fm / Hz), and the amplitude-squared susceptibility, | ˆ χ | (black, m / N ), as a functionof frequency. In all plots, the range shown for both G W and | ˆ χ | are adjusted to span exactly two decades; theranges of frequency shown are different. The approxi-mate positions of the peaks of the second mode (not de-tectable at the center of the beam) and the third mode are marked with arrows when in range. Several impor-tant observations can be made. It is evident that G W and | ˆ χ | show different frequency dependencies and peakpositions. It is precisely these variations that we will useto provide an experimental estimate of the PSD of theBrownian force noise. It is also clear from Fig. 2 that, asthe frequency of the fundamental mode increases for thedifferent devices of decreasing length, the overdampedresponse progressively turns underdamped at higher fre-quencies where the mass loading due to the fluid is re-duced [24].The continuous lines in Fig. 2 are from a theoreti-cal description that treats the beams as single-degree-of-freedom harmonic oscillators in a viscous fluid [24]. Wefirst express the linear response function in the familiargeneral form | ˆ χ ( ω ) | = 1[ k − m f ( ω ) ω ] + ω [ γ f ( ω )] . (1)The modal mass m f in fluid is a function of frequencydue to the mass of fluid that is moving in conjunctionwith m . In addition, the dissipation due to the vis-cous fluid γ f is frequency dependent. Both m f ( ω ) and γ f ( ω ) can be determined from fluid dynamics by approx-imating the beam as a long and slender blade (or cylin-der) oscillating perpendicular to its axis in a manner con-sistent with the fundamental mode amplitude profile ofthe beam [24, 31, 32]. This description yields γ f ( ω ) = m T ω Γ b (Re ω ) and m f ( ω ) = m (1 + T Γ b (Re ω )). The F o r c e N o i s e PS D G F [f N / H z × ] Frequency [MHz] m m m m m m m m (a) (b)(d) (e) m m (f) m m (c) w f /2 p w f /2 pw f /2 p w f /2 p w f /2 pw f /2 p w f /2 p w f /2 p FIG. 3. PSD of the Brownian force G F acting on each beam in water. The line is the prediction of Eq. (4) for an oscillatingblade in a viscous fluid. Arrows show the higher mode peak positions as in Fig. 2. The shading indicates the approximatefrequency region where the fundamental mode is dominant. hydrodynamic function of the blade Γ b is expressedas a function of the frequency-based Reynolds num-ber, Re ω = ωb ν f , where ν f is the kinematic viscosityof the fluid [24, 32]. Γ b is a complex valued function,Γ b (Re ω ) = Γ b (Re ω ) + i Γ b (Re ω ), that is determined asan O (1) correction to the hydrodynamic function of anoscillating cylinder in fluid [32, 33]. The mass loading pa-rameter, T = πρ f b ρ s h , is the ratio of the mass of a cylinderof fluid with diameter b to the mass of the beam, where ρ f and ρ s are fluid and solid densities, respectively. Us-ing these ideas, the amplitude-squared susceptibility canbe expressed as [24] | ˆ χ ( ω ) | = 1[ k − m (1 + T Γ b ) ω ] + ω [ m T ω Γ b ] (2)for the fundamental mode of the beam when measuredat x = l/
2. Comparing Eq. (2) with Eq. (1), one canclearly see how the oscillating blade solution provides theparameters for the single-degree-of-freedom harmonic os-cillator.The PSD of the position fluctuations of the fundamen-tal mode of the beam can be expressed as G W ( ω ) = | ˆ χ ( ω ) | G F ( ω ) . (3)It follows from the fluctuation-dissipation theorem [34,35] that the PSD of the Brownian force noise for an os- cillating blade or cylinder in fluid can be expressed as [19] G F ( ω ) = 4 k B T m T ω Γ b (Re ω ) . (4)In Fig. 2, the black lines use Eq. (2) and the blue linesuse Eqs. (3-4) where k and m are measured from exper-iment; the force is calibrated using k at zero frequency;and Γ b and T are calculated from the dimensions anddensity of the beam and the properties of water. In otherwords, there are no free fit parameters.With the displacement noise PSD and the suscepti-bility experimentally determined, we can find an esti-mate of the PSD of the Brownian force noise exertedon the beams by the surrounding liquid using G F ( ω ) = G W ( ω ) / | ˆ χ ( ω ) | . The symbols in Fig. 3(a)-(f) show theexperimentally-obtained PSDs of the Brownian force (inunits of fN / Hz) in water. The monotonically-increasingcontinuous lines are the theoretical predictions for an os-cillating blade in a viscous fluid given by Eq. (4). Thearrows indicate the approximate peak frequencies, ω f and ω f , of the higher modes as in Fig. 2. The experi-mental data in Fig. 3(a)-(f) increase with frequency for ω . ω f as predicted by the theory of an oscillating bladein fluid. However, the experiment begins to deviate fromtheory around ω f for all resonators (Fig. 3(a)-(f)). Aftermaking a dip, the data in Fig. 3(a) and (b) begin to in-crease again around ω f ; this feature remains out of themeasurement range in Fig. 3(c)-(f). We attribute thesedeviations from the theoretical prediction to the influence R e G ' ' Reynolds Number Re D i m en s i on l e ss D i ss i pa t i on m m 50 m m 40 m m 40 m m IPA 30 m m 20 m m15 m m Re G c '' Re w w G b '' w w FIG. 4. Semi-logarithmic and linear (inset) plots showing allthe data in the dimensionless dissipation form. Theoreticalpredictions are shown for an oscillating cylinder (dashed line)and for an oscillating blade (solid line). of the higher modes of oscillation in the driven responseof the beam [28], in qualitative agreement with [33].The general trends of the Brownian force can be madeclearer by plotting the experimental data nondimension-ally. From Eq. (4), the dimensionless variation of theforce PSD can be expressed as Re ω Γ (Re ω ). Figure 4shows Re ω Γ (Re ω ) for each data trace in Fig. 3 usingsemi-logarithmic (main) and linear (inset) plots. Tomake this plot, we have found experimental Γ valuesfrom Γ = G F / (4 k B T m T ω ) for each beam, with Re ω acting as a nondimensional frequency. Also shown inFig. 4 are the theoretical predictions for an oscillatingcylinder (dashed line) and blade (solid line) in a viscousfluid. In addition to the six data sets in water, we in-clude data taken in isopropyl alcohol (IPA) [28]. IPA,with its higher viscosity and lower density compared towater, allows us to extend our dimensionless parameterspace. The data in Fig. 4 extend over two decades in Re ω and follow the viscous dissipation of a blade (or cylinder)oscillating in the liquid over a range of frequencies.A more accurate description of the Brownian dynamicsof the continuous beams used in the experiments shouldinclude the contributions from the higher modes. If weuse a multimode lumped description, we can express thePSD of the displacement fluctuations as (cf. Eq. (3)) G W ( ω ) = ∞ X n =1 , ,... | ˆ χ n ( ω ) | G F,n ( ω ) , (5)where ˆ χ n ( ω ) represents the susceptibility of the n th mode and the even modes do not contribute since the measure-ment is at x = l/
2. Eq. (5) is the sum of the differentmodal contributions where the modes are assumed to beuncorrelated with one another. For small mode number n , a reasonable assumption is that G F,n ( ω ) is indepen-dent of n to yield G W ( ω ) ≈ G F ( ω ) P | ˆ χ n ( ω ) | , whereit is clear that the important quantity is the sum of thesquares of the susceptibilities of the individual modes.Similarly, the driven response should be described us-ing a multimode approach that accounts for the spatially-varying aspects of electrothermal drive. This analysisyields an approximate expression of the form W ω ≈ ¯ α ( πF ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)X n ψ n ˆ χ n ( ω ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (6)Here, F is the magnitude of the electrothermal force and ψ n accounts for the coupling of the drive to mode n [28].We have used the fact that magnitude of the odd modeshapes at the center of the beam are nearly constant inorder to factor out the constant ¯ α ≈ ( φ n (1 / − forsmall and odd n , where φ n (1 /
2) is the normalized modeshape evaluated at the center of the beam.Equations (5)-(6) provide some insight into the devia-tions of the experiment from theory observed in Figs. 3and 4. For the electrothermal drive applied at the distalends of the beam, the coefficients ψ n are not expectedto be constant and will result in non-trivial contribu-tions from the higher modes. Since, in the single modeapproximation, we estimate G F ( ω ) by dividing the dis-placement noise PSD at ω by the driven response at thesame ω , variations in the driven response due to ψ n resultin deviations from expected behavior for high frequencieswhere the influence of the higher modes are significant.We point out that the effects of ψ n cannot be simply de-convoluted or factored out. The electrothermal drivenresponse of the higher modes are entangled due to thelarge damping in the system, and the driven responsebecomes quite complicated as the frequency increases.We highlight that Eq. (6) contains the square of the sumwhereas Eq. (5) is the sum of the squares. As a result,Eq. (6) would contain complicated contributions due tothe cross terms even if the coefficients ψ n could be madenearly constant.This first direct measurement of the PSD of the Brow-nian force in a liquid employing nanomechanical res-onators is a remarkable manifestation of the fluctuation-dissipation theorem. Even a qualitative explanation ofthe experimental deviation from theory has required con-sideration of subtle aspects of the driven response of acontinuous system. In the near future, a transducer ca-pable of exerting forces at arbitrary positions with highspatial resolution [36] may allow for directly determin-ing | ˆ χ n ( ω ) | for several individual modes. This couldthen be used to extend the frequency range of the typeof measurements described here and would lead to fur-ther physical insights into the Brownian force acting ona continuous nanostructure.ABA and KLE acknowledge support from US NSFthrough Grant Nos. CBET-1604075 and CMMI-2001403. MRP acknowledges support from US NSF Grant No.CMMI-2001559. [1] D. Chandler, Introduction to Modern Statistical Mechan-ics , 1st ed. (Oxford University Press, New York, 1987).[2] R. Kubo, M. Toda, and N. Hashitsume,
StatisticalPhysics II: Nonequilibrium Statistical Mechanics , 2nded., Vol. 31 (Springer, Heidelberg, 1985).[3] T. Franosch, M. Grimm, M. Belushkin, F. M. Mor,G. Foffi, L. Forr´o, and S. Jeney, Resonances arising fromhydrodynamic memory in Brownian motion, Nature ,85 (2011).[4] A. Jannasch, M. Mahamdeh, and E. Sch¨affer, Inertialeffects of a small brownian particle cause a colored powerspectral density of thermal noise, Physical Review Letters , 228301 (2011).[5] P. F. Cohadon, A. Heidmann, and M. Pinard, Cooling ofa mirror by radiation pressure, Physical Review Letters , 3174 (1999).[6] A. Gillespie and F. Raab, Thermally excited vibrationsof the mirrors of laser interferometer gravitational-wavedetectors, Physical Review D , 577 (1995).[7] R. X. Adhikari, Gravitational radiation detection withlaser interferometry, Reviews of Modern Physics , 121(2014).[8] J. S. Bunch, A. M. Van Der Zande, S. S. Verbridge, I. W.Frank, D. M. Tanenbaum, J. M. Parpia, H. G. Craighead,and P. L. McEuen, Electromechanical resonators fromgraphene sheets, Science , 490 (2007).[9] A. W. Barnard, M. Zhang, G. S. Wiederhecker, M. Lip-son, and P. L. McEuen, Real-time vibrations of a carbonnanotube, Nature , 89 (2019).[10] A. N. Cleland and M. L. Roukes, A nanometre-scale me-chanical electrometer, Nature , 160 (1998).[11] K. L. Ekinci, X. M. Huang, and M. L. Roukes, Ultra-sensitive nanoelectromechanical mass detection, AppliedPhysics Letters , 4469 (2004).[12] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt,Cavity optomechanics, Reviews of Modern Physics ,1391 (2014).[13] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bial-czak, M. Lenander, E. Lucero, M. Neeley, D. Sank,H. Wang, M. Weides, J. Wenner, J. M. Martinis, andA. N. Cleland, Quantum ground state and single-phononcontrol of a mechanical resonator, Nature , 697(2010).[14] A. N. Cleland, Foundations of Nanomechanics (SpringerBerlin Heidelberg, Berlin, 2003).[15] P. R. Saulson, Thermal noise in mechanical experiments,Physical Review D , 2437 (1990).[16] A. N. Cleland and M. L. Roukes, Noise processes innanomechanical resonators, Journal of Applied Physics , 2758 (2002).[17] E. Majorana and Y. Ogawa, Mechanical thermal noise incoupled oscillators, Physics Letters, Section A: General,Atomic and Solid State Physics , 162 (1997).[18] C. Schwarz, B. Pigeau, L. Mercier De L´epinay, A. G.Kuhn, D. Kalita, N. Bendiab, L. Marty, V. Bouchiat, andO. Arcizet, Deviation from the Normal Mode Expansionin a Coupled Graphene-Nanomechanical System, Physi- cal Review Applied , 064021 (2016).[19] M. R. Paul and M. C. Cross, Stochastic dynamics ofnanoscale mechanical oscillators immersed in a viscousfluid, Physical Review Letters , 235501 (2004).[20] J. Mo, A. Simha, and M. G. Raizen, Broadband bound-ary effects on Brownian motion, Physical Review E - Sta-tistical, Nonlinear, and Soft Matter Physics , 062106(2015).[21] H. Miao, K. Srinivasan, and V. Aksyuk, A microelec-tromechanically controlled cavity optomechanical sens-ing system, New Journal of Physics , 10.1088/1367-2630/14/7/075015 (2012).[22] J. D. Teufel, T. Donner, M. A. Castellanos-Beltran, J. W.Harlow, and K. W. Lehnert, Nanomechanical motionmeasured with an imprecision below that at the standardquantum limit, Nature Nanotechnology , 820 (2009).[23] C. Doolin, P. H. Kim, B. D. Hauer, A. J. R. MacDonald,and J. P. Davis, Multidimensional optomechanical can-tilevers for high-frequency force sensing, New Journal ofPhysics , 10.1088/1367-2630/16/3/035001 (2014).[24] M. R. Paul, M. T. Clark, and M. C. Cross, The stochasticdynamics of micron and nanoscale elastic cantilevers influid: Fluctuations from dissipation, Nanotechnology ,4502 (2006), 0605035.[25] A. B. Ari, M. Cagatay Karakan, C. Yanik, I. I. Kaya,and M. Selim Hanay, Intermodal Coupling as a Probefor Detecting Nanomechanical Modes, Physical ReviewApplied , 034024 (2018).[26] I. Bargatin, I. Kozinsky, and M. L. Roukes, Efficient elec-trothermal actuation of multiple modes of high-frequencynanoelectromechanical resonators, Applied Physics Let-ters , 1 (2007).[27] V. Kara, Y. I. Sohn, H. Atikian, V. Yakhot, M. Lonˇcar,and K. L. Ekinci, Nanofluidics of Single-Crystal Di-amond Nanomechanical Resonators, Nano Letters ,8070 (2015).[28] See Supplemental Material for additional details on de-vice properties and theoretical derivations.[29] V. Kara, V. Yakhot, and K. L. Ekinci, Generalized Knud-sen Number for Unsteady Fluid Flow, Physical ReviewLetters , 10.1103/PhysRevLett.118.074505 (2017),arXiv:1702.07783.[30] J. Sethna, Statistical Mechanics: Entropy, Order Param-eters, and Complexity (OUP Oxford, Oxford, 2006).[31] L. Rosenhead,
Fluid Motion Memoirs: Laminar Bound-ary Layers , 1st ed. (Oxford University Press, Oxford,1963).[32] J. E. Sader, Frequency response of cantilever beams im-mersed in viscous fluids with applications to the atomicforce microscope, Journal of Applied Physics , 64(1998).[33] M. T. Clark, J. E. Sader, J. P. Cleveland, and M. R.Paul, Spectral properties of microcantilevers in viscousfluid, Physical Review E - Statistical, Nonlinear, and SoftMatter Physics , 046306 (2010).[34] H. B. Callen and T. A. Welton, Irreversibility and gen-eralized noise, Physical Review , 34 (1951). [35] H. B. Callen and R. F. Greene, On a theorem of ir-reversible thermodynamics, Physical Review , 702(1952).[36] A. Sampathkumar, T. W. Murray, and K. L. Ekinci, Pho- tothermal operation of high frequency nanoelectrome-chanical systems, Applied Physics Letters , 223104(2006). upplementary Material for“Nanomechanical Measurement of the Brownian Force Noise in a Viscous Liquid” Atakan B. Ari, M. Selim Hanay, Mark R. Paul, and Kamil L. Ekinci ∗ Department of Mechanical Engineering, Division of Materials Science and Engineering,and the Photonics Center, Boston University, Boston, Massachusetts 02215, USA Department of Mechanical Engineering, Bilkent University, 06800, Ankara, Turkey Department of Mechanical Engineering, Virginia Tech, Blacksburg, Virginia 24061, USA
Contents
I. Fabrication and Device Properties
II. Measurements
III. Theory
IV. Appendix: Additional Data References ∗ Electronic mail: [email protected] a r X i v : . [ c ond - m a t . o t h e r] D ec I. FABRICATION AND DEVICE PROPERTIESA. Material Properties
Our resonators are fabricated out of a 100-nm silicon nitride film that is deposited on top a 550- µ m-thick siliconhandle wafer [1] via low-pressure chemical vapor deposition (LPCVD). The material properties of silicon nitridedepend upon the specifics of the LPCVD process. The reported Young’s modulus and density values are in the range250 GPa ≤ E ≤
350 GPa and 2600 kg / m ≤ ρ s ≤ / m , respectively [2–12]. Furthermore, there is an unknownstress in the silicon nitride film, which results in a tension force S in the suspended beams.In Section I C 2, we show that beam theory and the measured eigenfrequencies of the beams allow us to estimate asingle set of optimal material properties that can be used to describe all of the devices used in our experiments. Wehave found that this approach is more insightful than trying to estimate the specific material properties of each beamindividually. The values we have found from this analysis are E = 300 GPa, ρ s = 2750 kg / m , σ = 86 . S = 7 . µ N, as listed in Tables S1 and S3.
B. Fabrication Process
A flowchart of the fabrication process is given in Fig. S1. Transducers and contact pads are patterned using electronbeam lithography (EBL). A 100-nm-thick gold film with a 5-nm chromium adhesion layer underneath is deposited bythermal evaporation. After lift-off, a second step of EBL is performed to define the beam structures. A 60-nm-thickcopper film is deposited via electron beam evaporation, which serves as a dry etch mask. Inductively coupled plasmaetching is used to etch the silicon nitride layer anisotropically; then, the beam is suspended by an isotropic etch stepthat removes the silicon beneath. After the removal of the copper etch mask, the devices are placed in a chip carrierwith 50-Ω strip lines and wirebonds are made between the contact pads and the striplines. Figure S2 shows SEMimages of a completed chip and one of the resonators, a 15- µ m-long beam. (a) (b) (c) (d))(b) (c) )(d) Silicon Silicon Nitride PMMA Gold Copper (e) (f) (g) (h)
FIG. S1: Flowchart for the fabrication process. (a) PMMA resist is spun on the silicon nitride layer, (b) then transducers arepatterned using EBL, followed by deposition of gold. (c) After lift-off, the transducers and contact pads are completed. (d)A second step of EBL is used to pattern the nanomechanical structures. (e) To mask the beam structures during dry etch, acopper etch mask is deposited. (f) Silicon nitride is etched anisotropically down towards the silicon layer, and then (g) siliconis etched isotropically to suspend the beams. (h) After the etch mask is removed, the fabrication is complete. (a) (b)
500 nm (c)
FIG. S2: (a) A 10 mm ×
10 mm chip that has 16 NEMS devices with lengths 15 µ m ≤ l ≤ µ m. (b) Close-up and tilted SEMimage of a 15- µ m-long beam. The gap between the suspended beam and the substrate is approximately 2 µ m. (c) Magnifiedtop view of a transducer. C. Device Properties
In this section, we provide more details on the device dimensions and mechanical properties and parameters. Thedimensions of the beams reported in Table S1 are determined by imaging them in SEM at high magnification. The in-plane dimensions are found to be close to the design dimensions with small deviations due to fabrication imperfections.The thickness of the silicon nitride beams has more uncertainty and is measured to be 92 . ± . .
1. Theoretical Description of the Natural Frequencies and Mode Shapes for a Beam with Tension
Due to the pre-stressed nature of the silicon nitride layer, the suspended beams are under tension. The tensioncauses the dynamics of the beams to deviate from that of a typical Euler-Bernoulli (EB) beam without tension andapproach that of a string under tension. The dynamics can be described by adding a tension term to the EB beamequation, which can be used to find the eigenfrequencies and mode shapes [13, 14] in the usual manner.We will follow the analytical approach provided by Bokaian [13] to describe the dynamics of a doubly clamped-beamunder tension in the absence of a fluid. The result of this analysis will be the eigenfrequencies (natural frequencies) ω n and mode shapes φ n of the beam. While we only need the mechanical parameters of the fundamental mode for thefluid experiments, we have used the higher mode frequencies to determine the material properties of silicon nitride,as we describe in Section I C 2. TABLE S1: Device parameters and experimentally-obtained mechanical properties of the devices used in this study. The mean-squared fluctuation amplitudes measured at the centers of the beams in air and water, (cid:10) W air (cid:11) and (cid:10) W water (cid:11) , respectively,are listed for comparison.Device l × b × h k m α ω / π U S (cid:10) W air (cid:11) (cid:10) W water (cid:11) ( µ m ) (N/m) (pg) (MHz) ( µ N) (m × − ) (m × − )60 µ m 60 × . × .
093 0.93 10.82 0.476 1.48 720 7.64 4.36 4.4150 µ m 50 × . × .
093 1.24 9.63 0.469 1.80 500 7.64 3.29 3.3340 µ m 40 × . × .
093 1.42 6.53 0.462 2.35 320 7.64 2.86 2.7530 µ m 30 × . × .
093 1.66 3.33 0.451 3.55 180 7.64 2.46 2.5620 µ m 20 × . × .
093 2.09 1.83 0.434 5.38 80 7.64 1.95 1.9415 µ m 15 × . × .
093 3.42 1.26 0.423 8.28 45 7.64 1.19 1.18
Our intention in the following is to provide only the necessary details needed to arrive at the expressions for ω n and φ n , and to clearly establish our notation and conventions. We refer the reader to Ref. [13] for more details. For along and slender doubly-clamped beam that is under tension and undergoing small flexural oscillations the dynamicscan be described by EIl ∂ W ( x ∗ , t ) ∂x ∗ − Sl ∂ W ( x ∗ , t ) ∂x ∗ + µ ∂ W ( x ∗ , t ) ∂t = 0 . (S1)The beam geometry is determined by its length l , width b , and thickness h . A long and slender beam indicates l (cid:29) b, h , and small deflections indicate that the maximum beam deflection is much smaller than the thickness of thebeam. Both of these conditions are satisfied by the nanomechanical devices that we explore here. In our notation, x ∗ = x/l is the dimensionless axial coordinate along the beam where x is the dimensional coordinate. In addition, W ( x ∗ , t ) is the transverse displacement (in the z direction) of the beam at axial position x ∗ and time t ; E is theYoung’s modulus, I = bh /
12 is the area moment of inertia, µ = ρ s bh is the mass per unit length of the beam, and S is the applied axial force. The minus sign on the second term indicates that the beam is under tension. The boundaryconditions for the doubly clamped beam are W (0 , t ) = W (1 , t ) = ∂ W (0 ,t ) ∂x ∗ = ∂ W (1 ,t ) ∂x ∗ = 0.In the main text, we focus upon experimental measurements at a specific axial location on the beam x ∗ . Tosimplify the notation for this case, we define the time varying oscillations at location x ∗ as W ( t ) = W ( x ∗ , t ), wherewe have used x ∗ = 1 /
2. However, in the following we will keep the discussion general and use W ( x ∗ , t ).We next assume that the oscillating solutions are of the form W ( x ∗ , t ) = Y n ( x ∗ ) e iω n t (S2)where n is the mode number, Y n ( x ∗ ) is the n th mode shape, and ω n is the n th natural frequency. Substituting thisexpression into Eq. (S1) yields EIl d Y n ( x ∗ ) dx ∗ − Sl d Y n ( x ∗ ) dx ∗ − µω n Y n ( x ∗ ) = 0 , (S3)whose solution can be expressed as Y n ( x ∗ ) = c ,n sinh( M n x ∗ ) + c ,n cosh( M n x ∗ ) + c ,n sin( N n x ∗ ) + c ,n cos( N n x ∗ ) , (S4)which is in terms of the constants M n = ( U + p U + Ω n ) / and N n = ( − U + p U + Ω n ) / . The amount oftension in the beam is described by the nondimensional tension parameter U , where U = S EI/l (S5)represents the ratio of the axial load on the beam to its rigidity. The nondimensional natural frequency for mode n is Ω n , where Ω n = ω n β/l (S6)and β = ( EI/µ ) / is a constant determined by material properties and linear dimensions.Boundary conditions, Y (0) = Y (1) = dY (0) /dx ∗ = dY (1) /dx ∗ = 0, allow for determining the values of theconstants as c ,n = 1 , (S7) c ,n = M n sin( N n ) − N n sinh( M n ) N n [cosh( M n ) − cos( N n )] , (S8) c ,n = − M n N n , (S9) c ,n = − c ,n . (S10)The characteristic equation relating U and Ω n isΩ n + U sinh (cid:18) U + q U + Ω n (cid:19) / sin (cid:18) − U + q U + Ω n (cid:19) / − Ω n cosh (cid:18) U + q U + Ω n (cid:19) / cos (cid:18) − U + q U + Ω n (cid:19) / = 0 . (S11) Beam Position x/l N o r m a li z ed M ode S hape FIG. S3: The first five normalized mode shapes φ n as a function of x/l for a doubly-clamped beam under tension. For thisfigure we have used U = 320 which corresponds to the parameters of the 40- µ m-long beam (Table S1 and S3). Given a value of U , all of the natural frequencies Ω n are found as solutions of Eq. (S11). In practice, we determinethe values of Ω n numerically as the successive roots of Eq. (S11). We will find it convenient to use the normalizedmode shapes φ n ( x ∗ ) = Y n ( x ∗ ) / ˜ N / n where ˜ N n is chosen such that Z φ n ( x ∗ ) φ m ( x ∗ ) dx ∗ = δ nm (S12)with δ nm = 1 for n = m and δ nm = 0 for n = m . Note that the orthonormal eigenfunctions φ n are dimensionless.
2. Determining Material Properties from Natural Frequencies
Since we do not know the precise values of S , E and ρ s , we use an error minimization approach to determine thesefrom theory and experimental observables of the beams in air, namely the eigenfrequencies. Theoretical predictionsof the eigenfrequencies in terms of Ω n and U are given by Eq. (S11) where Ω n depends on ρ s (Eq. (S6)) and U depends on S and E (Eq. (S5)). We take the following approach. We select trial values for E and ρ s . Using the lineardimensions of the beam, the measured ω , and the trial values of E and ρ s , we compute the tension S in the beam bynumerically solving Eq. (S11). We repeat this for all the fundamental frequencies of all the beams and find average S and U values. Once U is determined for these trial values of E and ρ s , we find all of the remaining eigenfrequenciesas the successive roots of Eq. (S11). Then, we form an error function defined as ε ( E, ρ s ) = X k | ω ( e ) k − ω ( t ) k | | ω ( e ) k | , (S13)where k runs through all (available) modes of all the beams, and ω ( e ) k and ω ( t ) k are the experimental and theoreticaleigenfrequencies. We repeat this process and calculate ε by sweeping the ρ s - E plane within the range of acceptablevalues from the literature which are 2600 kg / m ≤ ρ s ≤ / m and 250 GPa ≤ E ≤
350 GPa [2–12]. Wefind that E = 300 GPa and ρ s = 2750 kg / m values provide a global minimum for the function ε ( E, ρ s ) for theaverage tension value S ≈ . µ N. This tension value corresponds to a stress of σ ≈ . E , ρ s and the average value of S .The tension in the beams manifests itself in several ways. First, tension increases the natural frequencies of a beamresonator. Second, the relation between the resonance frequencies deviate from that of EB beam. Figure S4(a) showsthe ratio of ω n /ω for our beams as well as the two limiting cases. Orange symbols correspond to an EB beam where F r equen cy R a t i o ω n / ω F und . F r equen cy f [ M H z ] E ff e c t i v e M a ss m [ pg ] Length l [ m m] String EB Beam 60 um 50 um 40 um 20 um
10 20 30 40 50 602468 10 20 30 40 50 60024681012
Length l [ m m]Mode number n (a) (b) (c) FIG. S4: (a) Ratio of the n th mode frequency to the fundamental mode frequency, ω n /ω for four of the devices used in theexperiments. Theoretical ratios for an Euler-Bernoulli beam without tension (orange) and a string under tension (black) arealso shown. (b) Fundamental mode frequencies f = ω / π and (c) effective masses m for all the devices as a function of theirlengths. ω n /ω ≈ , . , . , . . . for n = 1 , , , . . . . The black symbols correspond to a string under tension where ω n /ω = n .The ratios ω n /ω for all of our experimentally measured beams lie in between these two limits and are shown by thedifferent open symbols. This behavior is due to the fact that U depends quadratically on the beam length as shownin Eq. (S5).Finally, the modeshapes also depend on S , but only weakly. The first five normalized mode shapes calculatedfor a 40- µ m-long beam are shown in Fig. S3. We note that it is experimentally not possible for us to measure themodeshapes with enough precision such that we can distinguish between an EB beam and a string under tension. Allof the relevant parameters used in finding the mode shapes are also given in Table S3.
3. The Spring Constant of the Fundamental Eigenmode
In the main text, we describe how the spring constant of the fundamental mode k for each beam is determined usingthe classical equipartition theorem. Figure S5(a) shows an example of how we calculate the mean-squared thermaldisplacement for the fundamental mode of a 60- µ m device from the power spectal density (PSD) of its Browniannoise in air. Since the thermal resonance is sharply peaked, we can find the area under the curve by integration,which provides the mean-squared displacement h W air i . In Fig. S5(a), this area is the shaded region under the curve(red). Then, using the classical equipartition theorem, we find the spring constant as k = k B T h W air i , where k B is theBoltzmann constant and T is the temperature. The variation of the measured spring constants k with length l forour devices are shown in Fig S5(b) by the data points. In addition, the values for k are listed in Table S1.As a comparison, the static spring constant of a doubly-clamped beam without tension when a point load is applied TABLE S2: Eigenfrequencies (natural frequencies) of the fundamental and higher modes of the devices used in this study.Values in parentheses are the theoretical estimations using the optimal material properties.Device ω / π ω / π ω / π ω / π ω / π (MHz) (MHz) (MHz) (MHz) (MHz)60 µ m 1.48 (1.56) 2.93 (3.16) 4.46 (4.82) 6.10 (6.57) 7.82 (8.44)50 µ m 1.80 (1.90) 3.78 (3.86) 5.78 (5.92) 7.89 (8.15) 10.14 (10.56)40 µ m 2.35 (2.42) 4.66 (4.96) 7.14 (7.70) 9.94 (10.74) 13.04 (14.13)30 µ m 3.55 (3.35) - (6.96) - (11.04) - (15.77) - (21.26)20 µ m 5.38 (5.42) 11.68 (11.68) 18.67 (19.42) - (28.94) - (53.97)15 µ m 8.28 (7.85) - (17.61) - (30.39) - (46.67) - (66.64) -29 -28 -27 -26 Frequency [MHz] (a)
10 20 30 40 50 6001234
Stress TheoryElastic TheoryExperiment S p r i ng C on s t an t k [ N / m ] G W [f m / H z ] D i s p . N o i s e PS D (b) Length l [ m m] FIG. S5: (a) Thermal noise peak of a 60- µ m-long nanomechanical beam resonator in air. The shaded area shows the integral ofthe PSD and is used in the equipartition theorem to calculate the spring constant of the device. (b) Experimentally measuredspring constants k as a function of the device length l (symbols). The blue dotted curve is the theoretical prediction in theabsence of tension given by Eq. (S14). The red solid line includes tension (Eq. (S15)). perpendicularly at the center of the beam, x = l/
2, is given by [15] k = 192 EIl . (S14)When the beam is under significant stress, Eq. (S14) becomes less accurate, since the axial load must be taken intoaccount. Lachut et al. [16] formulated a method to calculate the spring constant of a beam under stress, takinginto account the axial load. For a beam with stress of σ s , the ratio of the increase δk in the spring constant to thestress-free spring constant k is given by δkk = 310 (1 − ν p ) σ s Eh (cid:18) lh (cid:19) , (S15)where ν p is the Poisson ratio. The stress σ s is the result of the axial load S , where S = (1 − ν p ) σ s b . The springconstant is then given by k = k + δk . Figure S5(b) includes the predicted spring constants using Eq. (S14) (blue,dotted line) and using Eq. (S15) (red, solid line), using the linear dimensions, material properties, and tension valuesin Table S3 along with ν p = 0 .
4. The Effective Mass of the Fundamental Eigenmode
The effective mass m of the fundamental mode of a beam can be found using m = k /ω . In modal analysis,the effective mass (and the spring constant) can be defined only up to a constant. However, since we experimentallymeasure the resonance frequencies ω and determine the spring constants k using the equipartition theorem, thereis no ambiguity in our values, and we are able to present quantitative measurements for m . Figure S4(c) shows theeffective mass m as a function of the beam length l . In addition, the numerical values of measured m are listed inTable S1.In our experiments we use a point measurement at the optical spot to quantify the beam dynamics. In the following,we describe how this point measurement is connected with our lumped single-degree-of-freedom model. We equatethe kinetic energy of the fundamental mode as measured at position x = x of the beam to the kinetic energy of theentire beam with physical mass m = ρ s lbh that is oscillating at frequency ω as (cf. [15, 17]). This yields12 m ˙ W ( x , t ) = 12 Z l µ ˙ W ( x, t ) dx, (S16)where m is the effective mass of the fundamental mode when measured at position x . If we assume µ = m/l isconstant and describe the beam motion using the fundamental mode shape as W ( x, t ) = a Y ( x ) e iω t where a is anarbitrary constant determining the amplitude of the motion, this expression can be simplified to m m = 1 Y ( x ) l Z l Y ( x ) dx, (S17)which is equivalent to m m = 1 Y ( x ∗ ) Z Y ( x ∗ ) dx ∗ , (S18)where we now use the nondimensional coordinate x ∗ and x ∗ = x /l . Using the normalized dimensionless mode shape φ , this can be written as m m = 1 φ ( x ∗ ) Z φ ( x ∗ ) dx ∗ . (S19)Lastly, the integral is unity because of the orthogonality relationship for φ n , which gives the final result m m = 1 φ ( x ∗ ) . (S20)In the following, we will define α = 1 /φ ( x ∗ ) which yields the following useful relationship between the effectivemass of the fundamental mode when measured at location x ∗ and the actual mass of the beam: m = α m .Since all of our experimental measurements are at the center of the beams, we will always consider α to bedetermined at x ∗ = 1 /
2. For an Euler-Bernoulli doubly-clamped beam without tension, α = 0 . α = 1 /
2. For the devices used in this work, 0 . ≤ α ≤ . α for all our beams are provided in Table S3. These α values are calculated using the fundamental mode shapes of the beam with tension (Section I C 1) using the tensionvalue S = 7 . µ N.We can theoretically predict a value of m using m = α m , which requires the mode shape, beam geometry, andbeam material properties. This theoretical value is also included in Table S3 as the calculated effective mass. Thedifference between the experimentally measured mass and the calculated effective mass increases with beam length andis due to our limited knowledge of the precise device geometry and material properties. For example, the calculatedeffective mass is an estimation based on device dimensions and does not account for fabrication imperfections noradditional factors such as the masses of the metal films and the contributions of the overhangs at the base of theanchors (see Figure S2(b)-(c)). II. MEASUREMENTSA. Experimental Setup and Protocols
1. Optical Interferometer
A homodyne optical interferometer with path stabilization is used for detecting the nanomechanical motion of thebeams both in noise and driven measurements. A diagram of the optical setup can be found in Figure S6. A HeNelaser with a wavelength of λ = 632 . ∼
25 fm / √ Hz in the range 1-50 MHz with 40 µ W incident on the photodetector. In the noise measurements, we canfurther improve the displacement resolution to ∼ / √ Hz by averaging and numerically subtracting the backgroundnoise from the nanomechanical noise — as discussed below in Section II B. The optical spot is ∼
800 nm in diameter,and the typical power incident on a beam during a measurement is 20 µ W.
2. Experimental Protocols
The devices are glued to a chip carrier; the measurements are performed with the chip carrier placed on a piezoelec-tric
XY Z stage. The laser, shown in Fig. S6, is focused on a beam at the measurement location and the
XY Z stageis moved along the x axis for the measurements. The beam length is much larger than the optical spot radius, andthe measurements are essentially at a single point on the beam. For liquid experiments, we fill the chip carrier withthe liquid (deionized water or isopropyl alcohol) and then seal the carrier with a thin glass lid to prevent evaporation.We note that the depth of the liquid bath is much larger than any device dimensions. The liquid inside the cavitydoes not show any visible signs of evaporation even after several days. xz, W xz, W z, W FIG. S6: Optical detection setup used for the experiments. HWP: half wave plate, QWP: quarter wave plate, PBS: polarizedbeam splitter, PD: Photodetector, PID: Proportional-integral-derivative controller. For thermal fluctuation measurements, thedetection PD is connected to a spectrum analyzer; for driven measurements, it is connected to a lock-in amplifier. The secondPD with a PID controller and a movable mirror is used for path stabilization.
B. Noise Measurements
In noise measurements, a spectrum analyzer is used to detect the optical signals on the photodetector. We nu-merically subtract the background signal from the measurement signal to remove the low-frequency laser noise.Figure S7(a) shows an example of this background subtraction process. Brownian motion signals are recorded byfocusing the laser on a location of interest on the nanomechanical beam (green trace in Fig. S7(a)). We average ∼ data traces for each measurement at a 20-kHz resolution bandwidth. The measurements typically take several hourseach. We then repeat the procedure with identical experimental parameters on the anchor of the beam on top of -12 -11 -10 -9 -8 BackgroundSignalExtracted Signal st Mode2 nd Mode3 rd Mode
Frequency [MHz] S i gna l P o w e r [ W a tt] B ea m P o s i t i on x /l D i s p . N o i s e A m p . [ p m ] Frequency [MHz] Beam Position x/l (a) (b) (c)
FIG. S7: (a) Background subtraction on noise data of a 50- µ m beam in water. The background (black) is subtracted fromthe measured nanomechanical noise (green) to find the subtracted noise (red). This subtracted noise is then converted todisplacement noise by calibrating against the wavelength of the laser. Measurement bandwidth here is 20 kHz. (b) A color-map of the noise response of the of 50- µ m beam immersed in water, showing the noise amplitude as a function of frequencyand location of the measurement on the beam. In order to construct the color-map, we have collected noise measurements,such as those in (a), along the beam length at equal steps, starting from one anchor to the other. The data are interpolatedusing Matlab’s surf function. The color-bar corresponds to the root-mean squared displacement amplitude (within the 20 kHzbandwidth) in units of pm. (c) Experimentally measured mode shapes of the first three out-of plane modes for the 50- µ mbeam. These are cross-sections from the color plot along the beam position at frequencies of 0.23 MHz, 0.63 MHz, and 1.41MHz. VS W R Frequency [MHz] VS W R Frequency [MHz] (a) (b)
FIG. S8: (a) VSWR measurements for a 24.4 Ω electrode and a 50 Ω reference resistance in the range of 0 . − ≤
1% change in the dissipated power. the silicon nitride layer to obtain the background signal (black trace in Fig. S7(a)). The assumption is that noise isuncorrelated and therefore noise powers can be added (or subtracted). The subtracted signal is then calibrated usingthe laser wavelength to find the PSD in units of m / Hz.The frequency and position dependent thermal response of a 50- µ m-long nanomechanical beam is shown inFig. S7(b). For this color-map in water, the PSD of the thermal fluctuations of the beam is measured along thebeam length at equal intervals by moving the XY Z stage. The root-mean square (rms) amplitude at each pointis found by multiplying the PSD with the resolution bandwidth and then taking the square root. The color-mapis then generated using MATLAB’s surf function. Figure S7(c) shows the measured mode shapes of the first threethermal modes of the beam in water along the beam length. These are essentially slices along the frequency axis inthe color-map of Fig. S7(b) taken at frequencies of 0.23 MHz, 0.63 MHz, and 1.41 MHz. While noisy, the measuredmodeshapes match the eigenmodes of the undamped beam shown in Fig. S3 closely. We note that the displacementnoise measured in the experiment is reported as rms and always remains positive, in contrast to the theoretical modeshape.
C. Electrical Properties of the Electrothermal Actuators
Electrothermal actuators are used to excite the harmonic response and measure the susceptibility of the nanom-chanical beam resonators. Bargatin et al. [18] showed that the thermal time constant for such an actuator is on theorder of τ th − = 2 π ×
40 MHz. If one drives the actuators at shorter time scales than τ th , the actuation does notwork efficiently. In our experiments the maximum drive frequency is 5 MHz, well below this thermal cut-off frequency.Moreover, the above-mentioned time constant was determined in vacuum where the heat transfer from the actuatoris dominated by conduction through the solid. In our experiments, the devices are immersed in a liquid, where themedium enhances the heat transfer and facilitates more rapid cooling.The electrothermal actuators are designed so that their resistances are close to 50 Ω. This ensures that the actuationpower is dissipated in the resistor, making the actuation process less susceptible to parasitic circuit elements. Inorder to confirm that the actuation power (and hence the force) remains constant over the frequency range of ourexperiments, we have measured the voltage standing wave ratio (VSWR) of the electrothermal actuators. The VSWRis typically expressed as [19] VSWR = 1 + | S | − | S | , (S21)where S is the reflection coefficient and quantifies the power reflected from the load. The actuator electricalresistances throughout the devices are around 20 ± . . − ≈ .
19 to ≈ .
14 over this frequency range, which corresponds to less than a 1% change in the reflected power. Thus,the power dissipated on the transducer is, for all practical purposes, constant.The two-dimensional sheet resistance of the gold transducer layer is estimated to be R s = 0 .
24 Ω / (cid:3) , and theresistance of the u-shaped actuator wire that sits close to the anchor of the beam is calculated as ≈ . N o r m a li z ed A m p li t ude Beam Position x/l F o r c e F [ p N ] Voltage V [mV] N o r m a li z ed A m p li t ude Beam Position x/l
20 30 40 50 60 70 80
String EB Beam String EB Beam (c)(b)
Air Water (a)
FIG. S9: (a) Calculated drive force amplitude F vs. the drive voltage amplitude V for the 60- µ m beam in air (blue, squares)and in water (black, circles). Air calibration is based on the amplitude at the fundamental resonance frequency, whereaswater calibration is based on the amplitude at ω →
0. The dashed blue and solid black lines are fits to F = AV with A = 2 . × − N / V and A = 1 . × − N / V for air and water, respectively. (b) The experimentally measured modeshape of the fundamental mode of a 50- µ m beam in water (black, squares) that has been actuated at both ends of the beam.(c) The experimentally measured mode shape of the fundamental mode of the 20- µ m beam (violet, circles) driven by only oneactuator near x = 0. In panels (b)-(c), the red solid line is the theoretical mode shape of a string with tension and the bluedotted line is the theoretical mode shape of an Euler-Bernoulli beam with tension. top view of the actuator is shown in Fig. S2(c). The parasitic capacitance of the actuators have been estimated fromother measurements of similar chips to be around ∼
100 pF.
D. Driven Measurements
In the driven measurements, a signal generator drives the transducers with a sinusoidal voltage waveform. Thefrequency of the sinusoidal signal is swept while the amplitude is kept constant. To detect the signal, a lock-inamplifier is used in 2 f mode. The power dissipated on a resistor due to an applied oscillating voltage, V cos (cid:0) ω t (cid:1) ,is proportional to V + V cos ω t , where the constants V and ω are the drive voltage and angular frequency,respectively. The first term in the power dissipation expression causes an increase in the average temperature ofthe resistor and the second term causes the temperature to oscillate at twice the drive frequency. The temperatureoscillations result in nanomechanical oscillations, as described in [18].In the main text, we have calibrated the transducer in water using the response of the resonator at low frequency.This data set in water is shown in Fig. S9(a) as well as another data set that is taken on the same resonator in air. Inair, the drive data are calibrated using the amplitude at resonance — taking advantage of the sharp resonance peak(see Fig. 1(b) in the main text). The amplitude of the mode at the resonance frequency ω = ω can be expressed as W ω = F Q /k . Here, the quality factor Q of the peak has been found to be Q ≈
15 from fitting the resonancelineshape to a Lorentzian. With k , Q and W ω available from experiments, we determine the lumped force F actingon this mode for each drive voltage V . The results are shown in Fig. S9(a).Both force calibration data sets in Fig. S9(a) depend on the square of the applied voltage, as expected. Theefficiency of the transducer is slightly higher in air, i.e., the anchors heat up to a higher temperature in air comparedto water. This is reasonable because the thermal conductivity is higher in water. However, given that that there is nodramatic increase in transducer efficiency going from water to air, we conclude that the thermal conduction in thesestructures mostly takes place through the solid.Even though the structures are heating up to approximately the same temperature in air and water, the oscillationamplitudes in water are significantly smaller than those in air because of the viscous damping. We also note thatthe damping is larger for longer beams. In order to alleviate this, the longer beams (40, 50, and 60 µm ) are drivenusing both transducers on the two ends of the beam as shown in Fig. S2(b) [20]. Figure S9(b) shows the fundamentalmode shape of a 50- µ m-long beam measured when both actuators are driving the nanomechanical motion. Here, themode shape appears symmetric and close to the theoretical predictions. When only one transducer is used, the modeshape becomes slightly distorted toward the end of the beam where the drive transducer is located. An example ofa distorted fundamental mode shape of on a 20- µ m beam is shown in Fig. S9(c). This distortion is due to fact thatthe waves generated by the single transducer decay along the length of the beam because of the high viscosity beforethey can form the familiar standing wave patterns. While this mode shape distortion is not significant for shorter2 G W [f m / H z ] | χ | [ m ² / N ² ] ˆ A m p . - S q . S u sc . D i s p . N o i s e PS D F o r . N o i . PS D G F [f N / H z × ] Frequency [MHz] Frequency [MHz] (a) (b) | c | m m IPA ˆ G W FIG. S10: (a) Amplitude-squared susceptibility | ˆ χ ( ω ) | (black) and PSD of the Brownian fluctuations (blue) for a 40- µ m beamin IPA. Both quantities are measured at the center of the beam. The symbols show the experimental data, and the continuouslines are predictions of the single-mode theory with no free parameters. (b) PSD of the Brownian force acting on the beam inIPA. The continuous lines show the normalized experimental thermal (blue) and driven (black) responses for each beam from(a) as well as the theoretical PSD of the Brownian force (black). beams, we anticipate it to be a source of error when comparing the theory and experimental results. We discuss thisand other sources of errors further in Section III. E. Measurements in Isopropyl Alcohol (IPA)
In addition to water measurements, we have measured the Brownian force PSD for a 40- µ m beam immersed inisopropyl alcohol (IPA) as shown in Fig. S10(a). IPA is more viscous than water ( ν f = 2 . × − m / s) but hasa lower density ( ρ f = 786 kg / m ). Thus, the mechanical response of a given resonator is more damped in IPA ascompared to its response in water. This can be seen by comparing Fig. S10(a) with Fig. 2(c) in the main text. Theexperimental and theoretical curves for the PSD of the Brownian force G F are shown in Fig. S10(b). The theoreticalestimations are found using Eqs. (S45) and (S46) which are discussed further below. All of the necessary fit parametersare included in Table S3. F. Mass Loading and Frequency Shift
In the fluid-structure interaction problem discussed below in Section III B, one obtains an in-phase solution thatcouples to the mass term and gives mass loading. This added mass depends upon the oscillation frequency anddecreases the frequency of the peak in thermal and driven measurements in fluids. In Fig. S11, we show our mea-surements of the peak frequencies of different beams and modes in water, all obtained from driven responses. Thefrequency of each peak in water, ω nf , is normalized by the frequency of the same peak in air, ω n . The data arepresented as a function of the frequency-dependent Reynolds number Re ω , where Re ω is based on ω nf (see Eq. (S37)below). Because all these beams have the same width, Re ω is simply proportional to ω nf . In the region 0 < Re ω ≤ ω nf accurately. III. THEORY
In this section, we briefly derive the expressions that describe the dynamics of the beam with tension in a viscousfluid. We provide our results in terms of a lumped single-degree-of-freedom description. The theoretical descriptiondiscussed in detail in Refs. [17, 21] can be directly applied to our case with only a few modifications. In essence, thetype of beam and the presence of tension affect the mode shape and the frequency of oscillation. Once the mode shapeis known, it is possible to determine the fundamental natural frequency ω and the parameter α which can be usedto determine the effective mass m and spring constant k . With ω , m and k determined, the previously-derivedexpressions [17, 21] can be used to obtain the PSD of the displacement fluctuations of the beam and the susceptibility.3 w n f / w n F r equen cy R a t i o Re Reynolds Number w FIG. S11: Ratio of the n th mode frequency in water ( ω nf ) to that in air ( ω n ) for each beam as a function of frequency-dependentReynolds number, where Re ω = ω nf b ν f . A. A Theoretical Description of the Fundamental Mode
A point measurement made at a specific location on the beam can be described as a lumped single-degree-of-freedom system. In our work, we have made measurements in the middle ( x = l/
2) of the beam. In this description,the nanomechanical beam in its fundamental mode probed at its center can be described as a particle of effective mass m with modal displacement (or amplitude) W ( t ) that is attached to a linear spring of stiffness k . Furthermore, thebeam is exposed to an externally-applied driving force F d ( t ) and to a force due to the surrounding viscous fluid F f ( t ).This can be expressed as m ¨ W + k W = F f ( t ) + F d ( t ) , (S22)where k = m ω . In order to compute the susceptibility ˆ χ ( ω ), we use an unit impulse for the drive force F d ( t ) = δ ( t )and transform into Fourier space to obtain − m ω ˆ W ( ω ) + k ˆ W ( ω ) = ˆ F f ( ω ) + 1 , (S23)where we have used the conventions W ( t ) = 12 π Z ∞−∞ ˆ W ( ω ) e − iωt dω, (S24)ˆ W ( ω ) = Z ∞−∞ W ( t ) e iωt dt. (S25)If the fluid force is separated into its real and imaginary components, these will contribute to the mass and dampingterms, respectively, which we can represent as − m f ( ω ) ω ˆ W ( ω ) − iωγ f ( ω ) ˆ W ( ω ) + k ˆ W ( ω ) = 1 , (S26)where m f ( ω ) is the frequency-dependent mass that includes the contribution from the fluid and γ f ( ω ) is the frequency-dependent damping due to the fluid. Specific expressions for the oscillating cylinder and oscillating blade approxima-tions are discussed in Section III B. Solving for ˆ W ( ω ) and noting that this is the Fourier transform of the susceptibility,we have the desired result ˆ χ ( ω ) = 1 k − m f ( ω ) ω − iωγ f ( ω ) . (S27)4The squared magnitude of the susceptibility is then | ˆ χ ( ω ) | = 1[ k − m f ( ω ) ω ] + [ ωγ f ( ω )] . (S28)The PSD of the beam fluctuations when driven by Brownian motion is directly related to the susceptibility by [17] G W ( ω ) = 4 k B Tω ˆ χ ( ω ) , (S29)where ˆ χ ( ω ) is the imaginary component of ˆ χ ( ω ). This yields G W ( ω ) = 4 k B T γ f ( − m f ω + k ) + ( ωγ f ) . (S30)By the the fluctuation-dissipation theorem, the PSD of the force noise (Brownian force) driving the beam can beexpressed as G F ( ω ) = 4 k B T γ f ( ω ) , (S31)which results in G W ( ω ) = G F ( ω )( − m f ω + k ) + ( ωγ f ) . (S32)Lastly, we note that Eq. (S32) can also be expressed as G W ( ω ) = | ˆ χ ( ω ) | G F ( ω ) . (S33)Therefore, for the single mode description we have the result G F ( ω ) = G W ( ω ) | ˆ χ ( ω ) | . (S34)In the main text, we use the ratio of the experimentally measured G W ( ω ) and | ˆ χ ( ω ) | to explore the spectral propertiesof G F ( ω ). B. Treating the Beam as an Oscillating Blade in Fluid
In order to find the frequency dependence of the damping and mass from the last section, we look to fluid dynamics,i.e., we seek to evaluate the fluid force acting on the oscillating beam. We start by making the following assumptions:the fluid is incompressible; the beam is long and thin ( l (cid:29) b, h ); and its displacements are small compared to any ofits linear dimensions. With these assumptions, the beam can be approximated as an infinitely long cylinder oscillatingtransverse to its axis with the same displacement amplitude as the beam. The radius of the cylinder is taken as b/ ω ˆ W , and can be writtenas [17, 23] ˆ F f ( ω ) = m cyl,e ω Γ c (Re ω ) ˆ W ( ω ) . (S35)Here, Γ c (Re ω ) is the complex hydrodynamic function for the oscillating cylinder and has the formΓ c (Re ω ) = 1 + 4 iK ( − i √ i Re ω ) √ i Re ω K ( − i √ i Re ω ) , (S36)where K and K are zeroth and first order modified Bessel functions of the second kind. The argument of Γ c isexpressed in terms of the frequency-dependent Reynolds numberRe ω = ωb ν f , (S37)5 Blade
Cylinder Re w Reynolds Number R e w G '' D i m en s i on l e ss D i ss i pa t i on FIG. S12: Dimensionless PSD of the Brownian force as a function of the frequency-based Reynolds number Re ω . Eq. (S44) canbe manipulated to show that G F ∝ Re ω Γ (Re ω ). Theoretical predictions for the cylinder (red dashed line) and blade (blacksolid line) descriptions are shown. where ν f is the kinematic viscosity of the fluid and Re ω plays the role of a nondimensional frequency in the theory.The mass term m cyl,e is the effective mass of the oscillating cylinder of fluid of diameter b/ m cyl,e = α n m cyl with m cyl = ρ l π b l the actual mass of the fluid cylinder. The coefficient α n is described in Section I C 4. Since weare focusing on mode 1, we set n = 1 from here on.To account for the rectangular cross-section of the beams, we further apply a correction factor to Γ c (Re ω ) and findthe hydrodynamic function of an oscillating blade Γ b immersed in a liquid [23]Γ b (Re ω ) = ˜Ω(Re ω )Γ c (Re ω ) , (S38)where ˜Ω is a complex valued correction factor that is tabulated in Ref. [23].Separating ˆ F f ( ω ) to its real and imaginary parts, as we did in Eq. (S26), we express the frequency dependent massand damping as [17] m f ( ω ) = m + m cyl,e = m + α ρ f π b l Γ b (Re ω ) (S39)and γ f ( ω ) = α ρ f π b lω Γ b (Re ω ) . (S40)The mass loading parameter T is the ratio of the mass of a cylinder of fluid of diameter b to the actual mass of thebeam which can be expressed as T = m cyl m b = π ρ f bρ s h , (S41)where ρ f and ρ s are fluid and solid densities, respectively. Substituting into Eqs. (S39) and (S40), we write m f ( ω ) = m (1 + T Γ b (Re ω )) (S42)and γ f ( ω ) = m T ω Γ b (Re ω ) . (S43)The PSD of the Brownian force G F then becomes G F ( ω ) = 4 k B T m T ω Γ b (Re ω ) . (S44)6The PSD of the Brownian force is proportional to Re ω Γ b (Re ω ) for the cylinder and blade approximations and isshown in Fig. S12. Substituting the PSDs of the thermal force noise G F , frequency dependent mass m f , and damping γ f into Eq. (S30) and (S28), we arrive at explicit expressions for the PSD of Brownian fluctuations of the beam G W ( ω ) = 4 k B T m T ω Γ b (Re ω )[ k − m (1 + T Γ b (Re ω )) ω ] + ω [ m T ω Γ b (Re ω )] (S45)and for the squared-magnitude of the susceptibility | ˆ χ ( ω ) | = 1[ k − m (1 + T Γ b (Re ω )) ω ] + ω [ m T ω Γ b (Re ω )] . (S46)These expressions can be cast in a slightly different form as G W ( ω ) = 1 k k B T α m b T ω Γ b ( ω )[1 − ˜ ω (1 + T Γ b ( ω ))] + [˜ ω T Γ b ( ω )] . (S47)and | ˆ χ ( ω ) | = 1 k − ˜ ω (1 + T Γ b ( ω ))] + [˜ ω T Γ b ( ω )] (S48)where ˜ ω = ω/ω to yield the expressions provided in [17].In experiment, the beam is driven thermoelastically. The drive can be modeled as a sinusoidally varying forceapplied near the ends of the beam as F d ( x, t ) = ( F l sin( ω d t ) 0 ≤ x ≤ ξ L and ξ R ≤ x ≤
10 otherwise (S49)where ξ L and ξ R indicate the spatial extent over which the drive is applied. In our experiments we have used ξ L = 0 . ξ R = 0 .
95. Following the approach described in Ref. [21] the amplitude response can be expressed as W ω ( ω ) = (cid:18) πF ψ φ ( x ) (cid:19) k − ˜ ω (1 + T Γ b )) + (˜ ω T Γ b ) (S50)where ψ = Z ξ L φ ( x/l ) dx + Z ξ R φ ( x/l ) dx (S51)which captures the coupling of the driving force to the fundamental mode. For the single mode description we canexpress this result as W ω = (cid:18) πF ψ φ ( x /l ) (cid:19) | ˆ χ ( ω ) | , (S52)where it is apparent that the square of the amplitude is proportional to the magnitude of the susceptibility squared. C. Discussion and Possible Sources of Error
The quantitative values used in our study are given in Tables S1, S2, and S3. In Table S3, we have gatheredtogether many of the important parameters used in our study for the purposes of further discussion. The top sectionof Table S3 includes the material properties of the beam and the surrounding fluids that we have used. To recapitulate,the Young’s modulus E , the density of the silicon nitride beam ρ s , and the tension force in the beam S are foundusing the optimization approach described in Section I C 2. This provides a single set of values of ( E, ρ s , S ) that weuse to describe all of the beams in our study. The simplification of using a single set of optimized values to describeall of the beams is expected to result in some error in the parameter values for an individual beam.The middle section of Table S3 includes several important nondimensional parameters. The mass loading parameter T is evaluated using Eq. (S41). The mass loading parameter increases linearly with the density of the surrounding7fluid. This is evident by the increase in T when using water as the surrounding fluid when compared with the T value for IPA. The nondimensional tension parameter U is evaluated using Eq. (S5). Since S is a constant for ourstudy, the value of U increases quadratically with the length of the beam l . We have also included several valuesof the frequency based Reynolds number Re ω given by Eq. (S37). The Reynolds number provides insight into therelative importance of inertial and viscous effects. Each value of the Reynolds number is evaluated at the frequencyof the peak of the fundamental mode ω f of the PSD of the Brownian fluctuations for that fluid. For air we use ω since ω ≈ ω f for this case. In all cases, the Reynolds number remains small where 0 . ≤ Re ω ≤ m , natural frequency ω , and the effective spring constant k of the fundamental mode for the differentbeams. The experimental values of m , ω , and k are directly measured in experiment using air as the surroundingfluid. The frequency of the fundamental peak determines ω and the equipartition of energy is used to determine anexperimental value of k . Using these two measured value of ω and k , the experimental measurement of the effectivemass is determined from m = k /ω .The theoretical values of α , m , ω , and k are found using the theory of an Euler-Bernoulli beam with tensionthat is described in Section I C 1. We use the primed variables, m , ω , and k , to clearly distinguish between thesetheoretical predictions and the experimentally determined values, m , ω , and k . The theoretical value of α is foundusing φ (1 /
2) in Eq. (S20). A useful limiting case is that of a string under tension which yields α = 1 /
2. It is clearfrom our results that as the beam gets longer, the value of U increases, and therefore the value of α approachesa value of 1/2. Once α is determined, it is straight forward to compute the theoretical values m , ω , and k .The theoretical prediction for the effective mass is m = α m b ρ s lbh . The theoretical prediction for the fundamentalfrequency is given by Eq. (S6). Lastly, the prediction for the effective spring constant is k = α ρ s lbhω .It can be seen that there is some disagreement between the experimental and theoretical values of the fundamentalfrequency and a larger disagreement in the values of the effective mass and spring constant. The source of thisdisagreement depends upon several factors. An important contribution is our choice of the averaged material propertiesthat we use to describe all of the beams. For example, it is possible to have lower errors in m , ω , and k if we wereinstead to fit for the parameters for each individual beam while using an adjustable parameter, such as an effectivelength for the beam, and also allowing for Young’s modulus to have a value beyond what is expected for a siliconnitride beam. However, this complicates the description by creating parameter values for each beam and we have notfollowed this approach here. We note that m , ω , and k are not used in generating the theoretical curves shownin Fig. (2)-(3) of the main text, and they are computed here only to provide some insight into the ability of anEuler-Bernoulli beam with tension to capture the dynamics of our experimental beams. In addition, there are severalfactors that we have not included in our analysis. These include the stiffness and the mass of the metallic electrodesas well as the presence of fabrication imperfections such as undercuts.The theoretical predictions of the susceptibility squared | ˆ χ ( ω ) | and the PSD of the Brownian fluctuations G W ( ω )presented in Fig. (2) of the main text require as input m , ω , k , and T . Here, we have used the values for m , ω , and k that were determined from experiments conducted in air for each individual beam. The mass loadingparameter T depends upon the beam width and thickness, the density of the beam, and the density of the fluid. Ofall of these quantities, we are most uncertain of the density of the beam ρ s . We have used the value of ρ s that hasbeen determined using the approach described in Section I C 2.There are several additional possible sources of disagreement between the experimental and theoretical curves of | ˆ χ ( ω ) | and G W ( ω ) near the fundamental mode of oscillation that we would like to highlight. The first is the interactionof the oscillating viscous boundary layer and surrounding solid boundaries [24, 25]. This effect is often called squeezefilm damping. This additional source of damping causes a reduction in the quality factor of the oscillator. The effectof squeeze film damping can be quantified using the oscillating boundary layer thickness δ s that is generated by theoscillating beam (cf. [24, 26, 27]). The boundary layer thickness scales as δ s ∝ p ν f /ω , which captures the dependenceof δ s upon the kinematic viscosity ν f and the frequency of oscillation ω . An estimate shows that δ s ∼ µ m at 100kHz, which is of the same order as the gap distance between the beam and the closest solid surface. An importantaspect is that δ s varies as the inverse square root of the frequency. Therefore, the oscillations of the fundamentalmode will have a larger boundary layer thickness than the higher modes, and in addition, the longer beams will havethe largest boundary layers. For example, the boundary layer thickness decreases to δ s ∼
300 nm at 1 MHz. In termsof our results, the presence of this additional damping would cause the peak region near the fundamental mode forcurves of | ˆ χ ( ω ) | and G W ( ω ) to widen. This effect should be most prominent for our longest beam with l = 60 µ m.From Fig. 2(a) of the main text it is evident that the experimental curves are wider than the theory curves whichwould be expected due to the additional squeeze film damping in the experiment that is not included in the theory.Another source of error is in the mode shapes of the experimental beams. The experimentally measured mode shapesare distorted when compared with the theoretically predicted mode shapes. This is evident for the fundamental modeas shown in Fig. S9(b)-(c). Figure S9(b) shows the 50 µ m beam that has been driven at both ends of the beam. The8 TABLE S3: Parameters used for the calculation of the theoretical G W ( ω ) and | ˆ χ ( ω ) | .Parameter Symbol or Formula Units 60 µ m 50 µ m 40 µ m 30 µ m 20 µ m 15 µ mYoung’s Modulus E GPa 300 300 300 300 300 300Beam Density ρ s kg / m S µ
N 7.64 7.64 7.64 7.64 7.64 7.64Temperature T K 295.15 295.15 295.15 295.15 295.15 295.15Air Kinematic Viscosity ν f m / s × −
156 156 156 156 156 156Water Kinematic Viscosity ν f m / s × − ν f m / s × − - - 27.8 - - -Air Density ρ f kg / m - - 1.18 - - -Water Density ρ f kg / m ρ f kg / m - - 786 - - -Mass Loading Parameter in Water T - 2.91 2.91 2.91 2.91 2.91 2.91Mass Loading Parameter in IPA T - - - 2.29 - - -Tension Parameter U - 720 500 320 180 80 45Reynolds Number in Air Re ω ( ω ) - 0.13 0.16 0.21 0.32 0.49 0.75Reynolds Number in Water Re ω ( ω f ) - 0.34 0.55 0.68 1.36 2.52 3.78Reynolds Number in IPA Re ω ( ω f ) - - - 0.11 - - -Effective Mass (experiment) m pg 10.82 9.63 6.53 3.33 1.83 1.26Fundamental Frequency (experiment) ω π MHz 1.48 1.80 2.35 3.55 5.38 8.28Spring Constant (experiment) k N / m 0.93 1.24 1.42 1.66 2.09 3.42Effective Mass Coefficient (theory) α - 0.476 0.469 0.462 0.451 0.434 0.423Effective Mass (theory) m pg 6.93 5.70 4.49 3.29 2.11 1.54Fundamental Frequency (theory) ω π MHz 1.56 1.90 2.42 3.35 5.42 7.85Spring Constant (theory) k N/m 0.66 0.81 1.04 1.46 2.45 3.75 mode shape is significantly wider than what is predicted theoretically. Figure S9(c) shows the experimental modeshape of the 20 µ m beam which has been driven only at one end of the beam near x = 0. In this case, the fundamentalmode shape is skewed toward the edge of the beam near the electrothermal drive. In this case, it can be seen thatwhen a point measurement is made at the center, the measured amplitude will be ∼
20% smaller than the predictedone due to the distortion, resulting in a less accurate linear response measurement.9
IV. APPENDIX: ADDITIONAL DATA
In this appendix, we present additional data plots using different axes and non-dimensionalizations. All of therelevant information is included in the figure captions.
Frequency [MHz] G W [f m / H z ] | χ | [ m ² / N ² ] ˆ A m p . - S q . S u sc ep t i b ili t y D i s p . N o i s e PS D (a) (b)(d) (f)(e) (c) w f /2 p w f /2 p w f /2 p w f /2 p w f /2 p w f /2 pw f /2 p w f /2 p ˆ | c | m m m m m m m m m m m m G W FIG. S13: Figure 2 from main text plotted in linear scale. Amplitude-squared susceptibility | ˆ χ | (black) and PSD of thedisplacement fluctuations G W (blue) for each beam in water; beam length is indicated in the lower left of each sub-figure. Bothquantities are measured at the center of the beam. The arrows show the approximate positions of the second ( ω f / π ) andthird mode ( ω f / π ) peaks in fluid (when in range). m m m m m m Reduced Frequency ω / ω G W [f m / H z ] | χ | [ m ² / N ² ] ˆ A m p . - S q . S u sc ep t i b ili t y D i s p . N o i s e PS D (a) (b)(d) (f)(e) (c) w f / w w f / w ˆ w f / w w f / w w f / w w f / w w f / w w f / w m m m m G W | c | m m FIG. S14: Figure 2 from main text plotted as a function of reduced frequency, ω/ω , where ω is the fundamental frequencyin air. The arrows show the approximate positions of the second ( ω f /ω ) and third mode ( ω f /ω ) peaks in fluid (when inrange). m m m m m m m m Reduced Frequency ω / ω N o r m a li z ed R e s pon s e | | c ˆ (a) (b)(d) (e) (f) m m m m (c) G W FIG. S15: Normalized amplitude-squared susceptibility | ˆ χ | (black) and PSD of the displacement fluctuations G W (blue) as afunction of reduced frequency, ω/ω , for all the devices in water. m m m m m m Reduced Frequency ω / ω R e G '' D i m en s i on l e ss D i ss i pa t i on (a) (b)(d) (f) m m m m (c)(e) m m w FIG. S16: Semi-logarithmic plots showing the measured dimensionless dissipation for all the devices in water as a function ofreduced frequency, ω/ω . Reduced Frequency ω / ω m m (c) m m m m (f) (e) (d) m m (b) m m m m (a) R e G '' D i m en s i on l e ss D i ss i pa t i on w FIG. S17: Linear plots of the data in Figure S16. [1] University Wafer, 1918 Stoichiometric LPCVD low stress silicon nitride on silicon. , URL .[2] S. D. Senturia,
Microsystem Design (Kluwer Academic, Dordrecht, 2002), ISBN 0-7923-7246-8.[3] C. J. Drummond and T. J. Senden, Materials Science Forum , 107 (1995), ISSN 02555476.[4] S. Santucci, A. V. La Cecilia, A. DiGiacomo, R. A. Phani, and L. Lozzi, Journal of Non-Crystalline Solids , 228 (2001),ISSN 00223093.[5] A. Stoffel, A. Kov´acs, W. Kronast, and B. M¨uller, Journal of Micromechanics and Microengineering , 1 (1996), ISSN09601317.[6] T. Larsen, S. Schmid, L. G. Villanueva, and A. Boisen, ACS Nano , 6188 (2013), ISSN 19360851.[7] M. E. Fauver, D. L. Dunaway, D. H. Lilienfeld, H. G. Craighead, and G. H. Pollack, IEEE Transactions on BiomedicalEngineering , 891 (1998), ISSN 00189294.[8] K. E. Petersen, Proceedings of the IEEE , 420 (1982), ISSN 15582256.[9] R. L. Edwards, G. Coles, and W. N. Sharpe, Experimental Mechanics , 49 (2004), ISSN 00144851.[10] J. Yang, J. Gaspar, and O. Paul, Journal of Microelectromechanical Systems , 1120 (2008), ISSN 10577157.[11] J. Yang and O. Paul, in Sensors and Actuators, A: Physical (Elsevier, 2002), vol. 97-98, pp. 520–526, ISSN 09244247.[12] O. Tabata, K. Kawahata, S. Sugiyama, and I. Igarashi, Sensors and Actuators , 135 (1989), ISSN 02506874.[13] A. Bokaian, Journal of Sound and Vibration , 481 (1990), ISSN 10958568.[14] I. Stachiv, Journal of Applied Physics , 214310 (2014), ISSN 10897550.[15] A. N. Cleland, Foundations of Nanomechanics (Springer Berlin Heidelberg, Berlin, 2003), ISBN 978-3-540-43661-4.[16] M. J. Lachut and J. E. Sader, Physical Review B - Condensed Matter and Materials Physics , 085440 (2012), ISSN10980121.[17] M. R. Paul, M. T. Clark, and M. C. Cross, Nanotechnology , 4502 (2006), ISSN 09574484, 0605035.[18] I. Bargatin, I. Kozinsky, and M. L. Roukes, Applied Physics Letters , 1 (2007), ISSN 00036951.[19] J. Joines, W.T. and Palmer, W.D. and Bernhard, Microwave Transmission Line Circuits (Artech House, Norwood, 2012),ISBN 9781608075690.[20] V. Bello, A. B. Ari, M. Selim Hanay, and K. L. Ekinci, in (IEEE, 2020), pp. 1–5, ISBN 978-1-7281-4460-3.[21] M. T. Clark, J. E. Sader, J. P. Cleveland, and M. R. Paul, Physical Review E - Statistical, Nonlinear, and Soft MatterPhysics , 046306 (2010), ISSN 15393755.[22] L. Rosenhead, Fluid Motion Memoirs: Laminar Boundary Layers (Oxford University Press, Oxford, 1963), 1st ed., ISBN0198563159, 978-0198563150.[23] J. E. Sader, Journal of Applied Physics , 64 (1998), ISSN 00218979.[24] M. T. Clark and M. R. Paul, Journal of Applied Physics , 094910 (2008), ISSN 00218979, URL http://aip.scitation.org/doi/10.1063/1.2912989 .[25] C. Lissandrello, V. Yakhot, and K. L. Ekinci, Physical Review Letters (2012), ISSN 00319007.[26] C. P. Green and J. E. Sader, Physics of Fluids , 1 (2005), ISSN 10706631.[27] R. J. Clarke, O. E. Jensen, J. Billingham, A. P. Pearson, and P. M. Williams, Physical Review Letters96