Non-linear Frequency Transduction of Nano-mechanical Brownian Motion

TL;DR

This paper investigates how non-linear interactions in nano-mechanical resonators convert Brownian motion into frequency noise, revealing phase diffusion effects and fundamental limits on frequency resolution.

Contribution

It experimentally demonstrates the non-linear transduction of Brownian motion into frequency noise and validates a theoretical model without free parameters.

Findings

Observation of phase diffusion due to non-linear frequency transduction

Identification of two limit behaviors analogous to NMR phenomena

Matching experimental results with theoretical predictions over a wide range

Abstract

We report on experiments addressing the non-linear interaction between a nano-mechanical mode and position fluctuations. The Duffing non-linearity transduces the Brownian motion of the mode, and of other non-linearly coupled ones, into frequency noise. This mechanism, ubiquitous to all weakly-nonlinear resonators thermalized to a bath, results in a phase diffusion process altering the motion: two limit behaviors appear, analogous to motional narrowing and inhomogeneous broadening in NMR. Their crossover is found to depend non-trivially on the ratio of the frequency noise correlation time to its magnitude. Our measurements obtained over an unprecedented range covering the two limits match the theory of Y. Zhang and M. I. Dykman, Phys. Rev. B 92, 165419 (2015), with no free parameters. We finally discuss the fundamental bound on frequency resolution set by this mechanism, which is not marginal for bottom-up nanostructures.

Tables

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

	sping $k_n$	freq. ${\omega }_n/2\pi$	FWHH $\Delta {\omega }_n/2\pi$ at $0.84$T	non-lin. ${\beta }_{n,m}/2\pi$
$n=1$ device 300$\mu$m-n∘1	$0.40$N/m	$0.59$MHz	$100$Hz	${\beta }_{1,1}={\mathrm{1.35\hspace{0.17em}10}}^{15}$Hz/m2
$n=1$ device 300$\mu$m-n∘2	$0.45$N/m	$0.66$MHz	$140$Hz	${\beta }_{1,1}={\mathrm{1.70\hspace{0.17em}10}}^{15}$Hz/m2
$n=1$ device 15$\mu$m-n∘1	$1.2$N/m	$6.9$MHz	$550$Hz	${\beta }_{1,1}={\mathrm{5.0\hspace{0.17em}10}}^{19}$Hz/m2
$n=3$ device 300$\mu$m-n∘1	$3.4$N/m	$1.74$MHz	$40$Hz	${\beta }_{3,1}={\mathrm{2.75\hspace{0.17em}10}}^{15}$Hz/m2
$n=1$ device 250$\mu$m-n∘1	$0.45$N/m	$0.9$MHz	$190$Hz	${\beta }_{1,1}={\mathrm{8.5\hspace{0.17em}10}}^{15}$Hz/m2
$n=3$ device 250$\mu$m-n∘1	$4.0$N/m	$2.7$MHz	$25$Hz	${\beta }_{3,1}={\mathrm{1.6\hspace{0.17em}10}}^{16}$Hz/m2

Table 2. Table 2: Numerical results for estimates of Q 𝑄 Q factors of low-stress beams (left, Q n = Q 0 × [ 𝒬 n 0 + u ​ 𝒬 n 1 ] subscript 𝑄 𝑛 subscript 𝑄 0 delimited-[] superscript subscript 𝒬 𝑛 0 𝑢 superscript subscript 𝒬 𝑛 1 Q_{n}=Q_{0}\times[{\cal Q}_{n}^{0}+u\,{\cal Q}_{n}^{1}] ) and high-stress beams (right, Q n = Q 0 × [ 𝒬 n 0 + u ~ ​ 𝒬 n 1 + u ~ ​ 𝒬 n 2 ] / [ u ~ ​ 𝒬 n 1 + u ~ ​ 𝒬 n 2 ] subscript 𝑄 𝑛 subscript 𝑄 0 delimited-[] superscript subscript 𝒬 𝑛 0 ~ 𝑢 superscript subscript 𝒬 𝑛 1 ~ 𝑢 superscript subscript 𝒬 𝑛 2 delimited-[] ~ 𝑢 superscript subscript 𝒬 𝑛 1 ~ 𝑢 superscript subscript 𝒬 𝑛 2 Q_{n}=Q_{0}\times[{\cal Q}_{n}^{0}+\sqrt{\tilde{u}}\,{\cal Q}_{n}^{1}+\tilde{u}\,{\cal Q}_{n}^{2}]/[\sqrt{\tilde{u}}\,{\cal Q}_{n}^{1}+\tilde{u}\,{\cal Q}_{n}^{2}] ), assuming losses proportional to bending energy (see text). Tensile means u < 0 𝑢 0 u<0 , and u ~ > 0 ~ 𝑢 0 \tilde{u}>0 by construction. Q 0 subscript 𝑄 0 Q_{0} is a fit parameter corresponding to the bending-limited Q 𝑄 Q factor. Expansions are exact for high-stress, and the fit error to numerics is specified for low-stress.

$n$	${𝒬}_n^0$	${𝒬}_n^1$		${𝒬}_n^0$	${𝒬}_n^1$	${𝒬}_n^2$
1	$1$	$-0.03027\cdots$		$1$	$2$	$17.870\cdots$
2	$1$	$-0.015\cdots$		$1$	$2$	$47.480\cdots$
3	$1$	$-0.00760\cdots$		$\cdots$	$\cdots$	$\cdots$
4	$1$	$-0.0045\cdots$		$\cdots$	$\cdots$	$\cdots$
5	$1$	$-0.00300\cdots$		$\cdots$	$\cdots$	$\cdots$
$\cdots$	$\cdots$	$\cdots$		$\cdots$	$\cdots$	$\cdots$
$n$	$1$	$-\frac{0.0161\cdots +0.0322\cdots n+0.101\cdots n^2}{n^3\left(n+2.2\right)}$		$1$	2	$2\left[4+\frac{{\left(n\pi \right)}^2}{2}\right]$
		$\pm$1  % $n\ge 4$

Table 3. Table 3: Numerical results for low-stress beams (see text). The coefficients listed stand for ω n 0 = ( ℒ n 0 + u ​ ℒ n 1 ) 2 ​ E b ​ e ​ a ​ m ​ I z / L 3 / ( ρ b ​ e ​ a ​ m ​ A ​ L ) superscript subscript 𝜔 𝑛 0 superscript superscript subscript ℒ 𝑛 0 𝑢 superscript subscript ℒ 𝑛 1 2 subscript 𝐸 𝑏 𝑒 𝑎 𝑚 subscript 𝐼 𝑧 superscript 𝐿 3 subscript 𝜌 𝑏 𝑒 𝑎 𝑚 𝐴 𝐿 \omega_{n}^{0}=({\cal L}_{n}^{0}+u\,{\cal L}_{n}^{1})^{2}\sqrt{E_{beam}I_{z}/L^{3}/(\rho_{beam}AL)} , k n = ( E b ​ e ​ a ​ m ​ I z / L 3 ) ​ ( 𝒦 n 0 + u ​ 𝒦 n 1 ) subscript 𝑘 𝑛 subscript 𝐸 𝑏 𝑒 𝑎 𝑚 subscript 𝐼 𝑧 superscript 𝐿 3 superscript subscript 𝒦 𝑛 0 𝑢 superscript subscript 𝒦 𝑛 1 k_{n}=(E_{beam}I_{z}/L^{3})({\cal K}_{n}^{0}+u\,{\cal K}_{n}^{1}) , m n = ( ρ b ​ e ​ a ​ m ​ A ​ L ) ​ ( ℳ n 0 + u ​ ℳ n 1 ) subscript 𝑚 𝑛 subscript 𝜌 𝑏 𝑒 𝑎 𝑚 𝐴 𝐿 superscript subscript ℳ 𝑛 0 𝑢 superscript subscript ℳ 𝑛 1 m_{n}=(\rho_{beam}AL)({\cal M}_{n}^{0}+u\,{\cal M}_{n}^{1}) , and I n , n = ℐ n 0 + u ​ ℐ n 1 subscript 𝐼 𝑛 𝑛 superscript subscript ℐ 𝑛 0 𝑢 superscript subscript ℐ 𝑛 1 I_{n,n}={\cal I}_{n}^{0}+u\,{\cal I}_{n}^{1} . We remind that our definition of motion amplitude x n subscript 𝑥 𝑛 x_{n} is the maximum value along the beam; and tensile means u < 0 𝑢 0 u<0 . We give fits to these coefficients as a function of n 𝑛 n , with the stated accuracy in the last line.

$n$	${ℒ}_n^0$	${ℒ}_n^1$	${𝒦}_n^0$	${𝒦}_n^1$	${ℳ}_n^0$	${ℳ}_n^1$	${ℐ}_n^0$	${ℐ}_n^1$
1	$4.73004\cdots$	$-0.02906\cdots$	$198.46\cdots$	$-5.09494\cdots$	$0.396478\cdots$	$-0.0004340\cdots$	$4.8777\cdots$	$0.0006361\cdots$
2	$7.85320\cdots$	$-0.02377\cdots$	$\mathrm{1\hspace{0.17em}669.9}\cdots$	$-20.2173\cdots$	$0.439028\cdots$	$-0.0001160\cdots$	$20.217\cdots$	$0.003449\cdots$
3	$10.9956\cdots$	$-0.01860\cdots$	$\mathrm{7\hspace{0.17em}394.5}\cdots$	$-42.4997\cdots$	$0.505860\cdots$	$0.00051528\cdots$	$50.032\cdots$	$0.05501\cdots$
4	$14.1372\cdots$	$-0.015182\cdots$	$\mathrm{20\hspace{0.17em}082.8}\cdots$	$-77.6713\cdots$	$0.502777\cdots$	$0.00021525\cdots$	$86.270\cdots$	$0.03945\cdots$
5	$17.2788\cdots$	$-0.012794\cdots$	$\mathrm{44\hspace{0.17em}545.4}\cdots$	$-123.614\cdots$	$0.499750\cdots$	$0.00009326\cdots$	$131.93\cdots$	$0.02726\cdots$
6	$20.4203\cdots$	$-0.01104\cdots$	$\mathrm{86\hspace{0.17em}919.8}\cdots$	$-178.034\cdots$	$0.499880\cdots$	$0.00005748\cdots$	$189.12\cdots$	$0.02282\cdots$
7	$23.5619\cdots$	$-0.00971\cdots$	$\mathrm{154\hspace{0.17em}107.5}\cdots$	$-242.210\cdots$	$0.500011\cdots$	$0.00003834\cdots$	$254.02\cdots$	$0.02032\cdots$
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
$n$	$n\pi \left[1+\frac{1}{2n}\right]$	$-\frac{1}{4\pi n}\left[1-\frac{1}{n}\right]$	$\frac{{\left(n\pi \right)}^4}{2}\left[1+\frac{2.2\cdots }{n}\right]$	$-\frac{{\left(n\pi \right)}^2}{2}$	$\frac{1}{2}$	$0$	$\frac{{\left(n+\frac{1}{2}\right)}^2{\pi }^2-2\left(n+\frac{1}{2}\right)\pi }{2}$	$\frac{1}{2}\frac{\left(n+\frac{1}{2}\right)\pi -1}{{\left(n+\frac{1}{2}\right)}^2{\pi }^2}$
	$\pm$0.5  % $n\ge 1$	$\pm$1  % $n\ge 5$	$\pm$1.6  % $n\ge 5$	$\pm$1  % $n\ge 1$	$\pm$0.6  % $n\ge 4$	Negligible	$\pm$1  % $n\ge 3$	$\pm$2  % $n\ge 5$

Table 4. Table 4: Numerical results for high-stress beams (see text). The coefficients listed stand for ω n 0 = ( ℒ n 0 + u ~ ​ ℒ n 1 + u ~ ​ ℒ n 2 ) 2 ​ | T 0 | / L / ( ρ b ​ e ​ a ​ m ​ A ​ L ) superscript subscript 𝜔 𝑛 0 superscript superscript subscript ℒ 𝑛 0 ~ 𝑢 superscript subscript ℒ 𝑛 1 ~ 𝑢 superscript subscript ℒ 𝑛 2 2 subscript 𝑇 0 𝐿 subscript 𝜌 𝑏 𝑒 𝑎 𝑚 𝐴 𝐿 \omega_{n}^{0}=({\cal L}_{n}^{0}+\sqrt{\tilde{u}}\,{\cal L}_{n}^{1}+\tilde{u}\,{\cal L}_{n}^{2})^{2}\sqrt{\left|T_{0}\right|/L/(\rho_{beam}AL)} , k n = ( | T 0 | / L ) ​ ( 𝒦 n 0 + u ~ ​ 𝒦 n 1 + u ~ ​ 𝒦 n 2 ) subscript 𝑘 𝑛 subscript 𝑇 0 𝐿 superscript subscript 𝒦 𝑛 0 ~ 𝑢 superscript subscript 𝒦 𝑛 1 ~ 𝑢 superscript subscript 𝒦 𝑛 2 k_{n}=(\left|T_{0}\right|/L)({\cal K}_{n}^{0}+\sqrt{\tilde{u}}\,{\cal K}_{n}^{1}+\tilde{u}\,{\cal K}_{n}^{2}) , m n = ( ρ b ​ e ​ a ​ m ​ A ​ L ) ​ ( ℳ n 0 + u ~ ​ ℳ n 1 + u ~ ​ ℳ n 2 ) subscript 𝑚 𝑛 subscript 𝜌 𝑏 𝑒 𝑎 𝑚 𝐴 𝐿 superscript subscript ℳ 𝑛 0 ~ 𝑢 superscript subscript ℳ 𝑛 1 ~ 𝑢 superscript subscript ℳ 𝑛 2 m_{n}=(\rho_{beam}AL)({\cal M}_{n}^{0}+\sqrt{\tilde{u}}\,{\cal M}_{n}^{1}+\tilde{u}\,{\cal M}_{n}^{2}) , and I n , n = ℐ n 0 + u ~ ​ ℐ n 1 + u ~ ​ ℐ n 2 subscript 𝐼 𝑛 𝑛 superscript subscript ℐ 𝑛 0 ~ 𝑢 superscript subscript ℐ 𝑛 1 ~ 𝑢 superscript subscript ℐ 𝑛 2 I_{n,n}={\cal I}_{n}^{0}+\sqrt{\tilde{u}}\,{\cal I}_{n}^{1}+\tilde{u}\,{\cal I}_{n}^{2} . We remind that our definition of motion amplitude x n subscript 𝑥 𝑛 x_{n} is the maximum value along the beam; and u ~ > 0 ~ 𝑢 0 \tilde{u}>0 by construction. These expansions are exact.

$n$

${ℒ}_n^0$

${ℒ}_n^1$

${ℒ}_n^2$

${𝒦}_n^0$

${𝒦}_n^1$

${𝒦}_n^2$

${ℳ}_n^0$

${ℳ}_n^1$

${ℳ}_n^2$

${ℐ}_n^0$

${ℐ}_n^1$

${ℐ}_n^2$

1

$1.7724\cdots$

$1.7724\cdots$

$7.0320\cdots$

$4.9348\cdots$

$9.8696\cdots$

$68.4437\cdots$

$0.5$

$-1$

$0$

$4.9348\cdots$

$0$

$-19.739\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$\cdots$

$n$

$\sqrt{n\pi }$

$\sqrt{n\pi }$

$\sqrt{n\pi }\left[\frac{6+{\left(n\pi \right)}^2}{4}\right]$

$\frac{{\left(n\pi \right)}^2}{2}$

${\left(n\pi \right)}^2$

$2{\left(n\pi \right)}^2+\frac{{\left(n\pi \right)}^4}{2}$

$\frac{1}{2}$

$-1$

$0$

$\frac{{\left(n\pi \right)}^2}{2}$

0

$-2{\left(n\pi \right)}^2$

Table 5. Table 5: Numerical results for cross term I n , m subscript 𝐼 𝑛 𝑚 I_{n,m} ( n ≠ m 𝑛 𝑚 n\neq m ) for low-stress beams ( ℐ n , m 0 + u ​ ℐ n , m 1 superscript subscript ℐ 𝑛 𝑚 0 𝑢 superscript subscript ℐ 𝑛 𝑚 1 {\cal I}_{n,m}^{0}+u\,{\cal I}_{n,m}^{1} , left) and high-stress beams ( u ~ ​ ℐ n , m 1 + u ~ ​ ℐ n , m 2 ~ 𝑢 superscript subscript ℐ 𝑛 𝑚 1 ~ 𝑢 superscript subscript ℐ 𝑛 𝑚 2 \sqrt{\tilde{u}}\,{\cal I}_{n,m}^{1}+\tilde{u}\,{\cal I}_{n,m}^{2} , right). Note that I n , m = I m , n subscript 𝐼 𝑛 𝑚 subscript 𝐼 𝑚 𝑛 I_{n,m}=I_{m,n} . Expansions are exact for high-stress, and the fit error to numerics is specified for low-stress. C ​ [ x ] 𝐶 delimited-[] 𝑥 C[x] stands for the Integer part function (Ceil). Tensile means u < 0 𝑢 0 u<0 , and u ~ > 0 ~ 𝑢 0 \tilde{u}>0 by construction.

$n,m$	${ℐ}_{n,m}^0$	${ℐ}_{n,m}^1$		${ℐ}_{n,m}^1$	${ℐ}_{n,m}^2$
1,2	$0$	$0$		$0$	$0$
1,3	$4.35786\cdots$	$0.046108\cdots$		$29.6088\cdots$	$118.435\cdots$
1,4	$0$	$0$		$\cdots$	$\cdots$
1,5	$-3.38985\cdots$	$-0.017825\cdots$		$\cdots$	$\cdots$
1,6	$0$	$0$		$\cdots$	$\cdots$
1,7	$2.71958\cdots$	$0.0079437\cdots$		$\cdots$	$\cdots$
1,8	$0$	$0$		$\cdots$	$\cdots$
1,9	$-2.2558\cdots$	$-0.00374\cdots$		$\cdots$	$\cdots$
$\cdots$	$\cdots$	$\cdots$		$\cdots$	$\cdots$
2,3	$0$	$0$		$0$	$0$
2,4	$8.04755\cdots$	$0.02438\cdots$		$78.957\cdots$	$315.83\cdots$
2,5	$0$	$0$		$\cdots$	$\cdots$
2,6	$-7.11816\cdots$	$-0.005967\cdots$		$\cdots$	$\cdots$
$\cdots$	$\cdots$	$\cdots$		$\cdots$	$\cdots$
3,4	$0$	$0$		$0$	$0$
3,5	$12.243\cdots$	$0.0322\cdots$		$\cdots$	$\cdots$
3,6	$0$	$0$		$\cdots$	$\cdots$
3,7	$-11.5594\cdots$	$-0.01651\cdots$		$\cdots$	$\cdots$
$\cdots$	$\cdots$	$\cdots$		$\cdots$	$\cdots$
4,5	$0$	$0$		$0$	$0$
4,6	$15.6797\cdots$	Negligible		$\cdots$	$\cdots$
4,7	$0$	$0$		$\cdots$	$\cdots$
4,8	$-15.3318\cdots$	Negligible		$\cdots$	$\cdots$
$\cdots$	$\cdots$	$\cdots$		$\cdots$	$\cdots$
$n,m$	$-\frac{8\pi {\left(n+\frac{1}{2}\right)}^2{\left(m+\frac{1}{2}\right)}^2}{\left(n+m+1\right)\left[{\left(n+m+1\right)}^2+{\left(n-m\right)}^2\right]}\times$	$0$		$-nm{\pi }^2\times$	$-4nm{\pi }^2\times$
	$\frac{\left(1+{\left(-1\right)}^{n+m}\right)}{2}{\left(-1\right)}^{C\left[\frac{\left(n-1\right)}{2}\right]}{\left(-1\right)}^{C\left[\frac{\left(m-1\right)}{2}\right]}$			$\frac{\left(1+{\left(-1\right)}^{n+m}\right)}{2}{\left(-1\right)}^{C\left[\frac{\left(n-1\right)}{2}\right]}{\left(-1\right)}^{C\left[\frac{\left(m-1\right)}{2}\right]}$	$\frac{\left(1+{\left(-1\right)}^{n+m}\right)}{2}{\left(-1\right)}^{C\left[\frac{\left(n-1\right)}{2}\right]}{\left(-1\right)}^{C\left[\frac{\left(m-1\right)}{2}\right]}$
	$\pm$11  % $n=1,m\ge 9$	Negligible
	$\pm$6  % $n\ge 2,m>2$

Equations

$${\omega }_n$$

$${\omega }_m$$

$$\mathrm{FT}\left[\frac{\exp (+{\Gamma }_{n}t)}{\cosh (a_{n}t)+\frac{{\Gamma }_n}{a_n}\left(1+2i{\alpha }_n\right)\sinh (a_{n}t)}\right]\left(\omega \right),$$

$${\omega }_n$$

$$\Delta {\omega }_n$$

$$\frac{{\omega }_n-{\omega }_n^0}{{\omega }_n^0}$$

$$\frac{\Delta {\omega }_n-\Delta {\omega }_n^0}{\Delta {\omega }_n^0}$$

$$\delta U(t)$$

$$S^{c,c}(\omega )$$

$$S^{s,s}(\omega )$$

$$S^{c,s}(\omega )$$

$$S^{c,c}(\omega )$$

$${\rho }_{beam}A\frac{{\partial }^{2}u(z,t)}{\partial t^2}+\eta \frac{\partial \mathcal{L}\left[u(z,t)\right]}{\partial t}+E_{beam}{I}_z\frac{{\partial }^{4}u(z,t)}{\partial z^4}-\left[T_0+\frac{E_{beam}A}{2L}{\int }_0^L{\left(\frac{\partial u(z,t)}{\partial z}\right)}^{2}dz\right]\frac{{\partial }^{2}u(z,t)}{\partial z^2}=\frac{\partial F(z,t)}{\partial z},$$

$${\rho }_{beam}AL{J}_{n,n}{\ddot{x}}_n(t)+{\rho }_{beam}AL{J}_{n,m}{\ddot{x}}_m(t)+\eta K_{n,n}{\dot{x}}_n(t)+\eta K_{n,m}{\dot{x}}_m(t)+{\rho }_{beam}AL\left[{\left({\omega }_n^0\right)}^2{J}_{n,n}{\ddot{x}}_n(t)+{\left({\omega }_m^0\right)}^2{J}_{n,m}{\ddot{x}}_m(t)\right]$$

$$m_n{\ddot{x}}_n(t)+2{\lambda }_n{\dot{x}}_n(t)+m_n{\left({\omega }_n^0\right)}^2{x}_n(t)+k_{n,n}{x}_n(t{)}^3+k_{n,m}{x}_n(t)x_m(t{)}^2=F_n(t),$$

$${\ddot{x}}_n^0(t)+\Delta {\omega }_n{\dot{x}}_n^0(t)+{\left({\omega }_n^0\right)}^2{x}_n^0(t)+2{\omega }_n^0\frac{4}{3}{\beta }_{n,n}{x}_n^0(t{)}^3$$

$${\dot{\delta \tilde{x}}}_m(t)+\frac{\Delta {\omega }_m}{2}\delta {\tilde{x}}_m(t)=\frac{\delta F_m(t)e^{-i{\omega }_{m}t}}{2m_m{\omega }_m},$$

$${\dot{\tilde{x}}}_n^0(t)+\left(\frac{\Delta {\omega }_n}{2}+i[\omega -{\omega }_n-\delta \omega (t)]\right){\tilde{x}}_n^0(t)=\frac{F_n^0{e}^{i\varphi }}{4m_n{\omega }_n},$$

$${\tilde{x}}_n^0(t)={\int }_{-\infty }^t\frac{F_n^0{e}^{i\varphi }}{4m_n{\omega }_n}\exp \left(-\left[\frac{\Delta {\omega }_n}{2}+i(\omega -{\omega }_n)\right](t-t^{\prime })+i{\int }_{t^{\prime }}^t\delta \omega (t^{\prime \prime })\mathrm{d}{t}^{\prime \prime }\right)\mathrm{d}{t}^{\prime }.$$

$${\tilde{x}}_n^0(t)={\int }_{-\infty }^t\frac{F_n^0{e}^{i\varphi }}{4m_n{\omega }_n}{\chi }_{sl}^{\ast }(t,t^{\prime })\mathrm{d}{t}^{\prime },$$

$$\chi (\omega )=\frac{i}{2m_n{\omega }_n}{\int }_0^{+\infty }⟨{\chi }_{sl}(t,0)⟩\mathrm{d}t.$$

$$⟨{\chi }_{sl}(t,0)⟩=\exp \left(-\left[\frac{\Delta {\omega }_n}{2}+i(\omega -{\omega }_n)\right](t-t^{\prime })\right)⟨e^{i{\int }_0^t\delta \omega (t^{\prime \prime })\mathrm{d}{t}^{\prime \prime }}⟩.$$

$$S_X^n(\omega )$$

$${\tilde{x}}_m^0$$

$${\Gamma }_n$$

$${\Gamma }_m$$

$$a_n$$

$$a_m$$

$$a_{m,n}$$

$${\alpha }_n$$

$${\alpha }_m$$

$${\alpha }_{m,n}$$

$$S_X^n(\omega )$$

$${\tilde{x}}_m^0$$

$$\rho (\delta r)$$

$${\rho }_n(\delta \omega )$$

$${\rho }_m(\delta \omega )$$

Full text

**Non-linear Frequency Transduction of Nano-mechanical Brownian Motion

**

Olivier Maillet1, Xin Zhou1, Rasul Gazizulin1, Ana Maldonado Cid1, Martial Defoort1,2, Olivier Bourgeois1 & Eddy Collin1,∗

1*Université Grenoble Alpes, Institut Néel - CNRS, F-38042 Grenoble, France

2Now at: CEA, LETI, MINATEC Campus, 17 rue des Martyrs, 38054 Grenoble Cedex 9, France

*: Corresponding Author *

**We report on experiments addressing the non-linear interaction between a nano-mechanical mode and position fluctuations. The Duffing non-linearity transduces the Brownian motion of the mode, and of other non-linearly coupled ones, into frequency noise. This mechanism, ubiquitous to all weakly-nonlinear resonators thermalized to a bath, results in a phase diffusion process altering the motion: two limit behaviors appear, analogous to motional narrowing and inhomogeneous broadening in NMR. Their crossover is found to depend non-trivially on the ratio of the frequency noise correlation time to its magnitude. Our measurements obtained over an unprecedented range covering the two limits match the theory of Y. Zhang and M. I. Dykman, Phys. Rev. B 92, 165419 (2015), with no free parameters. We finally discuss the fundamental bound on frequency resolution set by this mechanism, which is not marginal for bottom-up nanostructures. **

1 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 first place, they serve as ultra-sensitive probes for force sensing [2] with applications e.g. to mass, charge, and even single electronic spin detection [3, 4, 5]. 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 [6], amplification [7] and bit storage [8]. 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 [9, 10, 11].

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 [12, 13]. This leads to a unique platform realizing the “ultimate force detector” foreseen by C. Caves in the 80’s [14]. 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 [15, 16, 17]. 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 [18, 19].

Essentially all applications require in the first 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 fluctuations in nano-mechanical devices becomes an essential technical topic [24, 2, 21, 22, 20, 23]. But in the first place, it is also a fundamental research goal: the measured frequency noise in actual devices is much larger than all expectations [23, 25, 26, 27], demonstrating even non-linear features for carbon-based systems [28, 29]. Thus, attempts have been made to model noise sources [30, 31], or to create model experiments experimentally demonstrating the underlying mechanisms [20, 34, 32, 35].

Clever driving schemes taking advantage of non-linearities have been devised to significantly suppress frequency noise [36, 37]. But what shall be an “ideally frequency-noise minimizing” nano-mechanical system in the first 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 order, these modes are weakly non-linear (so-called Duffing resonators) and are dispersively coupled one to the other [38, 39, 40]. 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 [20, 35, 23, 41, 30], and also on itself. Even in a system realized with ideal materials having no internal sources of noise, this built-in mechanism shall fix an ultimate limit to the mechanical resonance frequencies stability at $T\neq 0$. Only in the limit of $T\rightarrow 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-fluctuations of each of them: a sort of mechanical Lamb shift that dresses all the modes [42].

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 effective temperatures as high as 10^{9}\leavevmode\nobreak\K for this mode only. We extract the effect of this “artificial 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 effect on a nearby mode is measured with the sine-wave excitation scheme. The setup is carefully calibrated [43], 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 different 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” [32]. Based on Ref. [41] and simple expansions of Euler-Bernoulli theory (including non-linear coefficients [39, 40, 49]) 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$ [46]. For bottom-up structures like e.g. carbon nanotubes with high aspect ratio, this limit is not marginal [30].

2 Results

2.1 The nano-electro-mechanical systems

The devices under study are doubly-clamped silicon-nitride nano-beams having width $w=$300  nm and thickness $e_{SiN}=$100  nm. Two high-stress (1  GPa) beams of $L=$300$\leavevmode\nobreak\ \mu$m length have been measured plus a 250$\leavevmode\nobreak\ \mu$m one (samples 300$\leavevmode\nobreak\ \mu$m-n∘1, 300$\leavevmode\nobreak\ \mu$m-n∘2 and 250$\leavevmode\nobreak\ \mu$m-n∘1), together with one low-stress (100  MPa) $L=$15$\leavevmode\nobreak\ \mu$m shorter beam (15$\leavevmode\nobreak\ \mu$m-n∘1). A thin layer of aluminum ($e_{M}$ about 30  nm to 90  nm thick) has been added on top to create electrical contacts. The experiments are performed at 4.2  K in cryogenic vacuum (pressure $<10^{-6}\leavevmode\nobreak\ \mathrm{m}\bbar$).

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 Ref. [43]. 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 coefficients $\beta_{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 fit parameter is indeed an overall correction of the force noise not exceeding 15  % in amplitude (same order as in Ref. [11]). Actuation and detection are performed with the magnetomotive scheme [43, 47]. A drive current (composed of both the Gaussian noise component centered around resonance frequency $\omega_{1}$ and a sine-wave of frequency $\omega$ close to $\omega_{n}$, with $n=1$ or $n=3$) is injected in the NEMS metallic layer through a home-made adder and a $1\leavevmode\nobreak\ \mathrm{k}\mathrm{\SIUnitSymbolOhm}$ bias resistor. In an in-plane d.c. magnetic field orthogonal to the beams, this generates an out-of-plane driving force $F_{n}(t)$ with harmonic component $F_{n}^{0}\cos(\omega 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}(\omega)$ of the Brownian motion of mode $n=1$. Moreover, this enables to measure fluctuations on each of the two quadratures $X$, $Y$ independently (plus their cross-correlations). 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, signatures of squeezed statistics of motion [44] 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. [46].

2.2 Dispersive coupling driven by stochastic motion

Linear motion of thin nano-beams is very well described by the Euler-Bernoulli equation [48]. The basic ingredients involved are the inertia (through the density $\rho_{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$ [40, 39, 49]. This geometric non-linearity results in a tensioning $T_{0}+\delta T$ of the beam with $\delta T\propto x^{2}$, which can be incorporated into the beam equation [49]. 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:

[TABLE]

where we introduced $\omega_{n}^{0}$, $\omega_{m}^{0}$ the linear resonance frequencies, and the Duffing non-linear coefficients $\beta_{i,j}$ [40]. We remind for the interested reader the mathematical derivation of these expressions in S.M. [46].

Eqs. (1-2) can be adapted when one of the motions, say $x_{n}$, is a stochastic variable: $x_{n}=x_{n}^{0}+\delta x_{n}$, with $x_{n}^{0}$ the certain component and $\delta x_{n}$ the Gaussian and centered random component. In order to introduce the phenomenon, let us first consider the case depicted in Fig. 1 b), where only one mode $n$ is addressed. We apply on $n=1$ a Gaussian random force $\delta F_{n}(t)$ of spectrum $S^{n}_{F}(\omega)$ whose strength can be converted into an effective temperature $T_{ef\!f}$ through the Fluctuation-Dissipation Theorem $S^{n}_{F}=2k_{B}T_{ef\!f}m_{n}\Delta\omega_{n}$. $\Delta\omega_{n}$ is the linewidth of the resonance of mode $n$ (with $Q_{n}=\omega_{n}^{0}/\Delta\omega_{n}$ the quality factor), and $S^{n}_{F}(\omega)$ is white around the mode studied only (and negligible elsewhere). The mechanical mode thus experiences position fluctuations (Brownian motion) linked to $S^{n}_{F}$ through the mechanical susceptibility, whose spectrum $S_{X}^{n}(\omega)$ is peaked around $\omega_{n}^{0}$. Since T_{ef\!f}\gg 4.2\leavevmode\nobreak\K the experimental temperature, we safely neglect all other sources of fluctuations while enabling a thorough tuning of the Brownian motion amplitude of mode $n$ only.

2.3 Transduction mechanism

From the Duffing equations, the random motion $\delta x_{n}$ is transduced into a frequency noise $S_{\delta\omega}(\omega)$. 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\omega_{n}$. The high-frequency fluctuations are essentially filtered 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_{n}^{0}\cos(\omega t)$ weak enough to remain in the linear response limit, the motion $x_{n}^{0}$ will adiabatically follow the slow frequency fluctuations experiencing both a frequency shift and a spectral broadening [41]. The measurement scheme itself is always slow enough to ensure that all fluctuations are spanned while acquiring data. Note that the Brownian fluctuations 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 fluctuations $\tau_{c}=1/\Delta\omega_{n}$. Defining $\Sigma_{\delta\omega}=4\beta_{n,n}\Delta x_{n}^{2}$ a frequency noise amplitude parameter (essentially their standard deviation $\propto[\int\!S_{\delta\omega}d\omega]^{1/2}$), two regimes should be distinguished depending on the magnitude of the product $\tau_{c}\times\Sigma_{\delta\omega}$, see Fig. 1 c). The process can be understood in terms of phase-diffusion for the mode studied, the dynamics being averaged over all realizations of the fluctuating resonance frequency $\delta\omega$, namely $x_{n}^{0}(t)\propto\left\langle e^{\displaystyle i\int_{0}^{t}\delta\omega(t^{\prime\prime})\mathrm{d}t^{\prime\prime}}\right\rangle$ [41]. 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):

[TABLE]

FT meaning Fourier Transform, with $\Gamma_{n}=\Delta\omega_{n}/2$ the mode’s relaxation rate, $a_{n}=\Gamma_{n}\sqrt{1+4i\alpha_{n}}$ and $\alpha_{n}=\frac{\Sigma_{\delta\omega}}{2\Gamma_{n}}=\tau_{c}\times\Sigma_{\delta\omega}$ the motional narrowing parameter.

By analogy with Nuclear Magnetic Resonance, when $\tau_{c}\times\Sigma_{\delta\omega}\ll 1$ the certain component’s dynamics is said to be in the “motional narrowing” limit (MN), while for $\tau_{c}\times\Sigma_{\delta\omega}\gg 1$ it lies in the “inhomogeneous broadening” limit (IB). In the former case, the fluctuations are too fast to enable the resolution of the small frequency changes $\Sigma_{\delta\omega}$ [32, 33]: the random variable’s dynamics looses memory too quickly, and only a fraction of the frequency fluctuations impacts the driven motion. This leads to a certain frequency shift which is nothing but the average of the frequency fluctuations proportional to $\Delta x_{n}^{2}$, together with a (weaker, second order) symmetric broadening quadratic in $\Delta x_{n}^{2}$ (bottom-left “distribution” in Fig. 1 c). In the latter case, the fluctuations are slow enough so that the full range of frequency fluctuations can be explored by the $x_{n}^{0}$ sine-wave response [50, 34]: there is a large asymmetric broadening, which reflects the actual distribution of frequency fluctuations (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 $\Sigma_{\delta\omega}=2\beta_{m,n}\Delta x_{n}^{2}$ [41]. Besides, an equation similar to Eq. (3) holds for the direct calculation of non-linear Brownian spectra [41]. A brief description of the theoretical tools developed in Ref. [41] is given in S.M. [46].

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 Duffing non-linearity [20], and match it to the theory [41]. Since Brownian motions of two $m\neq 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.

2.4 Measured resonance properties

In Fig. 15 a) we present the direct measurement of the Brownian noise spectra $S_{X}^{n}(\omega)$ on mode $n=1$ for sample 300$\leavevmode\nobreak\ \mu$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 $\Delta x_{1}^{2}$ instead of $T_{ef\!f}$ (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 $\beta_{1,1}$ starts to impact the lineshape: the peak broadens and becomes asymmetric [20, 35]. 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. [41], computed with no free parameters: we call them “Duffing spectra” [46]. Note that no deviations from standard Gaussian statistics are measured in these conditions, as it should be for high-$Q$ devices [45]: spectra on the X quadratures are equivalent to the ones measured on Y, and no cross-correlations are detected [46].

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. [41] are compared in Fig. 15 b) with no free parameters. The X lineshapes look like the peaks obtained in the “Duffing spectrum” case, Fig. 15. The effect of the added force noise on the mode is again twofold: first, the resonance peak slightly shifts towards higher frequencies, and second it broadens (consequently flattens) and acquires an asymmetric shape. In Fig. 3 b) we summarize the characteristics of the measured resonance lines on device 300$\leavevmode\nobreak\ \mu$m-n∘2 (obtained from the X quadrature).

Measured resonance lines and calculations in the inter-mode 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 fluctuation levels (see Fig. 3 c). More data can be found in S.M. [46].

The three basic situations are compared in Fig. 3: we show the characteristics of the measured spectra and resonance lines on 300$\leavevmode\nobreak\ \mu$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$\leavevmode\nobreak\ \mu$m and 15$\leavevmode\nobreak\ \mu$m devices are also shown in S.M. [46]: since the non-linear coefficients depend strongly on the length $L$ of the structures, this demonstrates the robustness of the effect.

The global agreement between data and theory is remarkable. Essentially, “Duffing 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 first order effect is a global frequency shift proportional to $\Delta x_{n}^{2}$. 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 position-domain for NMR [34]). Further technical discussions of these two limits can be found in S.M. [46].

However, the theory of Ref. [41] 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 [44]. 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. [46], we show that the noise spectra measured on mode $n$ are indeed altered by the back-action of the sine-wave response $x_{n}^{0}$; the $X$ and $Y$ quadratures are not equivalent anymore, and cross-correlations are non-zero at some peculiar frequencies. Further work both theoretical and experimental is required to explore this new dynamical range.

2.5 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 [42], in the classical domain. Furthermore the linewidth of the resonance is impacted by a $T^{2}$ term, a “thermal decoherence” effect.

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

[TABLE]

The validity of these expansions has been experimentally verified 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\neq m$ is uncorrelated. For the sake of completeness, we also added the sum over the other family of transverse modes (in $\vec{y}$ direction), which coefficients are designed with a bar, and the index with a prime (the position standard deviation simply writes $\Delta y_{m^{\prime}}^{2}$). Indeed, the nonlinear coupling between flexural modes of different families has been studied recently [51]. 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 [40, 39, 49]. With the simple equipartition result $\Delta x_{n}^{2}=k_{B}T/k_{n}$, $\Delta y_{m^{\prime}}^{2}=k_{B}T/\mathrm{bar}{k}_{m^{\prime}}$ we can rewrite these expressions such that:

[TABLE]

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) [46]. Two key facts have to be highlighted: first, the prefactor in Eqs. (23-24) that gives the strength of the effect depends on materials properties and strongly on geometry. Second, increasing the stress in the structure does reduce the sensitivity to Brownian transduction.

3 Conclusion

By artificially heating a single mode of a NEMS structure, we have demonstrated experimentally the non-linear frequency transduction of the Brownian motion of this mode onto itself and onto a nearby one. Beyond harmonic mode-coupling [38, 39, 40], the correlation time $\tau_{c}$ of fluctuations impacts the dynamics. Two regimes are observed depending on the strength of the stochastic force applied: “motional narrowing” when the frequency fluctuations are small with respect to $1/\tau_{c}$, and “inhomogeneous broadening” when they are large. The data are compared to the theory from Ref. [41] that spans the whole range, and we demonstrate excellent agreement without free parameters. To our knowledge, the present work is the first one presenting a complete experimental analysis of this fundamental (classical) phenomenon, analogous to Nuclear Magnetic Resonance (quantum); effective temperatures up to 10^{9}\leavevmode\nobreak\K for the mechanical mode under study have been required to reach the “inhomogeneous broadening” limit.

When extending these results to the case of a family of modes thermalized to a bath at temperature $T$, we find 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 effect is foreseen to be limiting in sensing applications [46]. Besides, the certain frequency shift arising from the thermal dressing can also be a source of frequency instability if temperature $T$ is fluctuating. 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 different approach specific to nanotubes [30].

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.

4 Device, setup and magnetomotive scheme

In Fig. 5 we show two Scanning Electron Micrograph (SEM) pictures of typical 300$\leavevmode\nobreak\ \mu$m and 15$\leavevmode\nobreak\ \mu$m devices. The gates, used in other works (e.g. [34]) are not connected in the present work. Keeping them floating or grounded is found to be equivalent.

Actuation and detection are performed with the magnetomotive scheme [43, 47]. The force amplitude acting on mode $n$ in the sine-wave scheme writes $F^{0}_{n}=\zeta_{n}LBI_{0}$ with $I(t)=I_{0}\cos(\omega t)$ the drive current applied, $L$ the length of the beam and $B$ the magnetic field. $\zeta_{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 $v_{n}(t)$ of the mode: $V(t)=-\zeta_{n}LB\,v_{n}(t)$. We detect it in a standard lock-in scheme, synchronized to the drive frequency $\omega$. We choose for commodity our definitions 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 $v_{n}=i\omega\,x_{n}$ and displacement $x_{n}$ require an additional definition: we choose the position of maximum displacement of the beam (in the first mode, it is simply the motion at the center of the beam). The actual values of the $\zeta_{n}$ coefficients obviously depend on this definition.

In the magnetomotive scheme, the resonance is loaded by the impedance of the external circuit [47]. As a result, the linewidth measured presents a quadratic field dependence. We show this result for device 300$\leavevmode\nobreak\ \mu$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$\leavevmode\nobreak\ \mu$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$\Omega$ which corresponds to the bias resistor plus the NEMS metallic layer resistance.

The calibration scheme we use is described in Ref. [43]. The concept on which it is based is very simple: sending a non-resonant current at frequency $\omega^{\prime}$ 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 $\omega^{\prime}$, and taking the DC heating curve as reference and scaling all the others onto it, we can easily extract the transmission coefficient $\left|G(\omega^{\prime})\right|$ 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.

Measuring the height $x_{max}$ of the X peak we can easily obtain the spring constant $k_{n}$ of the mode through $x_{max}=F^{0}_{n}Q_{n}/k_{n}$ with $Q_{n}=\omega^{0}_{n}/\Delta\omega_{n}$ the quality factor of mode $n$. From $(\omega^{0}_{n})^{2}=k_{n}/m_{n}$ we deduce the mass associated to the mode.

The experimental definition of the $\beta_{i,j}$ coefficients needed is shown in Fig. 7, for device 300$\leavevmode\nobreak\ \mu$m-n∘1. $\beta_{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 $\beta_{m,m}$ can be neglected) [40, 39, 49]. The fit leads to coefficient $\beta_{m,n}$. Reversing the procedure, one obtains $\beta_{n,m}$.

The parameters $k_{n}$, $m_{n}$ and $\beta_{m,n}$ can be compared to the theoretical expectations given in Sec. 8. Typically, $k_{n}$ and $m_{n}$ match within at most \pm 10\leavevmode\nobreak\% while the $\beta_{i,j}$s match within \pm 20\leavevmode\nobreak\% (which accounts for a total error on calibrations of order 10  %). In some (rare) cases, discrepancies of about \pm 50\leavevmode\nobreak\% for the nonlinear coefficients have been seen between experiment and theory. Note also that device 300$\leavevmode\nobreak\ \mu$m-n∘1 has a frequency noticeably lower than 300$\leavevmode\nobreak\ \mu$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 \rho_{SiN}=3\leavevmode\nobreak\g/cm3, E_{SiN}=200\leavevmode\nobreak\GPa, and in-built stress 0.9\leavevmode\nobreak\GPa (high-stress sample), 120\leavevmode\nobreak\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 \rho_{Al}=2.7\leavevmode\nobreak\g/cm3, E_{Al}=50\leavevmode\nobreak\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. 1.

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 first check that the noise is indeed Gaussian, and centered, with a large bandwidth scope. We find the conversion factor from applied voltage (in Vpp) to standard deviation $\sigma$. This is performed in high-impedance mode, since the setup itself is high impedance (as opposed to 50  $\Omega$). We then verify with a spectrum analyser (and our lock-in technique, see below) that the spectrum is flat up to the generator cutoff frequency $\omega_{c}$ (about $2\pi\times$ 10  MHz for the Agilent and 150  MHz for the Tektronix). By definition, this cutoff is chosen such that the noise spectrum writes $S_{V}(\omega)=\pi\times\sigma^{2}/\omega_{c}$ in V2/(Rad/s) (the $\pi$ comes from standard Fourier Transform (FT) definitions). However, typically from DC to MegaHertz there is a tiny slope which accounts for about 5  % losses in amplitude. The home made filter is also measured experimentally using a sine-wave signal and both the scope and the lock-in (Fig. 9). It is essentially flat 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 filter [11]).

We can easily convert the noise voltage drive into a noise current $\delta I(t)$, and then in turn into a noise force $\delta F_{n}(t)=\zeta_{n}LB\delta I(t)$ onto mode $n$. The force noise correlator is then $C_{F}^{n}(\tau)=\left\langle\delta F_{n}(t)\delta F_{n}(t+\tau)\right\rangle$ and the spectrum $S_{F}^{n}(\omega)=\mathrm{FT}[C_{F}^{n}](\omega)$, thus $S_{F}^{n}(\omega)=(\zeta_{n}LB)^{2}\,S_{I}(\omega)$. By construction, we only need to evaluate $S_{I}(\omega)$ within the bandwidth of the filter, where mode $n=1$ is excited. Elsewhere, the current is almost zero thus $S_{F}^{m}(\omega)\approx 0$ for $m\neq n$.

In order to measure the Brownian motion, we designed a straightforward technique converting our lock-in into a (phase-resolved) 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 $\delta U(t)$ that arrives at the level of the lock-in is in fact composed of four terms: (1) the actual $\delta V(t)=-\zeta_{n}LB\,\delta v_{n}(t)$ to be studied generated by $\delta F_{n}(t)$, (2) an Ohmic component $R_{eq}\delta I(t)$ plus (3) a capacitive component $\int\delta I(t)/C_{eq}\,dt$, and finally (4) all other sources of noise $\delta U_{noise}(t)$ at the level of the detection which are not correlated to the current $\delta I(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 amplifier noise). $R_{eq}$ and $C_{eq}$ 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 $\delta I(t)=\delta I^{c}(t)\cos(\omega t)-\delta I^{s}(t)\sin(\omega t)$ and similarly for other variables. $\delta v_{n}(t)$ is readily calculated from the mechanical susceptibility relation $\delta x_{n}(t)=\chi_{n}(t)\ast\delta F_{n}(t)$. The total signal $\delta U(t)$ is thus:

[TABLE]

neglecting the derivatives of slow terms $\delta\dot{I}^{c}(t)$, $\delta\dot{I}^{s}(t)$ with respect to $\omega\delta I^{c}(t)$, $\omega\delta I^{s}(t)$. By definition the susceptibility terms write $Re[\chi_{n}(\omega)]=\frac{1}{m_{n}}(\omega_{n}^{2}-\omega^{2})/\left[(\omega_{n}^{2}-\omega^{2})^{2}+(\Delta\omega_{n}\,\omega)^{2}\right]$ and $Im[\chi_{n}(\omega)]=\frac{1}{m_{n}}(-\Delta\omega_{n}\,\omega)/\left[(\omega_{n}^{2}-\omega^{2})^{2}+(\Delta\omega_{n}\,\omega)^{2}\right]$ in the linear response limit, when the Duffing term can be neglected. Calling $\delta U^{c}(t)$ and $\delta U^{s}(t)$ the two components in parenthesis in Eq. (8), it is straightforward to compute the correlators $\left\langle\delta U^{c}(t)\delta U^{c}(t+\tau)\right\rangle$, $\left\langle\delta U^{s}(t)\delta U^{s}(t+\tau)\right\rangle$ and $\left\langle\delta U^{c}(t)\delta U^{s}(t+\tau)\right\rangle$ and take their Fourier Transform. One readily obtains:

[TABLE]

having used the properties of the injected noise current, of spectrum $S_{I}(\omega)$, white and uncorrelated between the two identical quadratures (the amplifier noise is also assumed regular, of spectrum $S_{\delta U_{noise}}[\omega]$). The sought information is in Eq. (9), first term, with $\left|\chi_{n}(\omega)\right|^{2}=Re[\chi_{n}(\omega)]^{2}+Im[\chi_{n}(\omega)]^{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:

[TABLE]

The first 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 field. Note that because of the loading effect which generates a peak width $\propto B^{2}$, the two first components tend towards a constant for large fields, with the main contribution having a stronger $B^{4}$ growth.

These components are fit on the data in Fig. 10. The fit is consistent with a model circuit of a 500$\leavevmode\nobreak\ \Omega$ 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 affects the background detection level. At large fields, the main contribution is thus essentially $\propto\left|\chi_{n}(\omega)\right|^{2}$, which defines the position spectra $S(\omega)$ 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 $\omega$ on a Ni-DAQ card which is much faster than the lock-in filter itself (typ. 100  kHz). We take rather long traces (typ. few thousand points $N$) from which we can calculate the auto-correlators $C_{X}(j)=<X_{i}\,X_{i+j}>$, $C_{Y}(j)=<Y_{i}\,Y_{i+j}>$, and the cross-correlator $C_{XY}(j)=<X_{i}\,Y_{i+j}>$ (with $<\cdots>=\sum_{i}\cdots/N$ using cyclic index notations when $i+j>N$). The Stanford SR 844 filter is used with the 24\leavevmode\nobreak\dB/oct. option, and from the chosen time-constant we compute the effective bandwidth $BW$ of the detector. By definition $S_{X}(\omega)=\pi\,C_{X}(0)/BW$, $S_{Y}(\omega)=\pi\,C_{Y}(0)/BW$ and $S_{XY}(\omega)=\pi\,C_{XY}(0)/BW$. When the statistics is regular, we find $S_{X}(\omega)=S_{Y}(\omega)$ and $S_{XY}(\omega)=0$ as we should (see Section 7).

The effect 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 $C_{X}(0)$ is. However, (2) when the bandwidth is large compared to the width of the mechanical resonance, then the true mechanical signal is “diluted” within the background. As a result, the measured peak height decreases with respect to the background, and broadens; this is shown on both Figs. 11 and 12. At the same time, 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 $\Delta\omega$. We see that the best trade-off is around $BW\approx 0.3\leavevmode\nobreak\ \Delta\omega_{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 Duffing non-linearity is negligible and the mechanical spectrum is a Lorentzian peak. Its area is proportional to Height\leavevmode\nobreak\ \times\leavevmode\nobreak\FWHH, and is directly related to the standard deviation $\Delta x_{n}^{2}$ of the mode position fluctuations. Furthermore, $\Delta x_{n}^{2}$ is directly linked to the applied force noise through $\Delta x_{n}^{2}=S_{F}^{n}\,Q_{n}\omega_{n}/(2k_{n}^{2})$. This relation is shown in Fig. 13, together with computed effective temperatures.

When the force noise becomes large, then the Duffing term is not negligible anymore: the spectrum becomes a “Duffing 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 [45]. Indeed, no anomalies are seen on the measured spectra (we still have $S_{X}(\omega)=S_{Y}(\omega)$ and $S_{XY}(\omega)=0$, see Section 7). Furthermore, since the measured peak is nothing but the convolution of the mechanical susceptibility with a “distribution” generated by the fluctuations, see Section 6, the area of this peak is preserved: thus the above proportionality relation between measured $\Delta x_{n}^{2}$ and force spectrum $S_{F}^{n}$ still applies for the distorted spectra. This can be seen experimentally at very large effective 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\leavevmode\nobreak\ \times\leavevmode\nobreak\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 \Delta x_{1}^{2}\approx 4.\,10^{-14}\leavevmode\nobreak\m2, 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 effect, which produces here only a tiny negative frequency shift which is calibrated in the first place [43].

Finally, let us comment the range of effective temperatures explored in this experiment. Even though the smallest stochastic drives already correspond to about 100\,000\leavevmode\nobreak\Kelvin, they lie within the “motional narrowing” limit (see analytical expansions in Section 6, Figs. 3 & 17 and final 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 fluctuations 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 fluctuations are present (see Section 7).

5 Non-linear mode coupling

The basic starting point is the Euler-Bernoulli equation in which the non-linear effect of tensioning has been incorporated [49]. It has been discussed extensively in the literature already [40, 39], and we just remind the maths here. The equation writes:

[TABLE]

with $\vec{z}$ the axis along the beam, $L$ the length, $A=we$ its cross section ($w$ is the width and $e$ the thickness), $I_{z}$ its second moment of area. $T_{0}$ is the stored tension, $\rho_{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), $\partial F(z,t)/\partial z$ is the overall external force per unit length applied to the beam, and $u(z,t)$ its motion in the $\vec{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 ${\cal 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. [39, 40, 38, 49]. Let $u(z,t)=x_{n}(t)\Psi_{n}(z)+x_{m}(t)\Psi_{m}(z)$ with $\Psi_{n}(z),\Psi_{m}(z)$ the mode shapes of mode $n$ and $m$. In this modal decomposition, we take the $\Psi_{n},\Psi_{m}$ functions to be normalized at 1 at their maximum. This defines our (time-dependent) amplitude parameters $x_{n}(t),x_{m}(t)$.

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

[TABLE]

having assumed $\partial F(z,t)/\partial z=F_{0}\cos(\omega t)$ independent of $z$, and used the modal relation $E_{beam}I_{z}\Psi_{i}^{{}^{\prime\prime\prime\prime}}-T_{0}\Psi_{i}^{{}^{\prime\prime}}=\rho_{beam}A\,\left(\omega_{i}^{0}\right)^{2}\Psi_{i}$. $\omega_{i}^{0}$ is the natural resonance frequency of mode $i$. The parameters introduced above are $\zeta_{n}=\int^{L}_{0}\Psi_{n}(z)dz/L$, $J_{n,m}=\int^{L}_{0}\Psi_{n}(z)\Psi_{m}(z)dz/L$, $K_{n,m}=\int^{L}_{0}\Psi_{n}(z){\cal L}\left[\Psi_{m}(z)\right]dz$ and $I_{n,m}=L\int^{L}_{0}\Psi_{n}(z)^{{}^{\prime}}\Psi_{m}(z)^{{}^{\prime}}dz$. Eq. (14) applies to the projection on mode $n$; the projection on $m$ is obtained by inverting $n\leftrightarrow m$. In the ideal case, the modes are orthogonal and $J_{n,m}=K_{n,m}=I_{n,m}=0$ when $n\neq m$. In the real case, orthogonality is not perfect; however, writing $x_{i}(t)=x_{i}^{0}\cos(\omega t+\phi)$ and using a RWT one immediately sees that many terms in Eq. (14) are off-resonant with respect to mode $n$. It is thus perfectly enough to recast the above result into:

[TABLE]

for mode $n$. An equivalent expression holds for mode $m$. We have defined the mode parameters $m_{n}=\rho_{beam}ALJ_{n,n}$, $\lambda_{n}=\eta K_{n,n}/2$ together with $k_{n}=m_{n}\left(\omega_{n}^{0}\right)^{2}$. The projection of the force on the mode is $F_{n}(t)=F_{n}^{0}\cos(\omega t)$ with $F_{n}^{0}=F_{0}L\zeta_{n}$. Finally, the non-linear coefficients write $k_{n,n}=\frac{E_{beam}\,A}{2L^{3}}I_{n,n}^{2}$ and $k_{n,m}=\frac{E_{beam}\,A}{2L^{3}}\left(I_{n,n}I_{m,m}+2I_{n,m}^{2}\right)$ for $n\neq m$. Eq. (15) is noting but the driven Duffing equation (for mode $n$), with an extra coupling term to mode $m$. This is the formalism developed in Refs. [39, 40, 38, 49].

We turn now to the situation where $F_{n}(t)=F_{n}^{0}\cos(\omega t)+\delta F_{n}(t)$ and $F_{m}(t)=F_{m}^{0}\cos(\omega t)+\delta F_{m}(t)$, introducing the stochastic variables $\delta F_{n}(t),\delta F_{m}(t)$. These are taken as being white, Gaussian, and uncorrelated. The obvious ansatz solving Eq. (15) consists in writing $x_{i}(t)=x_{i}^{0}(t)+\delta x_{i}(t)$ for $i=n,m$, with $x_{i}^{0}(t)$ a certain (and eventually sinusoidal) component and $\delta x_{i}(t)$ a random motion. Replacing in the above equation, we obtain:

[TABLE]

with $\Delta\omega_{n}=2\lambda_{n}/m_{n}$ the linewidth parameter associated to mode $n$. We introduced the Duffing parameters $\beta_{n,n}=\frac{3}{4}\frac{k_{n,n}/m_{n}}{2\omega_{n}^{0}}$ and $\beta_{n,m}=\frac{1}{2}\frac{k_{n,m}/m_{n}}{2\omega_{n}^{0}}$. Note that $k_{n,m}=k_{m,n}$ but $\beta_{n,m}\neq\beta_{m,n}$. The first line and the third line of Eq. (16) correspond to the Duffing equation written for the certain component and the stochastic one respectively. The first bracket is the non-linear coupling onto the certain dynamics: the first element $\propto x_{m}^{0}(t)^{2}x_{n}^{0}(t)$ is nothing but the standard mode-coupling term between mode $m$ and $n$, while the two others correspond to the coupling of the fluctuations to the certain position variable. Note the factor of 2 between the coupling to mode $m$ fluctuations, and to the same mode $n$. The second bracket is the equivalent mode-coupling term onto the random motion of mode $n$: the two first elements depend on the certain components of $m$ and $n$, while the last one is due to the fluctuations 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. [41].

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

The theoretical analysis we use is the one developed in Ref. [41]. Here, we quickly describe the maths leading to the analytic solution of the problem studied. The first 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 $\tilde{x}_{n}^{0}(t)$ writes $x_{n}^{0}(t)=(\tilde{x}_{n}^{0}(t)\mathrm{Exp}[i\omega_{n}^{0}t]+\tilde{x}_{n}^{0\,\ast}(t)\mathrm{Exp}[-i\omega_{n}^{0}t])/2$. Making a Rotating Wave Approximation (RWA), the component of the frequency fluctuations $\propto\delta x_{n}(t)^{2}$ at $2\omega_{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_{n}^{0}$ 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 fluctuations, which can be interpreted as phase diffusion for the resonator. As such, the final result of the averaging of this phase diffusion 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 briefly remind the reader of the approach followed by Zhang and Dykman in Ref. [41]. 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 fluctuations 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:

[TABLE]

[TABLE]

with $\delta\omega_{n}(t)=\beta_{n,m}|\delta\tilde{x}_{m}^{0}(t)|^{2}=\beta_{n,m}\left[\delta X_{m}(t)^{2}+\delta Y_{m}(t)^{2}\right]$, where $\delta X_{m},\delta Y_{m}$ are the quadratures of the fluctuating mode’s motion. Integrating Eq.(18), we obtain the response of the slow variable in the time domain, which is stochastic:

[TABLE]

The exponential term is nothing but a stochastic susceptibility in the time domain $\chi_{sl}(t,t^{\prime})$:

[TABLE]

such that the susceptibility in the frequency domain, i.e. the one that is measured, is:

[TABLE]

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:

[TABLE]

This averaging procedure is not trivial. The integral is made over a variable which is random at any time between [math] and $t$, with a finite 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 fluctuating mode between [math] and $t$. Assuming that the thermal bath is Markovian, the corresponding probability density functional is obtained through Eq. (17). The different phase paths interfere, with an efficiency set by the noisy mode dynamics that appears in the probability density. This so-called path integral approach is detailed in Ref. [41]. 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. [41]. A similar computation, also detailed in Ref. [41] leads to the averaging of the position spectrum.

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

[TABLE]

$FT$ meaning Fourier Transform, with:

[TABLE]

with $S_{X}^{n}$ the Brownian motion spectrum of mode $n$, and the measured quadratures $X=Im[\tilde{x}_{m}^{0}]$ and $Y=-Re[\tilde{x}_{m}^{0}]$ of mode $m$ using our definitions (matching the phase of the magnetomotive scheme detection). We write $S_{Lorentz}^{n}(\omega)$ and $X_{Lorentz}^{m}(\omega)$ 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” $\alpha_{\lambda}$ (with $\lambda=n$, $m$ or $n,m$). This index describes the different situations encountered, respectively: the effect of the Brownian motion on the spectrum itself (“Duffing 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$ flexure). 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 effect, the same measurements have been performed on different devices: some rather similar (two high-stress 300$\leavevmode\nobreak\ \mu$m beams, and a 250$\leavevmode\nobreak\ \mu$m) but one quite different (a low-stress 15$\leavevmode\nobreak\ \mu$m beam). Indeed, the parameter $\alpha_{\lambda}$ depends strongly on the non-linear coefficients $\beta_{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 effectively 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 fluctuation amplitude ($\Sigma_{\delta\omega}$) with the fluctuation correlation time ($1/\Delta\omega_{n}$, the position fluctuation’s correlation time being $2/\Delta\omega_{n}$): this is precisely the parameter $\alpha_{\lambda}$ called “motional narrowing parameter” in Ref. [41].

If this parameter is small, we are in the “motional narrowing” limit: the fluctuations 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-diffusion process is actually to shift the resonance frequency (by the average of the frequency fluctuations, $\propto\beta_{m,n}\Delta x_{n}^{2}$) at first order. At second order, it slightly broadens the resonance line without changing its shape: the effect is thus $\propto(\beta_{m,n}\Delta x_{n}^{2})^{2}$. The first term is simply the “dressing” of the mode by the interaction with all the others, while the second one is truly the “decoherence” effect; the first one is a certain quantity, while the second one quantifies by how much the frequency of the measured mode fluctuates.

If the motional narrowing parameter $\alpha_{\lambda}$ is large, then we are in the “inhomogeneous broadening” limit. This is the case of Ref. [34] for instance. Note that with an artificial telegraph frequency noise generated by a gate electrode, in Ref. [32] 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 fluctuations 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 fluctuations) [50]. In this case, the resonance expressions can also be simply described by a convolution with $\rho_{\lambda}(\delta\omega)$, $\lambda=n,m$ or $n,m$:

[TABLE]

with $\rho_{\lambda}(\delta\omega)$ the frequency distributions directly obtained from $\rho(\delta r)$ the distribution of position amplitudes. They simply write:

[TABLE]

with the standard deviations of the frequency noises defined by:

[TABLE]

Above, $\Theta[x]$ is the Heaviside step function. The position distribution is a 2D-Gaussian, while the frequency distributions are exponentials; only the couplings are different in the above expressions, and in the paper we simply quote $\Sigma_{\delta\omega}$ without index.

The simple averaging procedure works well for large $\alpha_{\lambda}$ [34]. 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 difference between the two approaches: this is shown in Fig. 18. Again, we demonstrate very good agreement with the exact theoretical calculation [41]: the “motional narrowing” effect reduces the impact of the asymmetry of the actual frequency-distribution.

7 Measured correlators and noise spectra

The theory applied from Ref. [41] 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 fluctuations, even though their spectra can be altered.

So, from the experimental point of view, it is important to check that this limit is satisfied. We thus first 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 confirm 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 [45].

But we expect new phenomena to show up if a strong sine-wave signal is applied, like e.g. noise squeezing [44]. 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” effect [39], but seen on a noise spectrum. From \beta_{1,3}=4.\,10^{16}\leavevmode\nobreak\Hz/m2, we compute a shift of the order of 300\leavevmode\nobreak\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 different. We interpret these features as signatures of noise squeezing, as measured in Ref. [44]. We illustrate empirically the effect 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 Duffing non-linearity of the mode and the altered statistics are responsible for the apparent saturation of the peak height. A new and difficult theoretical work would be needed to further investigate this very interesting regime.

8 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 at lowest order with the simple expansions, respectively, reproducing the results of Ref. [41]:

[TABLE]

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 verified 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\neq m$ is uncorrelated. For the sake of completeness, we also added the sum over the other family of transverse modes (in $\vec{y}$ direction), which coefficients are designed with a bar, and the index with a prime (the position standard deviation simply writes $\Delta y_{m^{\prime}}^{2}$). The nonlinear coupling between flexural modes of different family has been studied recently [D. Cadeddu et al. Nano Letters 16, 926 (2016)]. The calculation of Section 5 is easily adapted with $u(z,t)=x_{n}(t)\Psi_{n}(z)+y_{m^{\prime}}(t)\mathrm{bar}{\Psi}_{m^{\prime}}(z)$. Note that the mode shapes (and mode numbers) can be a bit different along the $\vec{x}$, $\vec{y}$ axes because the beam may not be a perfect square; the second moment of area $I_{z}$ is a bit different for the two families. Hence the bar notation introduced above.

Using the integrals $I_{n,m}$ over mode shapes introduced in Section 5, and the simple equipartition result $\Delta x_{n}^{2}=k_{B}T/k_{n}$, $\Delta y_{m^{\prime}}^{2}=k_{B}T/\mathrm{bar}{k}_{m^{\prime}}$ we can rewrite these expressions such that:

[TABLE]

having introduced equivalent integrals with bar and tilde applying to the other mode family (and cross term):

[TABLE]

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 effect from materials properties and geometry.

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 verified 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 $Q_{n}\propto(E_{bending}+E_{tension})/E_{bending}$. The formulas are summarized in Tab. 2, computed from the mode shapes $\Psi_{n}(z)$ in the two limits of interest: low-stress and high-stress (see discussion below). While these expressions of $Q_{n}$ are clearly not enough to fit perfectly the experimental findings (it overestimates the $Q$s at high $n$), it is a good starting point since it does not involve any fit 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 $(\omega_{n}^{0})^{2}=k_{n}/m_{n}$, and both $k_{n}$ and $m_{n}$ can be readily computed from the mode shapes $\Psi_{n}(z)$. The mass writes $m_{n}=\rho_{beam}A\int^{L}_{0}\Psi_{n}(z)^{2}dz$, and the spring constant $k_{n}=E_{beam}I_{z}\int^{L}_{0}(\partial^{2}\Psi_{n}[z]/\partial z^{2})^{2}dz+T_{0}\int^{L}_{0}(\partial^{2}\Psi_{n}[z]/\partial z^{2})\Psi_{n}(z)dz$ (first term stands for bending energy, second for tensioning; the same integrals are used for the estimate of the $Q$s). We limit the discussion to perfectly clamped beams. The mode shapes $\Psi_{n}(z)$ (and $\mathrm{bar}{\Psi}_{n}[z]$) are then obtained by solving the (linear) Euler-Bernoulli equation.

A full analytical solution of the Euler-Bernoulli equation does not exist. However, two limiting cases can be described: low-stress $u\ll 1$ and high-stress $\tilde{u}\ll 1$ where we define $u=T_{0}L^{2}/(E_{beam}I_{z})$ and $\tilde{u}=E_{beam}I_{z}/(\left|T_{0}\right|L^{2})$. Our definitions are $u<0$ for tensile, and $\tilde{u}>0$. In these regimes, analytic shapes can be found as Taylor expansions of $u$ (and $\sqrt{\tilde{u}}$), 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 fit 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 $\tilde{u}$. This is the limit that applies to all our SiN beams. Note that if all modes would be perfectly orthogonal, the $I_{n,m}$ terms would be 0 for $n\neq m$ (and $\mathrm{bar}{I}_{n^{\prime},m^{\prime}}$, $\tilde{I}_{n,m^{\prime}}$ as well). This is not strictly the case, as shown in Tab. 5.

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$ (first flexure). 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 $\mathrm{bar}{I}_{n^{\prime},m^{\prime}}$, $\tilde{I}_{n,m^{\prime}}$ should be computed correctly. Thanks to the denominator in the sums, they converge reasonably quickly (much before the high-frequency cut-off 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 first).

The exact numerical values obtained in Figs. 21 and 22 are of order $\approx 0.1$ for low-stress, but fall very quickly with increasing stress. Clearly, only for low-stress devices shall this effect be relevant. The “dressing” of the resonance frequency is a certain value if $T$ is fixed; 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 $\left(\frac{E_{beam}A}{2L^{3}}\right)^{2}\frac{(k_{B}T)^{2}}{(2k_{n}^{4})}Q_{n}^{2}$. It goes quadratically with temperature, and since $k_{n}\propto E_{beam}(we^{3}/L^{3})$ it depends very strongly on the aspect ratio $e/L$ of the structure (see inset of Fig. 22, $\propto(k_{B}T)^{2}Q_{n}^{2}/(E_{beam}V)^{2}\times(L/e)^{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 effect. 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 different means in Ref. [30].