Effects of noise on the internal resonance of a nonlinear oscillator

We numerically analyze the response to noise of a system formed by two coupled mechanical oscillators, one of them having Duffing and van der Pol nonlinearities, and being excited by a self–sustaining force proportional to its own velocity. This system models the internal resonance of two oscillation modes in a vibrating solid beam clamped at both ends. In applications to nano– and micromechanical devices, clamped–clamped beams are subjected to relatively large thermal and electronic noise, so that characterizing the fluctuations induced by these effects is an issue of both scientific and technological interest. We pay particular attention to the action of stochastic forces on the stability of internal–resonance motion, showing that resonant oscillations become more robust than other forms of periodic motion as the quality factor of the resonant mode increases. The dependence on other model parameters —in particular, on the coupling strength between the two oscillators— is also assessed.


Effects of noise on the internal resonance of a nonlinear oscillator Damián H. Zanette
We numerically analyze the response to noise of a system formed by two coupled mechanical oscillators, one of them having Duffing and van der Pol nonlinearities, and being excited by a selfsustaining force proportional to its own velocity. This system models the internal resonance of two oscillation modes in a vibrating solid beam clamped at both ends. In applications to nano-and micromechanical devices, clamped-clamped beams are subjected to relatively large thermal and electronic noise, so that characterizing the fluctuations induced by these effects is an issue of both scientific and technological interest. We pay particular attention to the action of stochastic forces on the stability of internal-resonance motion, showing that resonant oscillations become more robust than other forms of periodic motion as the quality factor of the resonant mode increases. The dependence on other model parameters -in particular, on the coupling strength between the two oscillators-is also assessed.
Vibrational dynamics of a solid beam is a classical problem in continuum mechanics -already dealt with by Euler and Bernoulli in the eighteenth century 1 -and has important applications in the construction of static and mobile structures, from machines to buildings of all kinds. Although the interest on this problem has hardly faded out and continues to attract the attention of researchers and designers 2 , its study has recently gained a renewed momentum in view of its relevance in the construction of micro-and nanomechanical systems 3 . In particular, time-keeping devices and sensors working at microscopic scales are expected to include oscillating elements consisting of tiny material beams, which are readily built during circuit printing and can be actuated by very small electric fields 4,5 .
A well-known phenomenon in the dynamics of vibrating beams is internal resonance, in which oscillations in a given mode (typically, the main oscillation mode) can excite other modes (typically, higher harmonics) if their frequencies are suitably syntonized with each other 6 . This is an effective mechanism for power transfer toward oscillation modes that might not be easy to excite by an external source, but which are able to efficiently couple with intermediary modes 7,8 . Recently, internal resonance has been proposed as a mechanism for frequency stabilization in nonlinear micromechanical oscillators formed by silica beams clamped at their two ends (clamped-clamped, or c-c, beams 9 ). Vibrations of a c-c beam are well described by the Duffing equation, i.e. as a nonlinear oscillator with a cubic term in its restoring force 10,11 . Nonlinearity implies that the oscillation frequency and amplitude depend on each other (a-f effect 12 ). This may hinder the application to time-keeping devices, as uncontrolled variations of amplitude -for instance, due to thermal or electronic noise-would bring about an undesired change in the frequency. However, it has been discovered that if the system operates within internal resonance, a regime develops in which the frequency becomes almost independent of the amplitude, and the a-f effect is thus suppressed 9 .
Micro-and nanomechanical devices usually work in highly fluctuating environments, where thermal and electronic noise is an unavoidable and significant dynamical ingredient. With respect to other forces involved in these systems, in fact, the relative effect of noise grows as the relevant length scales become smaller. Although the role of fluctuations in the operation of small-scale devices is widely acknowledged [12][13][14] , it does not seem that their action on the dynamics of micro-and nanomechanical oscillators has been systematically studied. In this paper, we numerically analyze the response of a vibrating c-c beam to noise, with emphasis in oscillations within the internal-resonance regime, although noise-induced transitions to other states are also considered. The effects of stochastic forces on a Duffing-like vibrating (macroscopic) beam has recently been characterized, both numerically and experimentally 15 , but the response of internal resonance to noise was not addressed in those previous studies. Here, we consider a c-c beam in a closed-loop configuration, whose oscillations are maintained by a self-sustaining force. This is the standard operational configuration in any time-keeping device 16 . The selfsustaining force is generated by conditioning a signal electronically read from the oscillation itself, and then reinjected as an "external" excitation. In our case, this force is proportional to the oscillation velocity delayed by a prescribed time. Such kind of feedback is able to generate self-sustained stationary oscillations if a van der Pol-like nonlinearity 17 is also present in the damping force 18,19 . The reinjection of the self-sustaining force allows the mechanical oscillator to attain stable periodic motion, characteristic of a wide class of dynamical systems spanning from plasma physics 20 to neurosciences 21 . The interaction between the resonant oscillation modes is described using a well-tested model 6,8,9 in which modes are represented as mutually coupled one-dimensional oscillators.
The main goal of our work is to identify the most important factors that control the response of the oscillator to stochastic forces, determining the dominant dependences of the fluctuations induced by noise on the several parameters that define the system. Many of these parameters, in fact, are not necessarily amenable to experimental control, and numerical analysis is therefore necessary to assess their role in the present problem. In the next section, we introduce the two-oscillator model and discuss the main aspects of its dynamics in the absence of noise. Then, we describe our numerical methodology and present results under the action of stochastic forces. The last section is devoted to a qualitative discussion of the results and to a few concluding remarks.

Two-oscillator model and noiseless dynamics
As a model for the interaction between two oscillation modes of a clamped-clamped (c-c) beam, we consider two linearly coupled one-dimensional oscillators, described by coordinates x 1 (t) and x 2 (t) 6 . The first oscillator is subjected to a self-sustaining force proportional to its own velocity at a delayed time, , and is affected by a cubic (Duffing 11 ) nonlinearity in the restoring force and a quadratic (van der Pol 17 ) nonlinearity in the damping force 18,19 . This oscillator represent the main oscillation mode of the c-c beam. The second oscillator, which represents a higher-harmonic mode, is linear. Although moderate Duffing-like nonlinear response in higher-harmonic c-c beam vibration has been reported in experiments 9,22 , we show below that significant response of the higher-harmonic mode under conditions of internal resonance is limited to a narrow frequency interval, where higher-harmonic amplitude-frequency interdependence has no room to manifest itself. A linear higher-harmonic mode is therefore expected to satisfactorily describe the phenomenology focused on in the present contribution. Recently, the effects of higher-harmonic nonlinearity on the noiseless dynamics of vibrating c-c beams have been analyzed theoretically, showing that internal resonance preserves its main qualitative features over wide parameter ranges 23 .
The equations of motion derived for the model are with γ 1,2 the damping coefficients per unit mass, ω 1,2 2 the natural frequencies, β the Duffing coefficient, and g the gain of the self-sustaining force. The coordinates x 1,2 have been rescaled in such a way that the coefficient of x 1 2 in the damping force of oscillator 1 has a prescribed value, and that the coupling strength j is the same for both oscillators. From now on, moreover, we fix ω 1 ≡ 1 by a suitable choice of time units. The stochastic forces ξ 1,2 are specified in the next section, at the stage of their numerical implementation.
In the noiseless case, ξ 1,2 ≡ 0, approximate stationary solutions to the equations of motion can be found by the usual technique of proposing harmonic oscillations, = Ω x t A t ( ) cos 1 1 and , and neglecting higher-harmonic contributions coming from nonlinear terms. For equations (1) and (2), this procedure yields with β β =  3 /4. These are equations for the amplitudes A 1,2 , the stationary frequency Ω, and the phase ψ. The main panel of Fig. 1A shows the resulting interdependence of A 1 vs Ω for γ 1 = 0.001, γ 2 = 0.003, ω 2 = 1.5, β = 0.04, g = 0.1, and j 2 = 0.001. The curve is parametrized by the time delay τ. We recognize the overall leaning shape of the Duffing resonance peak 11 . However, in contrast with the standard harmonically-forced Duffing oscillator, the present system does not have a low-amplitude solution for high frequencies 18 . This is a consequence of the linear dependence of the self-sustaining force on the velocity. Our main interest, nevertheless, lies in the zone where the frequency of stationary oscillations and the natural frequency of oscillator 2 are close to each other, Ω ≈ ω 2 . The lower-right inset in Fig. 1A is a close-up of that zone, where the resonance peak exhibits a gap in the frequency Ω 9,24 . Within the gap, no solution exists. The solutions at each side of the gap correspond to maximal oscillation amplitude for oscillator 2, as shown in the upper-left inset. These solutions represent the internal resonance of the two oscillation modes, for which the energy transfer from oscillator 1 to oscillator 2 is most efficient. For sufficiently large β, g, and small γ 1 -i.e. when the nonlinearity of oscillator 1 is well developed, with a markedly leaning, long, narrow Duffing peak-the width of the gap at Ω ≈ ω 2 is mainly controlled by the damping of oscillator 2 and the strength of coupling between the oscillators, through the combination γ Γ = − j 2 2 . For large values of this parameter, damping is too large and coupling is too weak to allow for resonance, and the gap is absent. As Γ becomes smaller, the gap is opened and initially widens but, when Γ decreases further, the width begins to decrease as well. For Γ small enough, the gap width is proportional to a power of this parameter. Taking the remaining parameters as in Fig. 1, the gap appears for γ ≈ . 0 0077 2 (Γ ≈ 7.7). However, for γ 2 = 0.003 -as in the figure-its width is already in the zone where it varies proportionally to a power of Γ. Figure 1B shows the solutions to equations (3) to (6) for the stationary frequency Ω as a function of the time delay τ, while Fig. 1C shows the corresponding solutions for A 1 vs τ. Full and dashed curves indicate stable and unstable solutions, respectively. We stress that negative values of τ -which, naturally, cannot be realized in an experiment-are to be interpreted, for stationary harmonic motion, as positive time delays just below the oscillation period, i.e. generating phase shifts  τ π Ω 2 . The almost horizontal branches within the frequency gap at Ω ≈ ω 2 , in Fig. 1B, demonstrate the phenomenon of frequency stabilization associated with the internal resonance, discussed in the Introduction. We begin next section studying the transitions induced by noise between stable states inside and outside the gap, for fixed delay, illustrated by the vertical double arrows at τ = 0.5 both in Fig. 1B and C.
The stability of stationary solutions has been determined using a variation of the multiple-scale method 6 , which assumes that the time scale associated with the relaxation of oscillation amplitudes is much longer than the oscillation period. This is in fact the case, in particular, for experiments with micromechanical oscillators, whose weak damping (high quality factor, ∼ − Q 10 10 4 5 ) insures a net separation between the two time scales. Around a stationary oscillation of frequency Ω, the amplitudes A 1,2 (t) and phases φ 1,2 (t) of the coordinates ( 1 s in ) c os( ) , As can be verified by direct inspection, the stationary solutions to these equations coincide with the solutions to equations (3) to (6) with ψ = φ 2 − φ 1 . Within the multiple-scale approximation, linearization of equations (7) to (10) makes it possible to decide on the stability of the stationary oscillations. Also, this analysis provides the relaxation rates toward the stationary state which, as shown in the next section, play a role in the response of the system to the effects of noise.

Response to stochastic forces
Inside the internal-resonance gap, the close proximity of a stable and an unstable solution for fixed τ, as illustrated by Fig. 1B (for instance, for τ = 0.5), naturally rises the question on how robust against the effects of noise the system is, if prepared to oscillate in the stable solution. Is the system able to resist stochastic fluctuations and remain inside the gap, or will it immediately migrate to the stable oscillation of higher frequency which, for the same value of τ, lies outside the gap? This question becomes particularly relevant if, as advanced in the Introduction, internal resonance is sought to be exploited as a way of stabilizing the frequency of the nonlinear oscillator against fluctuations of the oscillation amplitude.
To address this problem, we have solved equations (1) and (2) numerically, using a standard explicit secondorder Runge-Kutta scheme adapted to stochastic differential equations 25 . The stochastic forces ξ 1,2 (t) have been implemented by interpreting the equations for the velocities in the Ito form, dv i = F i dt + σ i dtW i , where F i are the deterministic forces and dW i are differential Wiener processes with 〈dW i 〉 = 0 and 〈 〉 = dW dt i 2 . The coefficients σ i weight the strength of noise. The stochastic differential equation is thus discretized (to the first order) as 26 . In our calculations, successive values of the stochastic variable ζ i have been drawn from a Gaussian distribution with zero mean and unitary variance, using the pseudorandom number generator RANDGEN 27 .
For the numerical runs, we have chosen initial conditions that led, in the absence of noise, to a selected stable solution -inside or outside the gap. We have first run the system without noise until it closely approached the selected solution and, then, the stochastic forces were turned on. Along each realization, we have recorded the successive local maxima of the coordinates x 1,2 (t), as a measure of the amplitudes A 1,2 affected by noise. In the presentation of our results, however, we focus on the behavior of A 1 , since we expect the amplitude of the main oscillation mode of the c-c beam -represented in our model by oscillator 1-to be most readily accessible to experimental measurements [7][8][9] . Figure 2 shows a few typical realizations for the amplitude of oscillator 1, with two values of σ 1 and σ 2 = 0. The other parameters are γ 1 = 0.001, γ 2 = 0.004, ω 2 = 1.5, β = 0.04, g = 0.1, and j 2 = 0.001, with time delay τ = 0.5. In Fig. 2A and B, σ 1 = 0.2. Figure 2A depicts two time series for A 1 , starting from initial conditions in the vicinity of either stable solution, inside (dark gray, lower series) and outside (light gray, upper series) the gap. Green lines show time averages over 300 successive values of A 1 . In Fig. 2B, normalized histograms of the values of A 1 in each series are shown. We see that for, this noise level, both solutions fluctuate around welldefined values, with a larger dispersion for the oscillation inside the gap. In Fig. 2C and D, the stochastic force on oscillator 1 has been increased to σ 1 = 0.3. Here, from any initial condition in the vicinity of the stable solutions, the amplitude A 1 exhibits occasional transitions between the two states. For these parameters, however, the signal is most of the time closer to the upper state. This is also evident from the corresponding histogram. A quantitative characterization of the fluctuations in the oscillation amplitude is given by the mean value and standard deviation of A 1 along the realizations described in the preceding paragraph. The leftmost panels in Fig. 3 show these averages for the same parameters of Fig. 2, as functions of the noise intensity σ 1 . Full and empty dots correspond to initial conditions that lead to the stable states inside and outside the internal resonance gap, respectively. For weak noise,  σ . 0 25, the corresponding mean values of A 1 are well separated from each other, while the standard deviations grow linearly with σ 1 , as demonstrated by the straight line between the two datasets. On the other hand, for sufficiently strong noise, σ > . ∼ 0 5, the two datasets collapse for both the mean amplitude and the standard deviation. Under the effect of such large fluctuations, the amplitude randomly visits a large interval around A 1 ≈ 7.7, comprising both stable states, with a standard deviation that still grows proportional to σ 1 . For intermediate noises, in the zone highlighted in Fig. 3 by the vertical stripe, the amplitude shows transitions between the two states, as illustrated in Fig. 2C, with a prevalence of the higher-amplitude solution. In the leftmost half of the stripe, it is still possible to measure clearly defined values for the mean amplitude and the standard deviation for each state, as the system spends well-separated periods in their respective vicinities (as illustrated in Fig. 2C). As σ 1 grows, however, the separation becomes less defined, and the two quantities merge into single values. Within this zone, the standard deviation clearly abandons its linear dependence on σ 1 .
The panels to the right of Fig. 3 show the same results as to the left but for a lower value of the damping coefficient of oscillator 2, γ 2 = 0.003. While the overall behavior is the same as for γ 2 = 0.004, now it is the state of lower amplitude which prevails over the other one during the transitions within the vertical stripe. Thus, the system spends more time inside the internal resonance gap. Note also that the vertical stripe is narrower than for larger γ 2 , and that the value of σ 1 at which the two datasets merge is slightly shifted to the right, from σ 1 ≈ 0.35 to 0.4. As γ 2 decreases further, the prevalence of the stationary state inside the gap grows. For γ 2 = 0.002 it is already impossible to find values of σ 1 where the amplitude shows transitions between the two states, and the threshold for which the signal merge drops to σ 1 ≈ 0.2. For γ 2 = 0.001 (=γ 1 ), although the oscillation outside the gap still is a stable solution to equations (7) to (10), very weak noise (  σ − 10 1 3 ) is able to induce the transition to the state inside the gap, which is later never abandoned. For small values of γ 2 , therefore, the resonant oscillation is virtually the only accessible stationary state under the action of noise.
In order to gain some insight on the origin of this behavior, we have measured the standard deviation of A 1 around the stationary states as a function of γ 2 . Figure 4A shows the results for a fixed noise intensity σ 1 = 0.01. Full dots corresponds to fluctuations around the state inside the gap. Deviations in the amplitude decrease as γ 2 becomes smaller, indicating that the resonant oscillation is more robust to noise for smaller damping of oscillator 2. The state outside the gap (empty dots), in contrast, has the opposite trend. In particular, the standard deviation  Fig. 2. Full and empty dots stand for initial conditions that, in the noiseless case, lead to stable states inside and outside the gap, respectively. Measurements correspond to averages over 2 × 10 5 to 5 × 10 6 time units. Dotted curves are spline interpolations, plotted as a guide to the eye. Straight segments between the two standard deviation datasets for small σ 1 have slope one. Vertical stripes correspond to the intervals of σ 1 where A 1 performs occasional transitions between the states inside and outside the gap, as in Fig. 2C. exhibits a sharp growth for γ > .
∼ 0 001 2 , which is compatible with the fact that the oscillators migrate to the resonant state even for very small noises.
In a naive picture, we expect the standard deviation of A 1 to be proportional to the noise strength, and to the square root of a typical time of relaxation toward equilibrium. These are in fact the two main ingredients that determine the size of fluctuations induced by noise on a system at equilibrium in a harmonic potential well (Ornstein-Uhlenbeck process 28 ). The lower panels of Fig. 3 make it clear that, for weak noise, the standard deviation is indeed proportional to σ 1 . The damping coefficient γ 2 , on the other hand, controls -along with all the other parameters-the eigenvalues of equations (7) to (10) at their stable fixed points, which can in turn be used to estimate the relaxation times toward equilibrium 17 . The eigenvalues relevant to this estimation have the form µ ν = −Λ ±i in,out in,out in,out where the subindices "in,out" refer to the states inside and outside the resonance gap. The respective relaxation times are given by minus the inverse of their real parts, Λ − in,out 1 . Figure 4B and C show the standard deviation of A 1 around each stable state as a function of the corresponding values of Λ − in,out 1/2 . While for the resonant state we find a linear relation, as expected on the basis of the above argument, for the state outside the gap the dependence is far from linear. This indicates that the naive picture of a response to noise controlled just by linear relaxation is not valid for the stationary state outside the resonance gap. Plausibly, nonlinear effects are involved in the process by which this state looses stability under the action of the stochastic force.
Dependence on other parameters. In the preceding paragraphs, we have analyzed the response to noise of the two-oscillator model as a function of the damping coefficient of oscillator 2 -which is one of the parameters that controls the width of the resonance gap-and of the strength of the stochastic force acting on oscillator 1. In this subsection, we briefly report on the effects of varying the remaining parameters.
In the first place, we have numerically solved equations (1) and (2) for two fixed values of σ 1 , varying the strength of the stochastic force on oscillator 2, σ 2 , with all the other parameters as in Figs 2-4 and γ 2 = 0.001. Naturally, in an experimental realization of the c-c oscillator, the noise amplitudes σ 1 and σ 2 are not independent parameters, but are related by the continuum dynamics of the beam. In our numerical simulations, however, we have the freedom of varying each of them separately. Figure 5A shows the standard deviation of the amplitude A 1 around the stable state inside the resonance gap for σ 1 = 0.01 and σ 1 = 0.001 as a function of σ 2 . In both cases, we find two well-defined regimes. For  σ σ 2 1 and below, the standard deviation does not depend on σ 2 , and its value coincides, up to small variations, with that obtained for σ 2 = 0. This indicates that, in this regime, the effect of noise is dominated by σ 1 . However, as the stochastic force on oscillator 2 grows stronger, σ σ > ∼ 10 2 1 and above, the fluctuations in the amplitude become proportional to σ 2 , indicating a prevalence with respect to the noise on oscillator 1. Note that, due to the linearity of the equation of motion for oscillator 2, its amplitude should fluctuate proportionally to σ 2 . The coordinate x 2 , in turn, enters the equation for oscillator 1 weighted by the coupling strength j. Therefore, we expect that the stochastic force on oscillator 2 contributes fluctuations proportional to the product σ j 2 to the amplitude of oscillator 1. We have already mentioned that the nonlinear coefficient β and the self-sustaining force gain g determine the form of the Duffing resonance peak but, as long as they are large enough, they do not affect the zone of the resonance gap. Therefore, they are not expected to influence the response to noise of the stable resonant oscillations. The same holds for the frequency ω 2 , which defines the position of the resonance gap, as long as it lies far enough from the end of the resonance peak. Moreover, the occurrence of internal resonance does not require that ω 2 has any special relation to ω 1 . On the other hand, the damping coefficient γ 1 controls the width of the peak along its whole length. To assess the effect of varying this coefficient, we have solved equations (1) and (2) with the same parameters as in Fig. 5A, fixing σ = .
0 01 1 and σ 2 = 0, for various values of γ 1 . In Fig. 5B we show the standard deviation of oscillator 1 in the stable state inside the gap as a function of Λ − in 1/2 (see Fig. 4B), with Λ in calculated as a function of γ 1 from the eigenvalues of equations (7) to (10). We verify that the linear dependence reported in Fig. 4C for varying γ 2 subsists upon variation of γ 1 .
Finally, we focus on the coupling strength j. In our discussion of the size of the resonance gap, we pointed out that its width is defined by the combination γ Γ = − j 2 2 . Conjecturing that this combination might also control the response to noise inside the gap, we have solved equations (1) and (2) with the same parameters as before, fixing γ 2 and varying j 2 . Figure 5C shows the standard deviation of oscillator 1 inside the gap as a function of Γ. The collapse of the results onto a single curve confirms that our conjecture is valid for these parameter sets.

Discussion
We have studied the response to stochastic forces of a system of two coupled mechanical oscillators, which model the internal resonance of a clamped-clamped solid beam. The coexistence of two stable oscillatory states -one of them representing the resonant oscillation-raises the question of their relative stability under the action of noise. Answering this question is relevant to possible applications of internal resonance, specifically, as a way of stabilizing the oscillation frequency in micromechanical nonlinear oscillators 9 .
We have shown that the relative stability of the two coexisting oscillations depends on the parameters that, in turn, control the width of the frequency gap induced by internal resonance -in particular, the damping coefficient of the resonant mode. Due to its dimensionality, it is not possible to formally describe our model as a system randomly evolving under the action of noise on a potential landscape 28 . However, from our numerical results, a qualitative picture emerges in which, as the damping coefficient of the resonant mode decreases (i.e. its quality factor grows), the resonant oscillation sharply strengthens its stability with respect to the other oscillation, in spite of the increasingly closer vicinity of an unstable state across the gap. This is schematically shown in Fig. 6, where stable and unstable oscillations are respectively pictured as the minima and maxima of a (fictitious) nonequilibrium potential. It is worth mentioning that, under experimentally controlled conditions, the damping coefficients of each mode can, to a certain extent, be varied independently from each other 29 . Here we have moreover considered variations of the coupling strength between the modes, which also influences the gap width. As a final comment, let us point out that our selection of a self-sustaining force proportional to the oscillation velocity is just but one of the possible feedbacks by which a micromechanical oscillator can be excited into stationary motion. In our case, the choice was motivated by convenience in the numerical implementation. Other experimentally tested feedbacks include self-sustaining forces with controlled maximal strength 9 , and proportional to the oscillation amplitude 19 . Although the zones of parameter space where stationary motion is stable or unstable may vary depending on the kind of excitation, internal resonance is expected to occur in all of them, so that our present conclusions should apply -at least, qualitatively-to those other realizations of the self-sustaining feedback. Figure 6. Qualitative picture, in terms of a fictitious nonequilibrium potential represented by the curves, of the relative stability of the stable oscillations inside and outside the resonance gap, when varying the damping coefficient γ 2 of oscillator 2. The two minima and the intermediate maximum stand, respectively, for stable and unstable oscillations.