Nanotube-Based NEMS: Control vs. Thermodynamic Fluctuations

Multi-scale simulations of nanotube-based nanoelectromechanical systems (NEMS) controlled by a nonuniform electric field are performed by an example of a gigahertz oscillator. Using molecular dynamics simulations, we obtain the friction coefficients and characteristics of the thermal noise associated with the relative motion of the nanotube walls. These results are used in a phenomenological one-dimensional oscillator model. The analysis based both on this model and the Fokker-Planck equation for the oscillation energy distribution function shows how thermodynamic fluctuations restrict the possibility of controlling NEMS operation for systems of small sizes. The parameters of the force for which control of the oscillator operation is possible are determined.

The gigahertz oscillator based on a double-walled carbon nanotube is one of the simplest nanotubebased NEMS. Thus, it is commonly used as a model system for studying the tribological behavior of the nanotube-based NEMS 23,24, , , , , , , , , , , , ,  25 26 27 28 29 30 31 32 33 34 35 36 37 and possible methods for controlling their operation 38,39 . The scheme, operational principles, and theory of the gigahertz oscillator based on the relative sliding of carbon nanotube walls were considered by Zheng et al. 14,15 . The oscillator was proposed to be used as a part of the device for surface profiling 40 where Q is the Q-factor corresponding to the oscillation with period .

II. MOLECULAR DYNAMICS SIMULATION OF GIGAHERTZ OSCILLATOR
To obtain the oscillator Q-factor, study the characteristics of the thermal noise, and demonstrate the possibility of controlling the motion of a functionalized nanotube wall by a nonuniform electric field, we performed MD simulations of the free and controlled operation of the ( The analysis of the free and controlled behavior of the gigahertz oscillator was performed using empirical interatomic potentials. Interaction between the inner and outer wall atoms i and j at distance was described by the Lennard-Jones 12-6 potential ij r ( ) where and are the total kinetic energies of the inner and outer walls, respectively, and are the total Brenner potentials of the walls, and the last term is the sum over all pairs of atoms i of the inner wall and atoms The value of the total van der Waals force, which retracts the telescopically extended inner wall back into the outer wall, was found to be W 1180 pN F = for the above parameters of the Lennard-Jones potential. For large-diameter nanotubes, the interwall van der Waals energy should be proportional to the overlap area between the walls. Therefore, the total van der Waals force is proportional to the nanotube diameter. The van der Waals force per unit nanotube diameter was found to be about 0.
where is the mass of the movable wall, and m υ is its center-of-mass velocity. The calculation of the = E E Q . There are two contributions to the Q-factor calculation error. On the one hand, the Q-factor depends on the oscillation amplitude. Thus, the damping of the oscillation leads to a change in the Q-factor. This effect becomes relevant at a long simulation time. On the other hand, there is also a stochastic contribution to the Q-factor calculation error related to thermodynamic fluctuations, which is high at a short simulation time and decreases with its increase. We took a simulation time of 500 ps, which minimizes the Q-factor calculation error in a single run and provides the accuracy of the Q-factor calculations within 20%.
Temperature T of the system was monitored based on the total thermal kinetic energy of the atoms in the system N where is Boltzmann's constant, B k is the number of atoms in the system. The total thermal kinetic energy was averaged over a time interval ps 1 = Δt , as this time interval must be much longer than the period of thermal vibrations (about 0.1 ps) and, on the other hand, much shorter than the period of the telescopic oscillation (about 10 ps). The temperature change over the simulation time was less than 9% at pre-heating temperatures of 50, 100, 150 and 300 K. At a pre-heating temperature of 0 K, temperature increased to 2 K within the simulation time. An estimation 24 showed that the influence of quantum effects on the dissipation rate can be neglected at temperatures above 10 K. The calculated Q-factor Q at different pre-heating temperatures is given in TABLE I. It was found that the Q-factor of the oscillator Q strongly increases with decreasing temperature (see The investigation of the tribological properties of the (5,5)@(10,10) nanotube-based oscillator has revealed that the populations of the nanotube vibrational levels are in thermal equilibrium, which justifies the applicability of the fluctuation-dissipation theorem 53,54 for this system. As the dissipation rate is proportional to the number of excited phonons, it was shown that the inverse Q-factor should linearly depend on temperature, in agreement with TABLE I. Since a certain energy exchange between the telescopic oscillation of the movable wall and the other degrees of freedom exists at the zero temperature even in the classical limit, the Q-factor remains finite at the zero pre-heating temperature.
According to the analysis performed within the framework of the fluctuation-dissipation theorem 53,54 , the relative deviation E δ of the relative energy change over a half-period of the oscillation is determined by the expression As the Q-factor of the considered gigahertz oscillator is almost inversely proportional to temperature (see TABLE I and also Ref. 24), it follows from Eq. (7) that E δ weakly depends on temperature.
However, this result is valid unless the applicability conditions of the fluctuation-dissipation theorem are not violated. One of these applicability conditions is that all degrees of freedom except for one should be in thermal equilibrium. We suppose that, at very low temperatures, the system is highly nonequilibrium. This results in the deviation from Eq. (7) and the great value of E δ for the zero pre-heating temperature (see TABLE I). Note that, according to Eq. (7), the significant fluctuations of the relative energy change over a half period of the oscillation should be attributed to the small oscillation energy, i.e. the small size of the system under consideration or the small amplitude.

E E / Δ
In the simulations of the controlled oscillations, the hydrogen functionalized inner wall was exposed to the harmonic electric field of a spherical capacitor. In this case, one more term was added to the where is the charge of the atom of the movable wall, and ( ) is electric field potential at the position of the atom i with the radius-vector i r r from the center of the spherical capacitor ( ) ( Here is the amplitude of the applied voltage, Here, = 3.1 nm is the wall length, L 0 34 R . δ = nm is the interwall distance, = 0.34 nm is the inner wall radius. As it is seen, the temperature difference between the walls is negligible. Moreover, as it was shown above, the Q-factor is inversely proportional to temperature. Therefore, oscillator operation at low temperatures requires a smaller amplitude of the control force. For instance, at the temperature of liquid helium (4.2 K), the necessary amplitude of the applied voltage is only several volts.
It should also be mentioned that, in some cases, even at relatively high amplitudes of the control force, the breakdown of the oscillation occurs (see FIG. 4c). The analysis below shows that this breakdown is induced by thermodynamic fluctuations. In fact, let us suppose that occasionally a significant negative fluctuation of the oscillation energy occurs. This means that the oscillation amplitude decreases. Since the oscillator frequency strongly depends on the oscillation amplitude 14,15 , the oscillator gets out of the resonance with the control force. This leads to a decrease of the work of the control force, which, in turn, results in a further decrease in the oscillation energy. As follows from this explanation, the stability of the oscillator operation might be improved with increasing the amplitude of the control force (above the critical value) or with decreasing the level of fluctuations in the system. It is seen that the breakdown of the oscillation induced by thermodynamic fluctuations should be expected for any strongly anharmonic oscillator.

III. PHENOMENOLOGICAL MODEL OF GIGAHERTZ OSCILLATOR
MD simulations allow the study of system behavior only at times of a few nanoseconds. To reach longer simulation times, a phenomenological model with the parameters derived from MD simulations is proposed. If there is no resonance between the telescopic oscillation and other vibrations in the nanotube (which is the case for the considered (5,5)@(10,10) nanotube), the oscillator dynamics may be roughly described by a one-dimensional equation of motion vdW fr 14 Here, x is the displacement of the movable wall; is its mass; is the interwall van der Waals force; x & is the friction force modeling energy dissipation; is the white noise representing random fluctuating forces; and is the control force. The initial conditions were As above, we considered the harmonic control force , where is control force amplitude, 0 F ω is the angular frequency of the control force, and ϕ is the initial phase shift between the control force and the velocity of the movable wall. At a relatively low oscillation energy, the dynamic friction force is proportional to the oscillator velocity 36,37 with the friction coefficient η , The friction coefficient η was chosen so as to reproduce the Q-factor of the oscillator at the initial oscillation amplitude 1 nm s = (see TABLE I). Since the phenomenological model was used here to study the qualitative behavior of the system, the dependence of the friction coefficient η on the oscillation amplitude was disregarded, i.e. the friction coefficient η was supposed to be constant. To determine the conditions of the stationary operation and study the relaxation of the oscillator to the stationary operation, the thermal noise was neglected ( ( ) 0 t ξ ≡ ).
Two approximations for the interwall van der Waals force were considered. The first one was the calculated dependence of the interwall force on the displacement between the walls for the can be assumed constant 14,15 vdW The magnitude of the interwall force was taken equal to that in the "short nanotube" case where is oscillation period corresponding to amplitude . s In this "long nanotube" case, motion equation (11) can be solved semi-analytically with the time step equal to the oscillation half-period. This leads to the decrease in the calculation time by a factor of 10 − 100 compared to that in the "short nanotube" case.
In the "long nanotube" case, the critical amplitude of the control force is determined by Eq. (1). In the "short nanotube" case, the interwall van der Waals force is not constant and has a smaller average value.
So, in the latter case the critical amplitude of the control force is less than the value given by Eq. (1).
The calculated critical amplitudes of the control force and voltage for the "short nanotube" case are presented in  . 4a and FIG. 4b) for the same nanotube within the accuracy of the Q-factor calculations.
where is the characteristic time required to reach the stationary operation mode, 16 The period of the frequency and amplitude oscillations is given by Note that in the case , the phase shift of the control force in the stationary operation mode is If the initial phase shift of the control force considerably differs from zero (ϕ~π ), the stationary operation mode is possible only in the case Here is the Dirac delta-function.
Integrating motion equation (11), we changed the value ( ) t ξ randomly with the time step ξ τ .
Therefore, the dispersion of ( ) t ξ was determined by Eq. (19) can be presented in the obvious form Let us first consider the free oscillation. The damping (or drift) of the oscillation energy is determined s QT Note that the product of the Q-factor and the oscillation period is independent of the oscillation amplitude for the considered phenomenological model with the constant friction coefficient η (see Eq. (13)).
Based on the fluctuation-dissipation theorem, it was shown that the thermal noise leads to the diffusion of the oscillation energy with the diffusion coefficient where is a numerical factor found below using the phenomenological model. 1 α ≈ Thus, the Fokker-Planck equation takes the form To solve Eq. (26), let us introduce the variables ( ) is the initial oscillation energy.
ζ χ = For these variables and t − , Eq. (26) is reduced to In the case 0 E ζ << , which means that the width of the oscillation energy distribution is small compared to the current average oscillation energy, one gets a simple one-dimensional diffusion From these considerations it can be seen that the solution of Eq. (26) for the initial condition can be approximated as Note that the diffusion coefficient still linearly depends on the oscillation energy.
We numerically calculated the energy distribution function for damping oscillations in the "long nanotube" case with 1 fs ( ) Therefore, the distribution function in the stationary operation mode, which is established within the time of about is given by Substituting expressions (27) for (with 0 75 .
It can be seen that the width of the oscillation energy distribution depends on the ratio of the stationary oscillation energy s E to the thermal kinetic energy and is independent of the Q-factor.
where the pre-exponential factor is chosen so as to provide s Substituting Eq. (41) into Eq. (40), the lifetime of the stationary oscillation is roughly given by c 0 2 0c  FIG. 6a,b. This is because, in the both cases, we consider the maximum relative c difference between the oscillation frequency and the frequency of the control force for which the oscillation is still sustained. Note that estimates of the oscillator stability regions similar to Eqs. (41) and (43) should be valid for any anharmonic oscillator with a strong dependence of the oscillation period on the oscillation energy ~. Therefore, the lifetime of the stationary operation mode should be finite for any strongly anharmonic oscillator.
the nanotube length at which the resonance between these modes is possible can be estimated as So for the (5,5)@(10,10) nanotube, the maximum length up to which the one-dimensional model is applicable is about . For longer nanotubes, coupling of the telescopic oscillation with the transverse acoustic modes should be taken into account. max 140 nm L ≈ In the present paper we only considered the case when there is no resonance between the telescopic oscillation and other nanotube vibrational modes. As it was shown by MD simulations, in the case of the resonance, not only the Q-factor of the oscillator decreases drastically but the level of thermodynamic fluctuations is also a few times greater than that under the conditions when the fluctuation-dissipation relation is satisfied (see Eq. (7)). Therefore, according to Eq. (42), the lifetime of the stationary operation mode in the case of the resonance should be significantly smaller compared to that for the case when the resonance is absent. So the operation of the oscillator in the case when there is a resonance with any nanotube vibrational mode should be avoided.
As for condition (2), the phenomenological parameters for the short nanotube of 3 nm length considered above were extracted from the MD simulations and validity of the model was based on matching to the MD results. For longer nanotubes, the phenomenological parameters of the model can be determined accurately on the basis of experimental measurements. Note that as the fluctuationdissipation theorem, which relates the thermal noise to the energy dissipation rate, is valid for these nanotubes, it is sufficient to measure only the Q-factor of the oscillator.
On the basis of the obtained results, let us discuss the restrictions imposed by thermodynamic fluctuations on sizes of the nanotube-based NEMS for which control over the NEMS operation is possible. From Eq. (42) and FIG. 10, it is seen that the lifetime of the stationary operation mode s τ can be increased by decreasing the level of fluctuations δ . According to Eq. (21), δ is proportional to the square root of the Q-factor divided by the oscillation energy. Therefore, a decrease of δ can be achieved by an increase of the oscillation amplitude. However, for the oscillation amplitude greater than about 30% of the inner wall length, the dissipation rate strongly increases (and, consequently, the Qfactor strongly decreases) with an increase of the oscillation amplitude due to the excitation of lowfrequency vibrational modes 26,35 . Thus, for an oscillator with a certain length the minimum level of fluctuations 35 δ and, therefore, the maximum lifetime of the stationary operation mode s τ correspond to some optimal oscillation amplitude about 30% of the oscillator length. From the MD simulations it follows that if the ratio of the initial telescopic extension of the inner wall to the oscillator length is maintained, the Q-factor is independent of the oscillator length (in the region , see Eq. (45)).
In the same case, the oscillation energy is proportional to the oscillator length. Hence, one gets the minimum value of max L L << δ for oscillators with the optimal oscillation amplitude to be inversely proportional to the square root of the oscillator length. In accordance with Eq. (42), this means that the lifetime of the stationary operation mode s τ for such oscillators increases exponentially with increasing their length.

V. CONCLUSION
In summary, we performed molecular dynamics (MD) simulations of the controlled operation of the nanotube-based gigahertz oscillator. The feasibility of control over a movable nanotube wall which is functionalized so that it has an electric dipole moment by a nonuniform electric field was demonstrated.
The MD simulations were used to obtain the oscillator Q-factor and the characteristics of the thermal noise. To study the possibility of the stationary operation mode (operation mode with a constant frequency) at a simulation time of 1 ms, a phenomenological one-dimensional model with the parameters derived from MD simulations was proposed. Using this model, the control force parameters