Olivier Maillet

Xin Zhou

Rasul Gazizulin

Ana Maldonado Cid

Martial Defoort

Olivier Bourgeois

Eddy Collin

Ana Maldonado Cid

& Eddy Collin

Non-linear Frequency Transduction of Nano-mechanical Brownian Motion

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Emerging from the tremendous development of micro/nano technologies, nano-electro-mechanical systems (NEMS) have opened unique capabilities to both engineers and physicists. In the rst place, they serve as ultra-sensitive probes for force sensing [1] with applications e.g. to mass, charge, and even single electronic spin detection [24]. In the second place, these objects are extremely fruitful (weakly) non-linear devices that are able to implement useful functions like e.g. mechanical frequency mixing [5], amplication [6] and bit storage [7].

On the fundamental level, high-quality NEMS structures can be thought of as model systems in which basic phenomena can be advantageously reproduced; one example being the ubiquitous bifurcation mechanism [810].

Ultimately, when coupled to a quantum-limited detection scheme such as a microwave cavity or a Single-Electron Transistor, their sensitivity can be brought to the quantum limit [11,12]. This leads to a unique platform realizing the ultimate force detector foreseen by C. Caves in the 80's [13]. Such moving structures that are macroscopic relative to the atomic scale but follow the laws of quantum mechanics are currently under development for tests of quantum foundations [1416]. Furthermore, they are thought to be a unique new quantum electronics component enabling e.g. coherent photon conversion from the microwave to the optical domain [17,18]. Essentially all applications require in the rst place the resonance frequency of the mechanical mode in use to be as stable as possible. As such, the understanding of the sources of frequency uctuations in nano-mechanical devices becomes an essential technical topic [1,1923].

But in the rst place, it is also a fundamental research goal: the measured frequency noise in actual devices is much larger than all expectations [22,2426], demonstrating even non-linear features for carbon-based systems [27,28]. Thus, attempts have been made to model noise sources [29,30], or to create model experiments experimentally demonstrating the underlying mechanisms [19,31,33,34].

Clever driving schemes taking advantage of nonlinearities have been devised to signicantly suppress frequency noise [35,36]. But what shall be an ideally frequency-noise minimizing nano-mechanical system in the rst place? We know that at lowest order, the dynamics of a mechanical structure can be described by a family of normal modes which are nothing but independent harmonic oscillators. Pushing to the next or-der, these modes are weakly non-linear (so-called Duing resonators) and are dispersively coupled one to the other [3739]. Since all of the modes are unavoidably coupled to a thermal reservoir (ideally the same one), Brownian motion of each of the modes will transduce into a frequency noise on all the others [19,22,29,34,40], and also on itself. Even in a system realized with ideal materials having no internal sources of noise, this built-in mechanism shall x an ultimate limit to the mechanical resonance frequencies stability at T = 0. Only in the limit of T → 0, when all the modes are in their quantum ground state, do the dispersive couplings lead to a simple frequency renormalization of the resonances through the zero-point-uctuations of each of them: a sort of mechanical Lamb shift that dresses all the modes [41].

In the present article, we report on a model experiment in which we use very high quality silicon-nitride NEMS cooled down to Kelvin temperatures. A single mode is driven by a stochastic force, leading to eective temperatures as high as 10 9 K for this mode only. We extract the eect of this articial out-of-equilibrium heating on the mode itself, both by measuring the spectrum of the motion and by measuring the simultaneous response of the same mode to a sine-wave excitation. The eect on a nearby mode is measured with the sine-wave excitation scheme. The setup is carefully calibrated [42], while the devices' characteristics are obtained by both measurements and calculations; the agreement with theory is obtained with no free parameters. Besides, the experiment is performed on dierent devices proving the reproducibility of the results.

We demonstrate experimentally the two regimes of the Brownian motion transduction, named after analog phenomena present in Nuclear Magnetic Resonance (NMR): motional narrowing and inhomogeneous broadening [31].

Based on Ref. [40] and simple expansions of Euler-Bernoulli theory (including non-linear coecients [38,39,[START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]) we give the analytic tools enabling the calculation of the ultimate frequency stability reached by any doubly-clamped device, depending on stress, dimensions and temperature T [45]. For bottom-up structures like e.g. carbon nanotubes with high aspect ratio, this limit is not marginal [29].

2 Results

2.1

The nano-electro-mechanical systems

The devices under study are doubly-clamped siliconnitride nano-beams having width w =300 nm and thick- Ref. [42]. We can thus infer forces F n and displacements

x n in S.I. units, and compute the devices characteristics (namely mass m n , spring constant k n , non-linear coecients β n,m ). These match the expected calculated values; note that a particular care has been taken in the calibration of the noise source. The only t parameter is indeed an overall correction of the force noise not exceeding 15 % in amplitude (same order as in Ref. [10]).

Actuation and detection are performed with the mag-netomotive scheme [42,46]. A drive current (composed of both the Gaussian noise component centered around resonance frequency ω 1 and a sine-wave of frequency ω close to ω n , with n = 1 or n = 3) is injected in the NEMS metallic layer through a home-made adder and a 1 kΩ bias resistor. In an in-plane d.c. magnetic eld orthogonal to the beams, this generates an out-of-plane driving force F n (t) with harmonic component F 0 n cos(ωt).

The motion is detected through the induced voltage by means of a standard lock-in detection. We obtain the two quadratures, in-phase (X) and out-of-phase (Y) with respect to the local oscillator.

In order to preserve our calibration capabilities, the lock-in has also been used for the spectral measurements S n X (ω) of the Brownian motion of mode n = 1. Moreover, this enables to measure uctuations on each of the two quadratures X, Y independently (plus their crosscorrelations). When the sinusoidal excitation is weak (or nonexistent), the spectra on X and Y are equivalent and no correlations are detected; this is the range of validity of the work presented here. However, sig-natures of squeezed statistics of motion [43] can be observed on measured spectra when the sinusoidal excitation is too large. Details on the measurement technique, calibrations and calculated parameters can be found in S.M. [45].

Dispersive coupling driven by stochastic motion

Linear motion of thin nano-beams is very well described by the Euler-Bernoulli equation [START_REF] Cleland | Foundations of nanomechanics[END_REF]. The basic ingredients involved are the inertia (through the density ρ beam ), the Young's modulus E beam and the tension T 0 generated by the in-built stress. For doubly-clamped beams, the non-linear behaviour is well understood: it arises from the stretching of the device under transverse motion x [38,39,[START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]. This geometric non-linearity results in a tensioning T 0 + δT of the beam with δT ∝ x 2 , which can be incorporated into the beam equation [START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]. This leads to a frequency shift of the modes that is proportional to the square of the displacement. When only two modes n, m are under study, it writes:

ω n = ω 0 n + β n,n x 2 n + β n,m x 2 m , (1)
ω m = ω 0 m + β m,m x 2 m + β m,n x 2 n , (2)
where we introduced ω 0 n , ω 0 m the linear resonance fre- quencies, and the Dung non-linear coecients β i,j [39].

We remind for the interested reader the mathematical derivation of these expressions in S.M. [45].

Eqs. (1-2) can be adapted when one of the motions, say x n , is a stochastic variable: (b) In-phase (X) and quadrature (Y) components measured for mode n = 1 while driving noise on the same n = 1 mode for sample 300 µm-n • 2 (intra-mode). (c) Same measurement performed on mode n = 3 while driving uctuations on n = 1 for sample 300 µm-n • 1 (inter-mode coupling). The standard deviation ∆x 2 1 (i.e Brownian motion level) is increased from left to right (essentially from MN to IB regime, see Fig. 3), and sinusoidal drives are kept in the linear regime. The grey data are the references obtained for very weak noise levels.

x n = x 0 n + δx n ,
The verticals are resonance position without Brownian transduction, and lines theoretical calculations (see text).

Transduction mechanism

From the Dung equations, the random motion δx n is transduced into a frequency noise S δω (ω). Since this dependence is quadratic, the frequency noise is neither Gaussian nor centered. Its spectrum depicted in Fig. 1 b) consists in a low frequency part and a high frequency component peaked around 2ω n . The high-frequency uctuations are essentially ltered out by the mode dynamics, as can be seen in a Rotating Wave Approximation (RWA). Thus, driving the mode with a sine wave force F 0 n cos(ωt) weak enough to remain in the linear response limit, the motion x 0 n will adiabatically follow the slow frequency uctuations experiencing both a frequency shift and a spectral broadening [40]. The measurement scheme itself is always slow enough to ensure that all uctuations are spanned while acquiring data.

Note that the Brownian uctuations do not need to be small for the theory to apply.

The phenomenon is non-trivial, and depends strongly on the correlation time of the uctuations τ c = 1/∆ω n . Dening Σ δω = 4β n,n ∆x 2 n a frequency noise ampli- tude parameter (essentially their standard deviation ∝ [ S δω dω] 1/2 ), two regimes should be distinguished depending on the magnitude of the product τ c × Σ δω , see Fig. 1 c). The process can be understood in terms of phase-diusion for the mode studied, the dynamics being averaged over all realizations of the uctuating resonance frequency δω, namely x 0 n (t) ∝ e i t 0 δω(t )dt [40].

The frequency-domain data can thus be described by a convolution of the linear response by a complex valued distribution of frequencies, as seen from the NEMS (bottom of Fig. 1 c):

FT exp(+Γ n t) cosh(a n t)

+ Γn an (1 + 2iα n ) sinh(a n t) (ω) , (3)
FT meaning Fourier Transform, with Γ n = ∆ω n /2 the mode's relaxation rate, a n = Γ n √ 1 + 4iα n and α n = Σ δω 2Γn = τ c × Σ δω the motional narrowing parameter.

By analogy with Nuclear Magnetic Resonance, when τ c × Σ δω 1 the certain component's dynamics is said to be in the motional narrowing limit (MN), while for τ c × Σ δω 1 it lies in the inhomogeneous broadening limit (IB). In the former case, the uctuations are too fast to enable the resolution of the small frequency changes Σ δω [31,32]: the random variable's dynamics looses memory too quickly, and only a fraction of the frequency uctuations impacts the driven motion. This leads to a certain frequency shift which is nothing but the average of the frequency uctuations proportional to ∆x 2 n , together with a (weaker, second order) sym- metric broadening quadratic in ∆x 2 n (bottom-left dis- tribution in Fig. 1 c). In the latter case, the uctuations are slow enough so that the full range of frequency uctuations can be explored by the x 0 n sine-wave response [33,[START_REF] Dykman | [END_REF]: there is a large asymmetric broadening, which reects the actual distribution of frequency uctuations (bottom-right in Fig. 1 c). When mode m = 3 is sine-wave driven and detected while force noise is still applied onto mode n = 1, the treatment is identical with the replacement Σ δω = 2β m,n ∆x 2 n [40]. Besides, an equation similar to Eq. ( 3) holds for the direct calculation of non-linear Brownian spectra [40]. A brief description of the theoretical tools developed in Ref. [40] is given in S.M. [45].

In the next Section, we present the experimental data and the theoretical calculations corresponding to these two situations. The displacement noise spectrum of mode n = 1 is also directly measured. We reach the limit where this spectrum itself is imprinted by the Duing non-linearity [19], and match it to the theory [40].

Since Brownian motions of two m = n distinct modes are not correlated, from these elementary measurements one can then deduce the generic situation where N thermalized modes of the same structure are coupled together.

N o i s e

F r e q . s h i f t ( H z )

∆x 2 1 ( m 2 ) 2 ( 4 β 1 , 1 ∆x 2 1 ) 2 / ∆ω 1 M o d e 1 b )
F r e q . s h i f t ( H z ) F r e q . s h i f t ( H z )

Measured resonance properties

In Fig. 15 a) we present the direct measurement of the Brownian noise spectra S n X (ω) on mode n = 1 for sample 300 µm-n • 2. No sine-wave excitation is applied, neither on n = 1 nor on m = 3 modes. The noise level is quoted in terms of standard deviation ∆x 2 1 instead of T eff (or force noise intensity) since this is the physical parameter of importance. For small Brownian excitations, the peak remains Lorentzian. However, when the amplitude of motion becomes large, the non-linear term β 1,1 starts to impact the lineshape: the peak broadens and becomes asymmetric [19,34]. As expected, the resonance peak globally shifts towards higher frequencies (see Fig. 3 a) for a summary of the spectrum characteristics). The lines are the exact theory from Ref. [40], computed with no free parameters: we call them Duing spectra [45]. Note that no deviations from standard Gaussian statistics are measured in these conditions, as it should be for high-Q devices [44]: spectra on the X quadratures are equivalent to the ones measured on Y, and no cross-correlations are detected [45].

We turn next to the case of the intra-mode coupling.

We still drive mode n = 1 with white noise, but we also measure and drive it with a sine-wave signal. Mode m = 3 is left unexcited. Data and theory from Ref. [40] are compared in Fig. 15 b) with no free parameters. The X lineshapes look like the peaks obtained in the Duing spectrum case, Fig. 15. The eect of the added force noise on the mode is again twofold: rst, the resonance peak slightly shifts towards higher frequencies, and second it broadens (consequently attens) and acquires an asymmetric shape. In Fig. 3 b) we summarize the characteristics of the measured resonance lines on device 300 µm-n • 2 (obtained from the X quadrature).

Measured resonance lines and calculations in the intermode case (sine-wave driving and measuring mode m = 3 while adding force noise on mode n = 1) are shown in Fig. 15 c). They resemble very much the intra-mode results of Fig. 15, even though the quality of the data did not enable to reach as high uctuation levels (see The global agreement between data and theory is remarkable. Essentially, Dung spectra, intra-mode and inter-mode Brownian frequency transduction display the same characteristic features. This highlights that the main ingredient is the dynamics of the noisy mode, not the one of the chosen probe. From Fig. 3, we see that we span the whole range of the phenomenon from motional narrowing to inhomogeneous broadening. In the motional narrowing limit, indeed the rst order eect is a global frequency shift proportional to ∆x 2 n . On the other hand, in the inhomogeneous broadening range the main feature is the asymmetric broadening which is nothing but the image of the frequency distribution (inhomogeneity in time-domain, as opposed to positiondomain for NMR [33]). Further technical discussions of these two limits can be found in S.M. [45].

However, the theory of Ref. [40] applies for sinusoidal excitation strengths lying within the linear response range. When the motion amplitude is increased beyond this limit, new phenomena are expected to take place like e.g. the parametric squeezing of the Brownian motion [43]. One signature obtained experimentally that fails to be reproduced by the theory is shown in Fig. 4: for large sine-wave excitations, the amplitude of the detected mechanical peak lies below the calculation, as if the impact of frequency noise was stronger than expected. In S.M.

[45], we show that the noise spectra measured on mode n are indeed altered by the back-action of the sine-wave response x 0 n ; the X and Y quadratures are not equivalent anymore, and cross-correlations are nonzero at some peculiar frequencies. Further work both theoretical and experimental is required to explore this new dynamical range.

Application to a thermalized family of modes

For a physical thermal bath, the device is always in the motional narrowing limit. In this case, the linear response of mode n to a weak sinusoidal drive remains Lorentzian, with a resonance frequency dressed by the Brownian motion of all modes (global frequency shift proportional to T ). This is essentially analogous to a mechanical Lamb shift [41], in the classical domain. Furthermore the linewidth of the resonance is impacted by a T 2 term, a thermal decoherence eect.

Reproducing results from Ref. [40], these can be written at lowest order in terms of simple expansions, respectively for mode n:

ω n = ω 0 n + 4β n,n ∆x 2 n + m =n 2β n,m ∆x 2 m + m 2 βn,m ∆y 2 m , (4)
∆ω n = ∆ω 0 n + 2 4β n,n ∆x 2 n 2 ∆ω 0 n + m =n 2 2β n,m ∆x 2 m 2 ∆ω 0 m + m 2 2 βn,m ∆y 2 m 2 ∆ω 0 m . (5)
The validity of these expansions has been experimentally veried in the present work for two modes only, Fig. 3.

They can be extended in this simple way to many modes since the Brownian motion between n = m is uncorrelated. For the sake of completeness, we also added the sum over the other family of transverse modes (in y direction), which coecients are designed with a bar, and the index with a prime (the position standard deviation simply writes ∆y 2 m ). Indeed, the nonlinear coupling be- tween exural modes of dierent families has been studied recently [50]. We shall not discuss the coupling to longitudinal and torsional modes, which is outside of the scope of beam mechanics; these depend directly on the Poisson's ratio, and shall be very weak.

Eqs. (4-5) can be easily evaluated for doubly-clamped beams by means of mode parameters calculated using the non-linear extension of Euler-Bernoulli beam theory [38,39,[START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]. With the simple equipartition result

∆x 2 n = k B T /k n , ∆y 2 m = k B T / km we can rewrite these
expressions such that:

ω n -ω 0 n ω 0 n ∝ E beam A 2L 3 (k B T ) (2k 2 n ) , (6)
∆ω n -∆ω 0 When extending these results to the case of a family of modes thermalized to a bath at temperature T , we nd that for typical high-stress top-down structures like the ones used here, the Brownian transduction phenomenon is clearly negligible. However, for much smaller low-stress structures with high aspect-ratio, the eect is foreseen to be limiting in sensing applications [45]. Besides, the certain frequency shift arising from the thermal dressing can also be a source of frequency instability if temperature T is uctuating. Finally when addressing the issue of intrinsic sources of decoherence, it is mandatory to control this phenomenon. Note that other Authors reached the same conclusions with a dierent approach specic to nanotubes [29].

n ∆ω 0 n ∝ E beam A 2L 3 2 (k B T ) 2 (2k 4 n ) Q 2 n , (7)

Supplementary Information for

Non-linear Frequency Transduction of Nano-mechanical Brownian Motion 4 Device, setup and magnetomotive scheme

In Fig. 5 we show two Scanning Electron Micrograph (SEM) pictures of typical 300 µm and 15 µm devices. The gates, used in other works (e.g. [33]) are not connected in the present work. Keeping them oating or grounded is found to be equivalent. Actuation and detection are performed with the magnetomotive scheme [42,46]. The force amplitude acting on mode n in the sine-wave scheme writes F 0 n = ζnLBI0 with I(t) = I0 cos(ωt) the drive current applied, L the length of the beam and B the magnetic eld. ζn is a mode-dependent parameter which can be readily computed from the mode shapes (see Section 5). The detected voltage V (t) is simply the image of the velocity vn(t) of the mode: V (t) = -ζnLB vn(t). We detect it in a standard lock-in scheme, synchronized to the drive frequency ω. We choose for commodity our denitions of the X and Y quadratures of the position to match the phase of the velocity, and thus the one of the detected signal. Note that velocity vn = iω xn and displacement xn require an additional denition: we choose the position of maximum displacement of the beam (in the rst mode, it is simply the motion at the center of the beam). The actual values of the ζn coecients obviously depend on this denition.

In the magnetomotive scheme, the resonance is loaded by the impedance of the external circuit [46]. As a result, the linewidth measured presents a quadratic eld dependence. We show this result for device 300 µm-n • 2 in Fig. 6. As such, the linewidths quoted in the paper are mostly due to the electric circuit: for instance device 300 µm-n • 2 has been measured mainly around 1 Tesla (at 0.84 T, vertical in Fig. 6). But data has also been acquired at 0.6 T and 0.42 T to verify that all features scale properly with the resonance linewidth. From the parabola in Fig. 6 one can deduce the circuit impedance as seen from the NEMS; we obtain about 2 kΩ which corresponds to the bias resistor plus the NEMS metallic layer resistance.

The calibration scheme we use is described in Ref. [42]. The concept on which it is based is very simple: sending a non-resonant current at frequency ω to the device we heat it and track the heating by means of the frequency shift of the resonance (which is T -dependent). Doing this procedure for various ω , and taking the DC heating curve as reference and scaling all the others onto it, we can easily extract the transmission coecient |G(ω )| of the injection line. This technique is very reliable: it is based on a genuinely local property, directly linked to the power injected in the nanomechanical beam.

Since the aluminum on the device is non-superconducting at these temperatures, by measuring its Ohmic contribution we can then infer the transmission of the detection line. The experimental denition of the βi,j coecients needed is shown in Fig. 7, for device 300 µm-n • 1. βn,n is obtained by measuring the position of the resonance maximum as a function of its amplitude squared, for mode n. Driving mode n with a strong sine-wave signal, we measure the quadratic frequency shift induced on mode m detected with a weak sine-wave drive (such that βm,m can be neglected) [38,39,[START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]. The t leads to coecient βm,n. Reversing the procedure, one obtains βn,m.

The parameters kn, mn and βm,n can be compared to the theoretical expectations given in Sec. 8. Typically, kn and mn match within at most ±10 % while the βi,js match within ±20 % (which accounts for a total error on calibrations of order 10 %). In some (rare) cases, discrepancies of about ±50 % for the nonlinear coecients have been seen between experiment and theory. Note also that device 300 µm-n • 1 has a frequency noticeably lower than 300 µm-n • 2: its parameters are consistent with a slightly larger mass, certainly due to an imperfection of the beam (e.g. etching residue of silicon, inhomogeneity in shape). Typical parameters taken for the theoretical evaluation are ρSiN = 3 g/cm 3 , ESiN = 200 GPa, and in-built stress 0.9 GPa (high-stress sample), 120 MPa (low-stress). Note that mathematically, both types of samples are in the so-called high-stress (or string) limit. We take for the metal ρ Al = 2.7 g/cm 3 , E Al = 50 GPa. The stess term includes thermal expansion mismatch. Computed average parameters for the composite beams are obtained from textbook results, e.g. Timoshenko's. We summarize the experimental modal parameters in Tab. A key in our work is the careful calibration of the noise level. We used two voltage AWG (Agilent HP34401A and Tektronix AFG3252), in noise mode, which were both calibrated. The case of the HP34401A is presented in Fig. 8. We rst check that the noise is indeed Gaussian, and centered, with a large bandwidth scope. We nd the conversion factor from applied voltage (in Vpp) to standard deviation σ. This is performed in high-impedance mode, since the setup itself is high impedance (as opposed to 50 Ω). We then verify with a spectrum analyser (and our lock-in technique, see below) that the spectrum is at up to the generator cuto frequency ωc (about 2π× 10 MHz for the Agilent and 150 MHz for the Tektronix). By denition, this cuto is chosen such that the noise spectrum writes SV (ω) = π × σ 2 /ωc in V 2 /(Rad/s) (the π comes from standard Fourier Transform (FT) denitions). However, typically from DC to MegaHertz there is a tiny slope which accounts for about 5 % losses in amplitude. The home made lter is also measured experimentally using a sine-wave signal and both the scope and the lock-in (Fig. 9). It is essentially at from 0.4 MHz to 0.75 MHz, with an insertion loss of -10 % in amplitude (for our 0.6 MHz devices; similar properties apply to the shorter beam resonating around 7 MHz using another lter [10]).

0 . 0 1 0 . 1 1 1 0 1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 1 F i l t e r t

r a n s m i s s i o n

F r e q u e n c y ( M H z ) We can easily convert the noise voltage drive into a noise current δI(t), and then in turn into a noise force δFn(t) = ζnLBδI(t) onto mode n. The force noise correlator is then C n F (τ ) = δFn(t)δFn(t + τ ) and the spectrum

S n F (ω) = FT[C n F ](ω), thus S n F (ω) = (ζnLB) 2 SI (ω)
. By construction, we only need to evaluate SI (ω) within the bandwidth of the lter, where mode n = 1 is excited. Elsewhere, the current is almost zero thus S m F (ω) ≈ 0 for m = n.

In order to measure the Brownian motion, we designed a straightforward technique converting our lock-in into a (phaseresolved) spectrum analyser. The aim of this technique is to be able to keep all our calibrations methods the same from sine-wave measurements to noise measurements, enabling us to be as quantitative as possible. Furthermore, this noise measurement being phase-resolved, it should enable to detect any deviation from perfectly Gaussian statistics. The signal δU (t) that arrives at the level of the lock-in is in fact composed of four terms: (1) the actual δV (t) = -ζnLB δvn(t) to be studied generated by δFn(t), (2)

neglecting the derivatives of slow terms δ İc (t), δ İs (t) with respect to ωδI c (t), ωδI s (t). By denition the susceptibility terms write Re[χn(ω)] = 1 mn (ω 2 n -ω 2 )/ (ω 2 n -ω 2 ) 2 + (∆ωn ω) 2 and Im[χn(ω)] = 1 mn (-∆ωn ω)/ (ω 2 n -ω 2 ) 2 + (∆ωn ω) 2 in the linear response limit, when the Dung term can be neglected. Calling δU c (t) and δU s (t) the two components in parenthesis in Eq. ( 8), it is straightforward to compute the correlators δU c (t)δU c (t + τ ) , δU s (t)δU s (t + τ ) and δU c (t)δU s (t + τ ) and take their Fourier Transform. One readily obtains:

S c,c (ω) = +(ζnLB) 4 ω 2 Re[χn(ω)] 2 + Im[χn(ω)] 2 SI (ω) + (ζnLB) 2 ω -2ReqIm[χn(ω)] + 2 1 Ceqω Re[χn(ω)] SI (ω) + (Req) 2 + 1 Ceqω 2 SI (ω) + S δU noise (ω), (9) S s,s (ω) = S c,c (ω), (10)
S c,s (ω) = 0, (11) having used the properties of the injected noise current, of spectrum SI (ω), white and uncorrelated between the two identical quadratures (the amplier noise is also assumed regular, of spectrum S δU noise [ω]). The sought information is in Eq. ( 9), rst term, with |χn(ω)| 2 = Re[χn(ω)] 2 + Im[χn(ω)] 2 . But the measurement contains also a cross-term between mechanics and background (second term), and a pure background component (last term). Eq. ( 9) can be rewritten:

S c,c (ω) = +(ζnLB) 4 ω 2 |χn(ω)| 2 -2Req(ζnLB) 2 ω Im[χn(ω)] SI (ω) + 2 1 Ceqω (ζnLB) 2 ω Re[χn(ω)] SI (ω) + (Req) 2 + 1 Ceqω 2 SI (ω) + S δU noise (ω). ( 12
)
The rst term is a peak, maximum at resonance, that we call the main term. The second is an imprint of the quadrature signal which is zero on resonance (we call it cross term). Finally the last term is a true background, non-mechanical, independent of eld. Note that because of the loading eect which generates a peak width ∝ B 2 , the two rst components tend towards a constant for large elds, with the main contribution having a stronger B 4 growth.

These components are t on the data in Fig. 10. The t is consistent with a model circuit of a 500 Ω NEMS series resistance in parallel with 500 pF of coaxial line capacitance. These values are within a factor of 2 from actual circuit parameters, as is the overall measured background magnitude. This discrepancy is expected in the present wiring of the setup, because of the direct cross-talk (distributed along the lines) from injection signal to detection signal that aects the background detection level. At large elds, the main contribution is thus essentially ∝ |χn(ω)| 2 , which denes the position spectra S(ω) presented in this work.

In practice, the lock-in measurement is simply equivalent to a phase-resolved spectrum analyser. As such, the bandwidth BW with which the data is acquired is a key parameter. We proceed by digitizing for each given settings (drive, frequency) both X and Y signals demodulated at ω on a Ni-DAQ card which is much faster than the lock-in lter itself (typ. 100 kHz). We take rather long traces (typ. few thousand points N ) from which we can calculate the auto-correlators CX (j) =< Xi Xi+j >, CY (j) =< Yi Yi+j >, and the cross-correlator CXY (j) =< Xi Yi+j > (with < • • • >= i • • • /N using cyclic index notations when i + j > N ). The Stanford SR 844 lter is used with the 24 dB/oct. option, and from the chosen time-constant we compute the eective bandwidth BW of the detector. By denition SX (ω) = π CX (0)/BW , SY (ω) = π CY (0)/BW and SXY (ω) = π CXY (0)/BW . When the statistics is regular, we nd SX (ω) = SY (ω) and SXY (ω) = 0 as we should (see Section 7).

The eect of the bandwidth BW is twofold: on one hand (1) the larger it is, the more noise is collected by the lock-in and thus the larger the detected signal CX (0) is. However, (2) when the bandwidth is large compared to the width of the the background simply scales as the bandwith, as it should for a white noise. On these graphs, we represent the heights and FWHH of the resonance peak normalized to their actual value with respect to the BW of the lock-in detection normalized to the mechanical resonance linewidth ∆ω. We see that the best trade-o is around BW ≈ 0.3 ∆ωn, where signal is rather large without a strong impact on the measured peak.

When the noise is not too large on mode n, the Dung non-linearity is negligible and the mechanical spectrum is a Lorentzian peak. Its area is proportional to Height × FWHH, and is directly related to the standard deviation ∆x 2 n of the mode position uctuations. Furthermore, ∆x 2 n is directly linked to the applied force noise through ∆x 2 n = S n F Qnωn/(2k 2 n ).

This relation is shown in Fig. 13, together with computed eective temperatures.

When the force noise becomes large, then the Dung term is not negligible anymore: the spectrum becomes a Dung spectrum (Fig. 2 of paper and Fig. 14 below). However, for high-Q devices one can show solving the Fokker-Planck equation that the statistics remains Gaussian [44]. Indeed, no anomalies are seen on the measured spectra (we still have SX (ω) = SY (ω) and SXY (ω) = 0, see Section 7). Furthermore, since the measured peak is nothing but the convolution of the mechanical susceptibility with a distribution generated by the uctuations, see Section 6, the area of this peak is preserved: thus the above proportionality relation between measured ∆x 2 n and force spectrum S n F still applies for the distorted spectra. This can be seen experimentally at very large eective temperatures in Fig. 13. Physically, this simply comes from the fact that the induced frequency noise by the non-linear term is purely dispersive, and no extra energy is transferred from the mode to the environment. Note that even though the line distorts, the Height × FWHH remains a rather good estimate of the area of the peak. For the sake of completeness, we should mention that for noise levels above about ∆x 2 1 ≈ 4. 10 -14 m 2 , we do see a Joule heating of the structure due to the large currents injected in the aluminum layer. It is easy to correct for this eect, which produces here only a tiny negative frequency shift which is calibrated in the rst place [42].

Finally, let us comment the range of eective temperatures explored in this experiment. Even though the smallest stochastic drives already correspond to about 100 000 Kelvin, they lie within the motional narrowing limit (see analytical expansions in Section 6, Figs.

3 & 17 and nal discussion on actual thermal bath Section 8). On the other hand, the highest levels of noise used are extremely high and lie in the inhomogeneous broadening range: nonetheless, the basic properties of position uctuations are not altered (Gaussianity, preserved stored energy in the resonance i.e. the expression of the equipartition theorem here). One needs quite peculiar conditions to destroy these properties, and e.g. create squeezed statistics. This is what can be achieved by adding a large sine-wave excitation force on the same mode where large uctuations are present (see Section 7).

Non-linear mode coupling

The basic starting point is the Euler-Bernoulli equation in which the non-linear eect of tensioning has been incorporated [START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF].

It has been discussed extensively in the literature already [38,39], and we just remind the maths here. The equation writes:

ρ beam A ∂ 2 u(z, t) ∂t 2 + η ∂L [u(z, t)] ∂t + E beam Iz ∂ 4 u(z, t) ∂z 4 -T0 + E beam A 2L L 0 ∂u(z, t) ∂z 2 dz ∂ 2 u(z, t) ∂z 2 = ∂F (z, t) ∂z , (13)
with z the axis along the beam, L the length, A = we its cross section (w is the width and e the thickness), Iz its second moment of area. T0 is the stored tension, ρ beam the mass density and E beam the Young's modulus. The parameters can be adapted to the case of a bilayer system (here Al on SiN, see e.g. Timoshenko's textbook).

In Eq. ( 13), ∂F (z, t)/∂z is the overall external force per unit length applied to the beam, and u(z, t) its motion in the x direction. We do not mention any explicit model of mechanical dissipation, and just assume that the friction force is proportional to the time-derivative of some linear functional L of the local displacement u(z, t). We will limit the discussion to two modes, the extension to more being straitghforward, following the works of Refs. [3739,[START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]. Let u(z, t) = xn(t)Ψn(z) + xm(t)Ψm(z) with Ψn(z), Ψm(z) the mode shapes of mode n and m. In this modal decomposition, we take the Ψn, Ψm functions to be normalized at 1 at their maximum. This denes our (time-dependent) amplitude parameters xn(t), xm(t).

Replacing u(z, t) into Eq. ( 13) and projecting on one of the modes (i.e. multiplying by Ψn(z) or Ψm(z) and integrating over the beam), we obtain:

ρ beam ALJn,n ẍn(t) + ρ beam ALJn,m ẍm(t) + ηKn,n ẋn(t) + ηKn,m ẋm(t) + ρ beam AL ω 0 n 2 Jn,n ẍn(t) + ω 0 m 2 Jn,m ẍm(t) + E beam A 2L 3 xn(t) 3 I 2 n,n + xm(t) 3 Im,mIn,m + xn(t) 2 xm(t) (3In,nIn,m) + xm(t) 2 xn(t) In,nIm,m + 2I 2 n,m = F0L ζn cos(ωt), (14)
having assumed ∂F (z, t)/∂z = F0 cos(ωt) independent of z, and used the modal relation 14) applies to the projection on mode n; the projection on m is obtained by inverting n ↔ m. In the ideal case, the modes are orthogonal and Jn,m = Kn,m = In,m = 0 when n = m. In the real case, orthogonality is not perfect; however, writing xi(t) = x 0 i cos(ωt + φ) and using a RWT one immediately sees that many terms in Eq. ( 14) are o-resonant with respect to mode n. It is thus perfectly enough to recast the above result into: mn ẍn(t) + 2λn ẋn(t) + mn ω 0 n 2 xn(t) + kn,n xn(t) 3 + kn,m xn(t)xm(t) 2 = Fn(t), (15) for mode n. An equivalent expression holds for mode m. We have dened the mode parameters mn = ρ beam ALJn,n, λn = ηKn,n/2 together with kn = mn ω 0 n 2 . The projection of the force on the mode is Fn(t) = F 0 n cos(ωt) with F 0 n = F0Lζn. Finally, the non-linear coecients write kn,n

E beam IzΨ i -T0Ψ i = ρ beam A ω 0 i 2 Ψi. ω 0 i is
= E beam A 2L 3 I 2 n,n and kn,m = E beam A 2L 3
In,nIm,m + 2I 2 n,m for n = m. Eq. ( 15) is noting but the driven Dung equation (for mode n), with an extra coupling term to mode m. This is the formalism developed in Refs. [3739,[START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF].

with δωn(t) = βn,m|δ x0 m (t)| 2 = βn,m δXm(t) 2 + δYm(t) 2 , where δXm, δYm are the quadratures of the uctuating mode's motion. Integrating Eq.( 18), we obtain the response of the slow variable in the time domain, which is stochastic:

x0 n (t) = t -∞ F 0 n e iφ 4mnωn exp - ∆ωn 2 + i(ω -ωn) (t -t ) + i t t δω(t )dt dt . (19)
The exponential term is nothing but a stochastic susceptibility in the time domain χ sl (t, t ):

x0 n (t) = t -∞ F 0 n e iφ 4mnωn χ * sl (t, t )dt , (20)
such that the susceptibility in the frequency domain, i.e. the one that is measured, is:

χ(ω) = i 2mnωn +∞ 0 χ sl (t, 0) dt. (21)
The time-domain susceptibility is made stochastic only by the accumulated random phase created by the nonlinearly coupled Brownian motion of the noisy mode. We are then left with averaging over the slow susceptibility:

χ sl (t, 0) = exp - ∆ωn 2 + i(ω -ωn) (t -t ) e i t 0 δω(t )dt . (22)
This averaging procedure is not trivial. The integral is made over a variable which is random at any time between 0 and t, with a nite correlation time, meaning that this frequency noise is highly structured. Thus, to take into account the correlations, one has to average over every path accessible for the accumulated phase in the quadrature space of the uctuating mode between 0 and t. Assuming that the thermal bath is Markovian, the corresponding probability density functional is obtained through Eq. ( 17). The dierent phase paths interfere, with an eciency set by the noisy mode dynamics that appears in the probability density. This so-called path integral approach is detailed in Ref. [40]. Using a discretization procedure for the time interval [0, t], the averaging is then reduced to a cumbersome yet fully analytical calculation of Gaussian integrals, which is performed in Ref. [40]. A similar computation, also detailed in Ref. [40] leads to the averaging of the position spectrum.

In frequency-space these calculations thus lead to convolutions, which are explicitly given from Ref. [40] as:

S n X (ω) = S n Lorentz (ω) * F T    exp(+2Γnt) cosh(ant) + Γn an (1 + 2iαn) sinh(ant) 2    (ω) Noise on n, measure noise on n x0 m = X m Lorentz (ω) * F T exp(+Γmt) cosh(amt) + Γm am (1 + 2iαm) sinh(amt) (ω) Noise on m, measure sine on m x0 m = X m Lorentz (ω) * F T exp(+Γnt) cosh(am,nt) + Γn am,n (1 + 2iαm,n) sinh(am,nt)
(ω) Noise on n, measure sine on m, F T meaning Fourier Transform, with:

Γn = ∆ωn/2 Γm = ∆ωm/2 an = Γn √ 1 + 4iαn am = Γm √ 1 + 4iαm am,n = Γn 1 + 4iαm,n αn = βn,n∆x 2 n Γn αm = 2βm,m∆x 2 m Γm αm,n = βm,n∆x 2 n Γn ,
with S n X the Brownian motion spectrum of mode n, and the measured quadratures X = Im[x 0 m ] and Y = -Re[x 0 m ] of mode m using our denitions (matching the phase of the magnetomotive scheme detection). We write S n Lorentz (ω) and X m Lorentz (ω) the standard Brownian motion spectrum and (complex-valued) response function when no noise is applied. In the RWA (high-Q limit) these are simple Lorentzian peaks. In these formula, the impact of the stochastic motion is given by the motional narrowing parameter α λ (with λ = n, m or n, m). This index describes the dierent situations encountered, respectively: the eect of the Brownian motion on the spectrum itself (Dung spectrum), the coupling of the noise on mode n to the sine-wave response on the same mode n (self-coupling) and the coupling of the noise on one mode n to the sine-wave response of another mode m (mode-coupling). These expressions are used in Fig. 2 of the main article to reproduce the data. We present in Fig. 14 a broader set of measured spectral peaks, and in Fig. 15 a broader set of intra-mode coupling data (putting noise and measuring with a sine wave the same n = 1 exure). Below in Fig. 16 more resonance lines are presented in the case of inter-mode coupling, having the noise on mode n = 1 and measuring with the sine-wave response mode m = 3. To prove the robustness of the eect, the same measurements have been performed on dierent devices: some rather similar (two high-stress 300 µm beams, and a 250 µm) but one quite dierent (a low-stress 15 µm beam). Indeed, the Calculations (full green line) are explained in the text.

parameter α λ depends strongly on the non-linear coecients βn,m (and these depend strongly on the length, see Tab. 1).

The agreement between data and theory in Fig. 17 using these extra devices is eectively as good as in Fig. 3. Note that the range explored in the mode-coupling case is smaller than for self-coupling for pure experimental reasons: the height of the detected peak is smaller, and does not allow to use as high noise intensities.

All data show the same behavior, which is driven simply by the product of the frequency uctuation amplitude (Σ δω )

with the uctuation correlation time (1/∆ωn, the position uctuation's correlation time being 2/∆ωn): this is precisely the parameter α λ called motional narrowing parameter in Ref. [40].

If this parameter is small, we are in the motional narrowing limit: the uctuations are too fast to resolve the full frequency distribution. As such, the function with which we have to convolve the linear response is not strictly speaking a distribution of frequency realization, since it is complex-valued. The result of the phase-diusion process is actually to shift the resonance frequency (by the average of the frequency uctuations, ∝ βm,n∆x 2 n ) at rst order. At second order, it slightly broadens the resonance line without changing its shape: the eect is thus ∝ (βm,n∆x 2 n ) 2 . The rst term is simply the dressing of the mode by the interaction with all the others, while the second one is truly the decoherence eect; the rst one is a certain quantity, while the second one quanties by how much the frequency of the measured mode uctuates.

If the motional narrowing parameter α λ is large, then we are in the inhomogeneous broadening limit. This is the case of Ref. [33] for instance. Note that with an articial telegraph frequency noise generated by a gate electrode, in Ref. [31] Chan and co-workers could switch from one limit to the other not by tuning the amplitude of the noise, but its correlation time instead. In this limit, the full range of uctuations is spanned by the oscillating mode, and one simply measures the averaged motion over all frequency realizations (assuming that the acquisition time is slow enough to capture all uctuations) [START_REF] Dykman | [END_REF].

In this case, the resonance expressions can also be simply described by a convolution with ρ λ (δω), λ = n, m or n, m: F r e q . s h i f t ( H z ) F r e q . s h i f t ( H z ) that home-made lter and setup characteristics were dierent, and had to be calibrated following the same procedure as for 300 µm devices.

write:

ρ(δr) = 1 ∆x 2 n δr exp - 1 2 ∆x 2 n δr 2 Θ [δr] , ρn(δω) = 1 Σ n δω exp - δω Σ n δω , ρm(δω) = 1 Σ m δω exp - δω Σ m δω , ρm,n(δω) = 1 Σ m,n δω exp - δω Σ m,n δω ,
with the standard deviations of the frequency noises dened by:

Σ n δω = 2βn,n∆x 2 n , Σ m δω = 4βm,m∆x 2 n , Σ m,n δω = 2βm,n∆x 2 n .
Above, Θ[x] is the Heaviside step function. The position distribution is a 2D-Gaussian, while the frequency distributions are exponentials; only the couplings are dierent in the above expressions, and in the paper we simply quote Σ δω without index.

The simple averaging procedure works well for large α λ [33]. However, even for moderate motional narrowing parameters it seems to reproduce not too badly the shapes measured (over estimating a bit the broadening); but it fails to capture the certain frequency shift, by construction. For very high Q and small noise levels, the resolution of the measurement can nonetheless be good enough to demonstrate the dierence between the two approaches: this is shown in Fig. 18. Again, we demonstrate very good agreement with the exact theoretical calculation [40]: the motional narrowing eect reduces the impact of the asymmetry of the actual frequency-distribution. F r e q u e n c y ( H Z )

Figure 18: Comparing ts. Data taken on sample 250 µm-n • 1, measuring mode m = 3 while driving uctuations on mode n = 1 (X quadrature on the left, and Y on the right). The applied eld is 1 T for a force F 0 3 = 0.13 pN; the actual mode n = 1 position uctuations is quoted in inset. The full green line is the exact theoretical t, while the dark grey corresponds to the simple average (inhomogeneous broadening theory, here peak shifted by about 15 Hz to match data, see text).

Measured correlators and noise spectra

The theory applied from Ref. [40] is valid for small sine wave excitations: one should remain in the linear response limit.

In this limit, there should be no back-action of the sinusoidal drive onto the statistics of the uctuations, even though their spectra can be altered. So, from the experimental point of view, it is important to check that this limit is satised. We thus rst show in Fig. 19 the raw spectra obtained with no sine wave excitation applied for both small and large noise amplitudes. We see that even if the spectrum is distorted at large Brownian motion levels, we still conrm that X and Y quadratures are equivalent; no cross correlations are detected either. The physical situation is perfectly normal, as it should be for a high Q device [44].

But we expect new phenomena to show up if a strong sine-wave signal is applied, like e.g. noise squeezing [43]. In Fig. 20 we show raw spectra obtained with a rather strong sinusoidal force applied onto mode m = 1 or m = 3, compared to smaller ones (or none). Of course, the data presented in the main part of the paper are obtained with the smallest possible drive levels. What we see is that when mode m = 3 is excited, the noise spectrum measured on n = 1 simply shifts with the amplitude of the sine wave motion. This is nothing but the usual mode-coupling eect [38], but seen on a noise spectrum. From β1,3 = 4. 10 16 Hz/m 2 , we compute a shift of the order of 300 Hz consistent with Fig. 20. No correlations between X and Y are detected, and the two spectra are equivalent.

The situation becomes more interesting when one drives strongly with a sine-wave the same mode where the noise is. We see a peak appearing in the XY correlation, and now spectra measured on X and Y are clearly dierent. We interpret these features as signatures of noise squeezing, as measured in Ref. [43]. We illustrate empirically the eect of large sinusoidal excitation levels on the measured response in Fig. 4 of the paper. The amplitude of the detected mechanical peak lies below the calculation, as if the impact of frequency noise was stronger than expected. Both the usual Dung non-linearity at lowest order with the simple expansions, respectively, reproducing the results of Ref. [40]: The terms in brackets in Eqs. (23)(24) have no dimensions and can be calculated from mode parameters. This is what we describe below. The prefactor gives the strength of the eect from materials properties and geometry. Qn ∝ (E bending + Etension)/E bending . The formulas are summarized in Tab. 2, computed from the mode shapes Ψn(z) in the two limits of interest: low-stress and high-stress (see discussion below). While these expressions of Qn are clearly not enough to t perfectly the experimental ndings (it overestimates the Qs at high n), it is a good starting point since it does not involve any t parameter in the computation of the bracket of Eq. ( 24), and leads to the proper tendency for the sums involved. This will at least produce a very reasonable upper bound for our estimates.

The resonance frequencies verify (ω 0 n ) 2 = kn/mn, and both kn and mn can be readily computed from the mode shapes Ψn(z). The mass writes mn = ρ beam A L 0 Ψn(z) 2 dz, and the spring constant kn = E beam Iz L 0 (∂ 2 Ψn[z]/∂z 2 ) 2 dz + T0 L 0 (∂ 2 Ψn[z]/∂z 2 )Ψn(z)dz (rst term stands for bending energy, second for tensioning; the same integrals are used for the estimate of the Qs). We limit the discussion to perfectly clamped beams. The mode shapes Ψn(z) (and Ψn[z]) are then obtained by solving the (linear) Euler-Bernoulli equation.

Figure 1 :

 1 Figure 1: (Color online). a) The nano-mechanical beams (left) are driven by means of a d.c. magnetic eld B and an a.c. current consisting of the sum of two components (top): a sine-wave which frequency is swept around a chosen mode (n = 1 or n = 3 here) and a Gaussian white noise ltered around mode n = 1. The motion is detected with a lock-in amplier through the induced voltage V , leading to the two quadratures X and Y for each mode n = 1 or n = 3 (right; lines are Lorentz ts and images Ansys R numerical simulations). Data corresponding to beam 300 µm-n • 1 in the linear regime. b) The Gaussian noise force applied onto the mode (here mode 1, center) is equivalent to an eective temperature T eff (right). The motion transduces into a frequency noise (spectrum on the left) because of the Dung non-linearity β1,1 due to tensioning.Only the low frequency part of these uctuations is relevant (in blue, with the d.c. average marked by an arrow), the high frequency term (red) is ltered-out by the dynamics of the mode (adiabatic picture in the rotating frame of the motion). c) Depending on the amplitude of frequency uctuations Σ δω (their standard deviation ∝ [ S δω dω] 1/2 ) with respect to their correlation time τc (here 1/∆ω1, with ∆ω1 the linewidth of the noisy mode) two regimes are distinguished: motional narrowing and inhomogeneous broadening. This is due to the underlying dynamics of phase diusion experienced by the mechanical mode, leading to the averaged frequency distribution depicted below the horizontal arrow (green and violet).

Fig. 1 a

 1 Fig. 1 a) shows a schematic of the setup. For each device and each mode n (or m) studied, we perform a careful calibration based on the technique developed in

with x 0 n

 0 the certain component and δx n the Gaussian and centered random component. In order to introduce the phenomenon, let us rst consider the case depicted in Fig.1 b), where only one mode n is addressed. We apply on n = 1 a Gaussian random force δF n (t) of spectrum S n F (ω) whose strength can be converted into an eective temperature T eff through the Fluctuation-Dissipation Theorem S n F = 2k B T eff m n ∆ω n . ∆ω n is the linewidth of the resonance of mode n (with Q n = ω 0 n /∆ω n the quality factor), and S n F (ω) is white around the mode studied only (and negligible elsewhere). The mechanical mode thus experiences position uctuations (Brownian motion) linked to S n F through the mechanical susceptibil- ity, whose spectrum S n X (ω) is peaked around ω 0 n . Since T eff 4.2 K the experimental temperature, we safely neglect all other sources of uctuations while enabling a thorough tuning of the Brownian motion amplitude of mode n only.

2 F

 2 r e q e n c y ( H z ) D i s p l . ( n m )

Figure 2 :

 2 Figure 2: (Color online). (a) Brownian motion spectra measured on mode n = 1 of sample 300 µm-n • 2 (Dung spectra).

∆ω 1 M o d e 3

 13

Figure 3 :

 3 Figure 3: (Color online). Frequency shift (left) and broadening (from FWHH, right) for (a) the (Dung) spectra measured on mode n = 1 for sample 300 µm-n • 2, (b) the sinewave excitation of mode n = 1, with Brownian motion on the same mode n = 1 for sample 300 µm-n • 2 (intra-mode case), and (c) similar result for sine-wave excitation of mode m = 3, with Brownian motion of mode n = 1 for sample 300 µm-n • 1 (inter-mode). The thin lines are the motional narrowing (MN) analytic expansions, with the dashed verticals corresponding to the cross-over towards inhomogeneous broadening (IB) when τc × Σ δω = 1. The full lines are from the complete theoretical model (see text).

Fig. 3 cFigure 4 :

 34 Fig. 3 c). More data can be found in S.M. [45]. The three basic situations are compared in Fig. 3: we show the characteristics of the measured spectra and resonance lines on 300 µm devices in terms of frequency shift (position of the maximum of the resonance peak) and broadening (measured from the Full Width at Half Height, FWHH). The same characteristics for 250 µm and 15 µm devices are also shown in S.M. [45]: since the non-linear coecients depend strongly on the length L of the structures, this demonstrates the robustness of

  with A = w e the cross-section. In S.M. we summarize the mode parameters obtained in the two extreme limits of Euler-Bernoulli: low-stress (beam) and high-stress (string)[45]. Two key facts have to be highlighted: rst, the prefactor in Eqs.(23)(24) that gives the strength of the eect depends on materials properties and strongly on geometry. Second, increasing the stress in the structure does reduce the sensitivity to Brownian transduction.3 ConclusionBy articially heating a single mode of a NEMS structure, we have demonstrated experimentally the nonlinear frequency transduction of the Brownian motion of this mode onto itself and onto a nearby one. Beyond harmonic mode-coupling [3739], the correlation time τ c of uctuations impacts the dynamics. Two regimes are observed depending on the strength of the stochastic force applied: motional narrowing when the frequency uctuations are small with respect to 1/τ c , and inhomogeneous broadening when they are large. The data are compared to the theory from Ref.[40] that spans the whole range, and we demonstrate excellent agreement without free parameters. To our knowledge, the present work is the rst one presenting a complete experimental analysis of this fundamental (classical) phenomenon, analogous to Nuclear Magnetic Resonance (quantum); eective temperatures up to 10 9 K for the mechanical mode under study have been required to reach the inhomogeneous broadening limit.

Figure 5 :

 5 Figure 5: Devices. SEM pictures of the devices. Left: 300 µm beam (in false colour, the dierent elements of the device are indicated). Right: 15 µm beam. The gates are not used in the present work.

Figure 6 : 2 F 2 F 31 x 2 D

 622312 Figure 6: Loading from electric circuit. Linewidth as a function of magnetic eld for device 300 µm-n • 2. The full line is a quadratic t (see text). The dashed vertical shows were most data have been taken.

Figure 7 :

 7 Figure 7: Nonlinear βi,j coecients. Measurements of the βi,j for device 300 µm-n • 1. The full lines correspond to quadratic ts enabling to extract the numerical values quoted in the paper.

Figure 8 :

 8 Figure 8: Gaussian noise. Characterization of Gaussian voltage noise. Histogram obtained with a large bandwidth scope;with 1 MΩ input. The apparatus was the HP34401A, with 0.5 Vpp/50 Ω settings (delivering thus 1 Vpp on 1 MΩ). The line is a Gaussian t with σ = 0.125 V.

Figure 9 :

 9 Figure 9: Home made lter. Transmission through our home-made lter. The curve is obtained from the ratio of input sine-wave drive amplitude to output sine-wave signal amplitude, in high-impedance mode. The dashed vertical represents the frequency of mode n = 1 of our 300 µm beams. Note that the mode linewidth ∆ω1 is much smaller than the bandwidth of the lter.

Figure 10 :Figure 11 :

 1011 Figure 10: Fit of spectral components vs eld. Data taken on sample 300 µm-n • 2, for a noise drive current of 2.8 10 -19 A 2 /Hz. Fit parameters are discussed below, and can be related to circuit parameters.

Figure 12 :Figure 13 :

 1213 Figure12: Measured peak width versus lock-in bandwidth. FWHH of the motion spectra normalized to its minimum value, with respect to the bandwidth of the lock-in detection normalized to the mechanical linewidth (same conditions as Fig.11). Conventions are the same as in the previous graph.

0000

  the natural resonance frequency of mode i. The parameters introduced above are ζn = L Ψn(z)dz/L, Jn,m = L Ψn(z)Ψm(z)dz/L, Kn,m = L Ψn(z)L [Ψm(z)] dz and In,m = L L Ψn(z) Ψm(z) dz. Eq. (

  q u e n c y ( H z )

Figure 14 :

 14 Figure 14: Brownian spectra in non-linear regime (Color online). Measured displacement spectra S n X (ω) for mode n = (rst exure, no sine-wave drive) on sample 300 µm-n • 2. The standard deviation ∆x 2 1 (i.e Brownian motion level) is increased from top-left to bottom-right, and the lineshape distorts from a pure Lorentzian peak to a Dung spectrum . Lines are theoretical calculations (see text).

Figure 15 :

 15 Figure 15: Intra-coupling ts (Color online). In phase (X) and quadrature (Y) components measured for mode n = 1 (rst exure) on sample 300 µm-n • 2. The sinusoidal force is kept at F1 = 91 fN, while the standard deviation ∆x 2 1 (i.e Brownian motion level) is increased (from top to bottom). Lines are theoretical calculations (see text).

Figure 16 :

 16 Figure 16: Inter-coupling ts. Data taken on sample 300 µm-n • 1, measuring mode m = 3 while driving uctuations on mode n = 1. The applied eld is 0.84 T for a force F 0 3 = 0.6 pN; the actual mode 1 position uctuations are quoted in inset.

  S n X (ω) = S n Lorentz (ω) * ρn(ω) Noise on n, measure noise on n x0 m = X m Lorentz (ω) * ρm(ω) Noise on m, measure sine on m x0 m = X m Lorentz (ω) * ρm,n(ω) Noise on n, measure sine on m, with ρ λ (δω) the frequency distributions directly obtained from ρ(δr) the distribution of position amplitudes. They simply 20

Figure 17 :

 17 Figure 17: Characteristics for 15 µm and 250 µm. Same data as Fig. 3 of the paper, for 2 other devices (see tabular 1). a) self-coupling on 250 µm device. b) self-coupling on 15 µm device. c) mode-coupling on 250 µm device. Note

.

  These are Eqs. (4) and(5) of the main paper, on which the thermal bath discussion is based. The validity of these expansions has been experimentally veried in the present work for two modes only, Figs.3 & 17. They can be extended in this simple way to many modes since the Brownian motion between n = m is uncorrelated. For the sake of completeness, we also added the sum over the other family of transverse modes (in y direction), which coecients are designed with a bar, and the index with a prime (the position standard deviation simply writes ∆y 2 m ). The nonlinear coupling between exural modes of dierent family has been studied recently [D. Cadeddu et al. Nano Letters16, 926 (2016)]. The calculation of Section 5 is easily adapted with u(z, t) = xn(t)Ψn(z) + y m (t) Ψm (z). Note that the mode shapes (and mode numbers) can be a bit dierent along the x, y axes because the beam may not be a perfect square; the second moment of area Iz is a bit dierent for the two families. Hence the bar notation introduced above.Using the integrals In,m over mode shapes introduced in Section 5, and the simple equipartition result ∆x 2 n = kBT /kn, ∆y 2 m = kBT / km we can rewrite these expressions such that:ωn -nIm,m + 2I 2 n,m ) (km/kn) + m (In,n Īm ,m + 2 Ĩ2 n,m ) ( km /kn) , integrals with bar and tilde applying to the other mode family (and cross term): Ψm (z) dz.

Table 1 :

 1 1. sping k n freq. ω n /2π FWHH ∆ω n /2π at 0.84 T non-lin. β n,m /2π n = 1 device 300 µm-n • 1 0.40 N/m 0.59 MHz 100 Hz β 1,1 = 1.35 10 15 Hz/m 2 n = 1 device 300 µm-n • 2 0.45 N/m 0.66 MHz 140 Hz β 1,1 = 1.70 10 15 Hz/m 2 n = 1 device 15 µm-n • 1 1.2 N/m 6.9 MHz 550 Hz β 1,1 = 5.0 10 19 Hz/m 2 n = 3 device 300 µm-n • 1 3.4 N/m 1.74 MHz 40 Hz β 3,1 = 2.75 10 15 Hz/m 2 n = 1 device 250 µm-n • 1 0.45 N/m 0.9 MHz 190 Hz β 1,1 = 8.5 10 15 Hz/m 2 n = 3 device 250 µm-n • 1 4.0 N/m 2.7 MHz 25 Hz β 3,1 = 1.6 10 16 Hz/m 2 Measured

mode parameters relevant to the data presented here. All agree fairly well with analytic computation (see text).

  an Ohmic component ReqδI(t) plus (3) a capacitive component δI(t)/Ceq dt, and nally (4) all other sources of noise δUnoise(t) at the level of the detection which are not correlated to the current δI(t). The total signal δU (t) is thus:

	δU (t) =	-(ζnLB) 2 [+ωIm[χn(ω)]δI c (t) -ωRe[χn(ω)]δI s (t)] + ReqδI c (t) +	δI s (t) Ceqω	+ δU c noise cos(ωt)
	--(ζnLB) 2 [+ωRe[χn(ω)]δI c (t) + ωIm[χn(ω)]δI s (t)] + ReqδI s (t) -	δI c (t) Ceqω	+ δU s noise sin(ωt),

The three last components of the signal are called backgrounds, and arise from the wiring of the experiment (the resistive metallic layer, the transmission line itself from 4.2 K to room temperature and the amplier noise). Req and Ceq are the series transform of the actual setup as seen from the lock-in, and can be estimated from known circuit parameters. Making a RWT on the noise current, we write δI(t) = δI c (t) cos(ωt) -δI s (t) sin(ωt) and similarly for other variables. δvn(t) is readily calculated from the mechanical susceptibility relation δxn(t) = χn(t) * δFn(t).

  Up to date, there is no universal microscopic model explaining nano-mechanical damping, even though many theories exist. They are however phenomenological descriptions of experimental results on SiN beams, e.g. Quirin P. Unterreithmeier et al., PRL 105, 027205 (2010) and A. Suhel et al. APL 100, 173111 (2012). For our devices, we have veried that these approaches apply rather well, see M. Defoort, PhD thesis Université de Grenoble (2014). To make our estimates of thermal decoherence, we thus use the lowest order expression representing a damping proportional to the bending energy

Acknowledgements We thank J. Minet and C. Guttin for help in setting up the experiment, as well as J.-F. Motte, S. Dufresnes and T. Crozes from facility Nanofab for help in the device fabrication. We also warmly thank Mark Dykman for very useful discussions. We acknowledge support from the ANR grant MajoranaPRO No. ANR-13-BS04-0009-01 and the ERC CoG grant ULT-NEMS No. 647917. This work has been performed in the framework of the European Microkelvin Platform (EMP) collaboration.

We turn now to the situation where Fn(t) = F 0 n cos(ωt) + δFn(t) and Fm(t) = F 0 m cos(ωt) + δFm(t), introducing the stochastic variables δFn(t), δFm(t). These are taken as being white, Gaussian, and uncorrelated. The obvious ansatz solving Eq. ( 15) consists in writing xi(t) = x 0 i (t) + δxi(t) for i = n, m, with x 0 i (t) a certain (and eventually sinusoidal) component and δxi(t) a random motion. Replacing in the above equation, we obtain:

with ∆ωn = 2λn/mn the linewidth parameter associated to mode n. We introduced the Dung parameters βn,n = 3 . Note that kn,m = km,n but βn,m = βm,n. The rst line and the third line of Eq. ( 16) correspond to the Dung equation written for the certain component and the stochastic one respectively. The rst bracket is the non-linear coupling onto the certain dynamics: the rst element ∝ x 0 m (t) 2 x 0 n (t) is nothing but the standard mode-coupling term between mode m and n, while the two others correspond to the coupling of the uctuations to the certain position variable. Note the factor of 2 between the coupling to mode m uctuations, and to the same mode n. The second bracket is the equivalent mode-coupling term onto the random motion of mode n: the two rst elements depend on the certain components of m and n, while the last one is due to the uctuations of mode m. Note again the factor of 2 between inter-mode coupling and self-coupling (or intra-mode coupling). In the next section, we explain how to solve Eq. ( 16) for small sine-wave drives and give the explicit theoretical results relevant to the present work, taken from Ref. [40].

6 Adiabatic description of the noise coupling: from motional narrowing to inhomogeneous broadening

The theoretical analysis we use is the one developed in Ref. [40]. Here, we quickly describe the maths leading to the analytic solution of the problem studied. The rst step is to write the dynamics equation Eq. (16) in the frame rotating at the speed of the mode's oscillation. The transformed motion variable x0 n (t) writes x 0 n (t) = (x 0 n (t)Exp[iω 0 n t] + x0 * n (t)Exp[-iω 0 n t])/2. Making a Rotating Wave Approximation (RWA), the component of the frequency uctuations ∝ δxn(t) 2 at 2ωn is removed, keeping only its slow dynamics (and static average). Note that this high-frequency term corresponds somehow to a random parametric pumping for the certain dynamics x 0 n of mode n; but one can show nonetheless that its impact is negligible, and can be safely neglected.

Thus, the dynamics of the certain component adiabatically follows the slow uctuations, which can be interpreted as phase diusion for the resonator. As such, the nal result of the averaging of this phase diusion can be written in time-domain as a multiplication of the mode's decay function by an expression which Fourier Transform we call a frequency distribution (in frequency-domain). This is not strictly speaking a frequency distribution, since in the motional narrowing limit it is complex-valued.

Let us briey remind the reader of the approach followed by Zhang and Dykman in Ref. [40]. To simplify the discussion, we neglect for the time being the intra-mode non-linearity. Nonetheless, the following approach can be extended without loss of generality to the intra-mode case, provided that the driven motion is not too high. Under this assumption, the uctuations of the driven mode and its certain component can be linearly separated. The linearized NEMS' dynamics equations for modes n (driven, i.e. the probe mode) and m (noisy, undriven) in the RWA, are then: S p e c t r a ( V 2 )

F r e q u e n c y ( H z )

F r e q u e n c y ( H z )

Figure 20: Raw spectra, with sine-wave drive. Raw spectra of Brownian motion of mode n = 1 measured with a sinusoidal drive on the third mode m = 3 (left), and on the rst m = 1 (right). Data from sample 250 µm-n • 1 taken at 1 T, with X quadrature in red, and Y in blue. The green inset is the cross-correlation spectrum. The dark color lines are the references with small (or none) sine-wave drive, while the light color corresponds to the strongly driven case. F 0 3 = 0 pN and F 0 3 = 3.1 pN (left), and F 0 1 = 2.7 pN, F 0 1 = 11 pN (right) for small and large settings respectively. For all graphs the motional noise level was ∆x 2 1 = 2.2 10 -17 m 2 . Acquisition bandwidth 78 Hz.

of the mode and the altered statistics are responsible for the apparent saturation of the peak height. A new and dicult theoretical work would be needed to further investigate this very interesting regime.

Extension to thermally induced Brownian noise

When the stochastic driving force is a real thermal bath, the system is always in the motional narrowing limit. In this case, the response is Lorentzian with a dressed frequency, and an additional thermal decoherence. These can be written 23)

for a single family of modes. Here, numerics corresponding to mode n = 1. The overall prefactor is given in inset (same denitions as for Fig. 22). The two limits are depicted (low-stress LS and high-stress HS), with their asymptotic behaviors (dashed lines).

A full analytical solution of the Euler-Bernoulli equation does not exist. However, two limiting cases can be described: low-stress u 1 and high-stress ũ 1 where we dene u = T0L 2 /(E beam Iz) and ũ = E beam Iz/(|T0| L 2 ). Our denitions are u < 0 for tensile, and ũ > 0. In these regimes, analytic shapes can be found as Taylor expansions of u (and √ ũ), and integrated to obtain all the required parameters. We summarize the results in Tabs. 3 and 4 respectively (and Tab. 2 for Q factors). In the low-stress case, we t analytic functions that reproduce rather well the numerical results, and have the proper asymptotic dependence for large n. For high-stress, the expansion is exact at the lowest orders in ũ. This is the limit that applies to all our SiN beams. Note that if all modes would be perfectly orthogonal, the In,m terms would be 0 for n = m (and Īn ,m , Ĩn,m as well). This is not strictly the case, as shown in Tab. 5.

Numerical results for estimates of Q factors of low-stress beams (left,

), assuming losses proportional to bending energy (see text). Tensile means u < 0, and ũ > 0 by construction. Q 0 is a t parameter corresponding to the bending-limited Q factor. Expansions are exact for high-stress, and the t error to numerics is specied for low-stress.

In Figs. 21 and 22 we show the u-dependence of the brackets in Eqs. (23)(24) calculated from the Tabulars for mode n = 1 (rst exure). We compute only a single sum, the one on the same family of modes; for a perfectly squared (monolithic) beam, taking into account the second sum simply amounts to multiply the result by 2. In the more generic case of a bilayer and rectangular beam, the Īn ,m , Ĩn,m should be computed correctly. Thanks to the denominator in the sums, they converge reasonably quickly (much before the high-frequency cut-o that should delimit the validity range of the Euler-Bernoulli equation; namely the atomic size for the wavelength associated to high-frequency modes or the phonon correlation time, whichever comes rst).

Figure 22: Thermal dephasing coecient vs in-built stress. Calculated sum in the bracket of Eq. ( 24) for a single family of modes. Here, numerics corresponding to mode n = 1. The overall prefactor is given in inset (see text). The two limits are depicted (low-stress LS and high-stress HS), with their asymptotic behaviors (dashed lines).

The exact numerical values obtained in Figs. 21 and 22 are of order ≈ 0.1 for low-stress, but fall very quickly with increasing stress. Clearly, only for low-stress devices shall this eect be relevant. The dressing of the resonance frequency is a certain value if T is xed; any temperature instability will translate into a frequency noise which can be calculated from Fig. 21. Furthermore, the strength of thermal decoherence for low-stress devices is essentially given by the prefactor

It goes quadratically with temperature, and since kn ∝ E beam (we 3 /L 3 ) it depends very strongly on the aspect ratio e/L of the structure (see inset of Fig. 22,∝ 8 with V = weL the volume of the beam). It becomes obvious that low-stress bottom-up structures, especially with high-Q, shall be very sensitive to this eect. In practice, it means that an unstressed doubly-clamped nanotube can display a poor spectral Q even if its intrinsic one is high, limiting thus sensing applications. A conclusion also reached by dierent means in Ref. [29].

We remind that our denition of motion amplitude x n is the maximum value along the beam; and tensile means u < 0. We give ts to these coecients as a function of n, with the stated accuracy in the last line.

We remind that our denition of motion amplitude x n is the maximum value along the beam; and ũ > 0 by construction. These expansions are exact.

n, m