Nanotube-based nanoelectromechanical systems: Control versus thermodynamic fluctuations

Multiscale 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.


I. INTRODUCTION
The ability of the free relative sliding and rotation of carbon nanotube walls 1-3 and their excellent "wearproof" characteristics 3 allow the use of carbon nanotube walls as movable elements in nanoelectromechanical systems ͑NEMS͒. A number of devices offering great promise for applications in NEMS and based on the use of the relative motion of carbon nanotube walls have been proposed recently. These devices include rotational 4,5 and plain 2 nanobearings, nanogears, 6 electromechanical nanoswitchs, 7 nanoactuators, 8,9 Brownian motors, 10 nanobolt-nanonut pairs, [11][12][13] and gigahertz oscillators. 14,15 Furthermore, nanomotors based on the relative rotation of carbon nanotube walls [16][17][18][19] and memory cells based on the relative sliding of carbon nanotube walls 20 have been implemented.
The crucial issue in nanotechnology is the actuation of the NEMS components in a controllable way. A new method for controlling the motion of NEMS based on carbon nanotubes 21 was proposed lately. Namely, the wall of a nanotube has an electric dipole moment if electron donors or/and acceptors are adsorbed at the ends of the wall. The motion of such a functionalized wall with an electric dipole moment can be controlled by a nonuniform electric field. Here, molecular dynamics ͑MD͒ simulations are performed to demonstrate the feasibility of this control method.
As compared to microelectromechanical systems, the principal feature of NEMS related to a small number of atoms in the system is the significance of thermodynamic fluctuations in NEMS. Here, we perform multiscale simulations to study the influence of these fluctuations on the NEMS operation by the example of a nanotube-based gigahertz oscillator. We believe that the approach developed and the re-sults of this study can be also useful for other types of NEMS, such as artificial molecular rotors 22 and so on.
The gigahertz oscillator based on a double-walled carbon nanotube is one of the simplest nanotube-based NEMS. Thus, it is commonly used as a model system for studying the tribological behavior of the nanotube-based NEMS ͑Refs. 23-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 Upon the telescopic extension of the inner wall outside the outer wall, the van der Waals force F W turns the inner wall back into the outer wall, thereby makes this NEMS oscillate. However, such oscillations are anharmonic and dissipative with the Q-factor Q Ϸ 10-1000 ͑Refs. 24, 26-28, and 32͒ ͑Q-factor of the system is the ratio of the total oscillation energy to the energy loss per one oscillation period͒. The frequency of the damping oscillations increases with decreasing the oscillation amplitude. 14 Thus, this frequency increases with time. 27,30 Consequently, to provide the stationary operation of the gigahertz oscillator, that is, to keep its frequency constant, it is necessary to compensate the energy dissipation by the work of an external force. To obtain the critical value of this force ͑i.e., the minimum value at which the work of the external force can compensate the energy dissipation͒, the energy balance in the controlled operation of the gigahertz oscillator was analyzed. 21 It was shown that the critical amplitude F 0c of the harmonic control force F͑t͒ = F 0 cos͑2t / T s ͒ applied to the movable inner wall attains a minimum value for the oscillator with walls of equal lengths and is given by the expression where Q is the Q-factor corresponding to the oscillation with period T s . Here we study the characteristics of the actuation of a functionalized nanotube wall by a nonuniform electric field through multiscale simulations of the controlled operation of the ͑5,5͒@͑10,10͒ nanotube-based gigahertz oscillator. MD simulations are performed to demonstrate the feasibility of this control method. To estimate the critical amplitude of the control force, we calculate the Q-factor of the gigahertz oscillator using MD simulations of the damping oscillations. However, the MD technique does not allow studying systems containing a large number of atoms within a long simulation time. So, the operation of the gigahertz oscillator is analyzed within the framework of a phenomenological onedimensional model with the parameters derived from MD simulations. This analysis is used to investigate the system behavior in the course of its relaxation to the stationary operation and to determine the conditions under which the stationary operation of this NEMS is possible.
By now, thermodynamic fluctuations in NEMS were considered and investigated using MD simulations only to demonstrate the possibility of the relative diffusion of the NEMS components 41 and regarding the equilibrium rotational 42 and translational 37 dynamics of walls in carbon nanotubes. As for the nanotube-based NEMS, such diffusion can be used in Brownian motors. 10 However, diffusion 43 or displacement 9 of the NEMS components due to thermodynamic fluctuations can disturb the NEMS operation. 9,43 For example, it was shown that thermodynamic fluctuations restrict the minimum sizes of the electromechanical nanothermometer based on the interaction of nanotube walls 43 and the electromechanical nanorelay based on the relative motion of nanotube walls. 9 The MD simulations performed here revealed substantial fluctuations in the nanotube-based gigahertz oscillator. The analysis of the motion equation of this NEMS indicates the critical influence of these fluctuations on the possibility of controlling its operation. In this way, the principal restrictions imposed by thermodynamic fluctuations on the possibility of controlling the NEMS operation are studied within the framework of the phenomenological one-dimensional model and analyzed on the basis of the Fokker-Planck equation for the oscillation energy distribution function.
The paper is organized in the following way. The results of the MD simulations of the operation of the gigahertz oscillator are given in Sec. II. In Sec. III, we examine the influence of the characteristics of the control force on the operation of the gigahertz oscillator within the framework of the phenomenological model with the parameters derived from the MD simulations. The effect of thermodynamic fluctuations on the NEMS operation is analyzed in Sec. IV. Our conclusions are summarized in Sec. V.

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 con-trolling the motion of a functionalized nanotube wall by a nonuniform electric field, we performed MD simulations of the free and controlled operation of the ͑5,5͒@͑10,10͒ nanotube-based oscillator ͑see Fig. 1͒. For this double-walled nanotube, there is no resonance between the telescopic oscillation and other nanotube vibrations, and the populations of the vibrational levels are in thermal equilibrium. 24 Both nanotube walls were taken equal to 3.1 nm in length. One end of the inner wall was capped and the other end was open and terminated with hydrogen atoms ͑see Fig. 1͒. Both ends of the outer wall were open and not functionalized. The charge distribution in the inner wall was calculated by the semiempirical method of molecular orbitals with the PM3 parameterization of the Hamiltonian. 44 The adequacy of the PM3 parameterization in the case under consideration was demonstrated by calculations of bond lengths in the C 60 fullerene with the symmetry I h : the calculated bond lengths agree with their experimental counterparts within 10 −4 nm. 45 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 r ij was described by the Lennard-Jones 12-6 potential with the parameters CC = 3.73 meV, CC = 3.40 Å and CH = 0.65 meV, CH = 2.59 Å for the carbon-carbon and carbon-hydrogen interactions, respectively, taken from the AMBER database 46 for aromatic carbon and hydrogen bonded to aromatic carbon. The parameters provide a consistent description of the pairwise carbon-carbon and carbon-hydrogen interactions. The cut-off distance of the Lennard-Jones potential was taken equal to 12 Å. The covalent carbon-carbon and carbon-hydrogen interactions inside the walls were described by the empirical Brenner potential, 47 which was shown to correctly reproduce the vibrational spectra of carbon nanotubes. 48 So the Hamiltonian of the oscillator based on a double-walled nanotube was given by where K ͑out͒ and K ͑in͒ are the total kinetic energies of the inner and outer walls, respectively, V B ͑out͒ and V B ͑in͒ 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 j of the outer wall. The value of the total van der Waals force, 24 which retracts the telescopically extended inner wall back into the outer wall, was found to be F W = 1180 pN for the above parameters of the Lennard-Jones potential. For largediameter 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.1-0.2 N/m in the experiments of Cumings and Zettl, 2 and Kis et al. 3 We obtained the value of the van der Waals force per unit nanotube diameter of about 0.3 N/m. However, one should take into account that the measurements 2,3 were performed for largediameter multiwalled nanotubes, whereas a nonlinear dependence of the van der Waals force on the nanotube diameter should be expected for small-diameter nanotubes. The barrier to relative sliding of the walls was calculated to be 0.008 meV per atom of the outer wall, in agreement with the recent result obtained for the ͑5,5͒@͑10,10͒ nanotube through local-density-approximation calculations. 49 The Lennard-Jones potential was also shown to provide the interlayer interaction energy in graphite of about 62 meV/atom, which is consistent with the experimental value 52Ϯ 5 meV/ atom obtained recently from the experiments on the thermal desorption of polyaromatic hydrocarbons from a graphitic surface 50 but greater than the values reported earlier ͑see Ref. 51 and references therein͒.
An in-house MD-KMC code 52 was implemented. The code used the velocity Verlet algorithm and neighbor lists to improve the computing performance. The time step was 0.2 fs, which is about two orders of magnitude shorter than the period of thermal vibrations of hydrogen atoms. The initial configuration of the nanotube was optimized at zero temperature. The initial velocities of the atoms were distributed according to the Maxwell-Boltzmann distribution with a doubled preheating temperature 2T 0 . This provided that during the MD simulation, the initial distribution of velocities reached the Maxwell-Boltzmann temperature corresponding to the preheating temperature T 0 in less than 0.02 ps. To start the oscillation, the inner wall was pulled out along the nanotube axis at the distance s =1 nm ͑about 30% of its length͒ and released with the zero center-of-mass velocity. The outer wall was fixed at three atoms. The relative fluctuations of the total energy of the system caused by numerical errors were less than 0.3% of the interwall van der Waals energy.
To get the temperature dependence of the oscillator Q-factor, we performed microcanonical MD simulations of free oscillations ͑see Fig. 2͒ at preheating temperatures of 0, 50, 100, 150, and 300 K. The oscillation energy E is given by the sum of the center-of-mass kinetic energy of the movable wall and the excess of the interwall van der Waals energy V͑x͒ associated with the displacement x of the movable wall from the position corresponding to the minimum of the interwall van der Waals energy ͓V͑0͒ =0͔ where m is the mass of the movable wall and is its centerof-mass velocity. The calculation of the oscillation energy was performed at the moments when the movable wall crossed the position x = 0 corresponding to the minimum of the interwall van der Waals energy. So the oscillation energy was found at these moments as E = m 2 / 2. Since the oscillation energy loss ⌬E over a half period of the oscillation is much smaller than the oscillation energy E ͑⌬E Ӷ E͒, the Q-factor was estimated using the ratio ⌬E / E averaged over a single simulation run

͑5͒
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 3 2 where k B is Boltzmann's constant and N is the number of atoms in the system. The total thermal kinetic energy was averaged over a time interval ⌬t = 1 ps, 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 preheating temperatures of 50, 100, 150, and 300 K. At a preheating 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 preheating temperatures is given in Table I.
To estimate the level of fluctuations in the system, we also calculated the root-mean-square deviation of the relative energy change ⌬E / E over every half period of the oscillation. Significant fluctuations of the quantity ⌬E / E were revealed ͑see Fig. 3͒. The calculated relative root-mean-square deviation ␦ E of ⌬E / E at different temperatures is given in Table I.
It was found that the Q-factor of the oscillator Q strongly increases with decreasing temperature ͑see Table I͒, in agreement with the other papers. 24,[34][35][36][37] Nevertheless, the Q-factor does not go to infinity at the zero preheating temperature.
The relative deviation ␦ E of ⌬E / E weakly depends on temperature in the temperature range of 50-300 K. However, the considerable increase in ␦ E is observed for the zero preheating temperature. Explanations of these results are presented below.
The investigation of the tribological properties of the ͑5,5͒@͑10,10͒ nanotube-based oscillator 24 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, 24 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 preheating temperature.
According to the analysis 24 performed within the framework of the fluctuation-dissipation theorem, 53,54 the relative deviation ␦ E of the relative energy change ⌬E / E over a halfperiod 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 preheating temperature ͑see Table I͒. Note that, according to Eq. ͑7͒, the significant fluctuations of the relative energy change ⌬E / E 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.
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 Hamiltonian ͑3͒ of the system where q i is the charge of the atom i of the movable wall and ͑r ជ i ͒ is electric field potential at the position of the atom i with the radius vector r ជ i from the center of the spherical capacitor Here U 0 is the amplitude of the applied voltage, is the angular frequency of the electric field, is the initial phase shift between the electric field and the velocity of the movable wall, R 1 and R 2 are the radii of the inner and outer plates of the capacitor, respectively. The angular frequency of the electric field equaled the oscillation angular frequency 0 corresponding to the initial oscillation amplitude s =1 nm ͑ = 0 ͒. The initial phase shift between the electric field and the inner wall velocity equaled zero ͑ =0͒. The radii of the plates of the spherical capacitor were taken equal to R 1 = 100 nm and R 2 = 110 nm. Furthermore, in these simulations, the temperature of the outer wall was maintained by periodically rescaling atomic velocities every 0.1 ps ͑the Berendsen thermostat 55 ͒. Though only the outer wall was kept in contact with the thermostat, our estimation demonstrates that the temperatures of both walls are almost equal. In fact, the dissipation rate at a temperature of 300 K was found to be W Ϸ 10 −9 W. The thermal conductivity between the nanotube walls can be assumed equal to that between graphite layers 56 Ќ Ϸ 5 W/ ͑m K͒. Under stationary conditions, the dissipation rate should be equal the rate of heat transfer between the walls. Supposing that the heat flux is identical through all coaxial cylindrical surfaces between the inner and outer walls, one gets the temperature difference between the nanotube walls Here, L = 3.1 nm is the wall length, ␦R = 0.34 nm is the interwall distance, R in = 0.34 nm is the inner wall radius. As it is seen, the temperature difference between the walls is negligible.
The results of the MD simulations of the controlled oscillations are presented in Fig. 4. The amplitude of the applied voltage was 60.5 and 61.4 V at a temperature of 300 K and 10.9 V at 50 K. The oscillations shown in Fig. 4 correspond to the control force amplitude greater than the critical value for the considered oscillator. In this case, the variation in the oscillation amplitude is observed, in agreement with the results obtained below using the phenomenological onedimensional model.
Recently, a nanomotor based on the relative rotation of nanotube walls in a multiwalled carbon nanotube was implemented. 17 In this experiment, the voltage applied between the nanotube wall and the control electrode reached 100 V. This value is of the order of the voltage amplitude needed to sustain the oscillation at a constant amplitude in the performed MD simulations. Thus, the MD simulations have demonstrated that the method under consideration can be used to control the motion of the nanotube-based NEMS.
The amplitude of the applied voltage needed to sustain the oscillation at a constant amplitude can be decreased if the dipole moment of the inner wall is increased. This can be achieved, for example, through the adsorption of hydrogen and fluorine atoms at opposite open ends of the nanotube. 21 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. 4͑c͔͒. 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 in 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 24 ͔, the oscillator dynamics may be roughly described by a one-dimensional equation of motion Here, x is the displacement of the movable wall; m is its mass; F vdW ͑x͒ is the interwall van der Waals force; F fr ͑ẋ͒ is the friction force modeling energy dissipation; ͑t͒ is the white noise representing random fluctuating forces; and F͑t͒ is the control force. The initial conditions were x͑0͒ = s = 1 nm, ẋ͑0͒ =0. As above, we considered the harmonic control force F͑t͒ = F 0 cos͑t + ͒, where F 0 is control force amplitude, 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 , F fr ͓ẋ͑t͔͒ =−ẋ͑t͒. The friction coefficient was chosen so as to reproduce the Q-factor of the oscillator at the initial oscillation amplitude s =1 nm ͑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 ͓͑t͒ϵ0͔.
Two approximations for the interwall van der Waals force F vdW ͑x͒ were considered. The first one was the calculated dependence of the interwall force on the displacement between the walls for the ͑5,5͒@͑10,10͒ nanotube 3.1 nm in length, which was used in our MD simulations ͑see Fig. 5͒. To calculate this dependence, the walls were separately relaxed and the inner wall was then rigidly shifted along the wall common axis. The interwall force was found as the partial derivative of the total interwall van der Waals energy with respect to the inner wall displacement. Since the nanotube under consideration is rather short, the region in which the interwall van der Waals force increases from zero to the nearly constant value F W Ϸ 1180 pN ͑x Ͻ 0.5 nm͒ is comparable to the oscillation amplitude ͑see Fig. 5͒. Thus, Eq. ͑11͒ with this approximation of the interwall van der Waals force represents a "short nanotube" case of the proposed model. Equation ͑11͒ was solved numerically with a time step of 1 fs and the friction coefficient for the given temperature was fitted so that the model reproduced the values of the Q-factor obtained in the MD simulations.
For a long nanotube, one can neglect the short region x Ͻ 0.5 nm in the dependence of the interwall force on the displacement between the walls, where the force increases from zero to the nearly constant value. Thus, the interwall van der Waals force F vdW ͑x͒ can be assumed constant 14,15 The magnitude of the interwall force F W was taken equal to that in the short nanotube case F W = 1180 pN ͑see Fig. 5͒. Approximation ͑12͒ of the interwall van der Waals force provides the following expression for the friction coefficient: is oscillation period corresponding to amplitude s. In this "long nanotube" case, motion equation ͑11͒ can be solved semianalytically 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 Table I. The critical voltage amplitudes agree with the results of the MD simulations ͓see Figs. 4͑a͒ and 4͑b͔͒ for the same nanotube within the accuracy of the Q-factor calculations.
Let us consider the conditions of switching the control force for the free oscillation with an amplitude s = 1 nm. The possibility of the stationary operation mode for the given F 0 / F 0c ratio is determined by the initial phase shift between the control force and relative velocity of the walls and the angular frequency of the control force. We examined the ranges of normalized parameters of the control force F 0 / F 0c , / 0 ͑where 0 is the angular frequency corresponding to the oscillation with the amplitude s =1 nm͒ and / where the stationary mode is possible for both cases of the proposed model. As it can be seen in Figs. 6 and 7, on increasing the amplitude of the control force F 0 or decreasing the Q-factor, the ranges of and which correspond to the stationary operation mode become wider. At the high amplitude of the control force F 0 = 10-100F 0c , switching of the control force is reasonable almost at an arbitrary moment ͑see Figs. 6 and 7͒.
In both short nanotube and long nanotube cases of the model we obtained similar results ͑see Fig. 6͒. So, it was demonstrated that the choice of the potential for the model has negligible influence on the behavior of the considered oscillator and the simpler long nanotube case ͑disregarding the shape of the dependence of the interwall van der Waals energy on the displacement of the movable wall͒ can be used for qualitative studies of the oscillator operation for nanotubes longer than 3 nm and oscillation amplitudes above 1 nm. This statement is also confirmed by the results obtained in Sec. IV with account of the thermodynamic fluctuations.
The relaxation to the stationary operation mode was studied for = 0 and = 0 . If the amplitude of the control force exceeds the critical value, the stationary operation mode is attained in 10-100 ns. During the relaxation to the stationary operation mode, the amplitude and frequency of the NEMS oscillate and tend to their stationary values. The frequency of the gigahertz oscillator approaches the frequency of the control force. The time dependence of the frequency of the gigahertz oscillator can be approximated by FIG. 6. Calculated parameters of the control force F 0 / F 0c , / 0 , and ͑shown in black͒ for which the stationary operation mode is possible at Q = 500. ͓͑a͒, ͑c͒, and ͑e͔͒ Short nanotube case and ͓͑b͒, ͑d͒, and ͑f͔͒ long nanotube case of the phenomenological model. ͓͑a͒ and ͑b͔͒ =0; ͓͑c͒ and ͑d͔͒ / 0 = 1; and ͓͑e͒ and ͑f͔͒ F 0 / F 0c =2.
where rel is the characteristic time required to reach the stationary operation mode, ␦ is the parameter characterizing the magnitude of the frequency oscillations, T osc is the period of the frequency and amplitude oscillations, and is the fitting parameter. A similar time dependence can be written for the amplitude of the gigahertz oscillator. Note that the variation in the oscillation amplitude revealed by the phenomenological model conforms to the results obtained in the MD simulations ͓see Figs. 4͑a͒ and 4͑b͔͒. The time parameters characterizing relaxation to the stationary operation mode as functions of the Q-factor and the amplitude of the control force are shown in Figs. 8 and 9.
Note that the period of the frequency and amplitude oscillations exceeds the oscillation period of the gigahertz oscillator by several orders of magnitude. As it is seen in Fig. 8, an increase in the Q-factor leads to an increase in the characteristic time rel ϰ Q 1.00Ϯ0.01 and the period T osc ϰ Q 0.50Ϯ0.01 , and also in a decrease in the magnitude of the frequency oscillations ␦ ϰ Q −0.50Ϯ0.01 . Near the critical value of the control force amplitude F 0 = F 0c ͓see Eq. ͑1͔͒, the parameter ␦ characterizing the magnitude of the frequency oscillations and the period T osc of the frequency and amplitude oscillations as functions of ⌬F = F 0 − F 0c are described by "universal" dependences: ␦ ϰ⌬F 0.75Ϯ0.01 and T osc ϰ⌬F −0.25Ϯ0.03 ͑see Fig.   9͒. In the stationary operation mode, the frequency and amplitude of the gigahertz oscillator and also the phase shift between the control force and relative velocity of the walls are constant.
The obtained dependences of ␦ and T osc on the Q-factor and the control force amplitude can be explained in the following way. In the stationary oscillation mode, the phase shift 0 between the control force and the velocity of the movable wall is given by cos 0 = F 0c / F 0 . So if ⌬F / F 0c Ӷ 1, 0 Ϸ ͱ 2⌬F / F 0c . If the initial phase shift of the control force equals zero ͑ =0͒, the number of the oscillation periods needed to reach the phase shift 0 is N = 0 / ͑2␦͒. Dur-  ing this time, the energy perturbation ␦E =2F W s␦ / is compensated by the excessive work of the control force N⌬Fs. So, using Eq. ͑1͒ for F W , one gets The period of the frequency and amplitude oscillations is given by Note that in the case ⌬F / F 0c ӷ 1, the phase shift of the control force in the stationary operation mode is 0 Ϸ / 2, and If the initial phase shift of the control force considerably differs from zero ͑ ϳ ͒, the stationary operation mode is possible only in the case ⌬F / F 0c ӷ 1, and the dependences in Eq. ͑18͒ are still valid. Note that estimates in Eqs. ͑16͒-͑18͒ not only give the same qualitative dependences as presented in Figs. 8 and 9 but also the numerical factors with an accuracy of 10-50 %.

IV. THERMODYNAMIC FLUCTUATIONS IN GIGAHERTZ OSCILLATOR
As it was shown in Sec. II using the MD simulations, thermodynamic fluctuations lead to the breakdown of the stationary oscillation. To examine the effect of fluctuations with the help of the one-dimensional model, one should introduce noise ͑t͒ into motion equation ͑11͒. For the ͑5,5͒@͑10,10͒ nanotube-based oscillator, the relative deviation ␦ E of the relative energy change ⌬E / E over the halfperiod of the oscillation calculated through the MD simulations was shown to be in good agreement with the analytical results obtained on the basis of the fluctuation-dissipation theorem. 24 So, thermal noise at temperature T should be a Gaussian white noise of zero mean ͗͑t͒͘ = 0 which satisfies the fluctuation-dissipation relation 53,54 ͗͑t͒͑t − ͒͘ =2k B T␦͑͒. 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 can be presented in the obvious form where T s is the period of the oscillation with the amplitude s, ͉͗F fr ͉͘ s =4s / T s is the average friction force for the oscillation with the amplitude s, and ␦ = ͱ T s k B T / ͑4s 2 ͒ is the coefficient characterizing the ratio of the noise to the friction force. In the long nanotube case, as it follows from expression ͑13͒ for the friction coefficient and expression ͑14͒ for the oscillation period, ␦ is given by the equation Comparing Eqs. ͑7͒ and ͑21͒, one can see that ␦ / ␦ E Ϸ 1. The point is that ␦ is determined by the ratio of the noise to the friction force averaged over time, while ␦ E is given by the ratio of the same quantities averaged over the displacements of the movable wall. So, ␦ and ␦ E should differ only by a numerical factor of an order of 1.
If the fluctuation-dissipation relation is not satisfied, one can suppose that the noise ͑t͒ is white and introduce it according to Eq. ͑20͒. However, in this case, the dispersion of ͑t͒ is not related to the Q factor and temperature and can be extracted, for example, from the values of ␦ E obtained by MD simulations. 24 In our simulations, we changed the value ͑t͒ randomly with the time step equal to the simulation time step 1 fs in the numerical calculations for the long nanotube and short nanotube cases and equal to the oscillation half-period in the semianalytical calculations for the long nanotube case. The dispersion of ͑t͒ was assumed to be constant. We characterize this dispersion by the value of ␦ in Eq. ͑20͒ corresponding to the initial oscillation amplitude 1 nm. Note that the values of ␦ were taken of the order of ␦ E calculated through the MD simulations ͑see Table I͒. The breakdown of the stationary mode was defined as the moment when the oscillation amplitude becomes less than 0.4 nm. The lifetime of the stationary mode s ͑i.e., the time to the breakdown of the controlled oscillation͒ for given F 0 / F 0c and ␦ was averaged over 100-400 numerical solutions of Eq. ͑11͒.
In both short nanotube and long nanotube cases, the model predicts that the fluctuations actually lead to the breakdown of oscillations, i.e., the average lifetime of the stationary mode is finite. The calculations with the onedimensional model have shown that the lifetime s increases with decreasing ␦ ͑see Fig. 10͒ and approaches infinity at ␦ → 0. The lifetime of the stationary mode can also be increased with increasing the amplitude of the control force ͑see Fig. 10͒. More specifically, the lifetime almost exponentially depends on the ratio F 0 / F 0c and, for ␦ = 1.5, becomes longer than 1 ms at F 0 / F 0c Ͼ 10. The explanation of the dependences obtained is given below. Note that the model predicts the average lifetime of the stationary mode to be about 1 ns at the amplitude of the control force slightly above the critical value, in good agreement with the breakdown of the stationary mode observed in the MD simulations ͓see Fig.  4͑c͔͒. The curves obtained in the short nanotube and long nanotube cases of the model are very close to each other, thus proving once again the applicability of the long nanotube case of the model ͑see Fig. 10͒. The semianalytical approach with the time step equal to the oscillation halfperiod provides shorter lifetimes compared to the calculations with the time step =1 fs ͑see Fig. 10͒. However, the qualitative dependences of the lifetime on F 0 / F 0c and ␦ are identical.
To explain the calculated dependence of the lifetime of the stationary operation mode of the gigahertz oscillator on F 0 / F 0c and ␦, we performed the analysis of the oscillation energy distribution function in the long nanotube case. Our analysis was based on the following assumptions: ͑1͒ the fluctuation-dissipation relation is satisfied and ͑2͒ friction coefficient is constant and independent of the oscillation amplitude. However, we believe that the qualitative conclusions made in the analysis below are valid for wider conditions of NEMS operation. The oscillation energy distribution function f͑E , t͒ characterizes the probability for the oscillator to have oscillation energy E at time instance t. The evolution of the oscillation energy distribution function in the gigahertz oscillator is determined by the Fokker-Planck equation Let us first consider the free oscillation. The damping ͑or drift͒ of the oscillation energy is determined as Note that the product of the Q-factor and the oscillation period QT s 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 24 that the thermal noise leads to the diffusion of the oscillation energy with the diffusion coefficient 24 where ␣ Ϸ 1 is a numerical factor found below using the phenomenological model. Thus, the Fokker-Planck equation takes the form To solve Eq. ͑26͒, let us introduce the variables where E 0 is the initial oscillation energy. For these variables and F͑ , ͒ = f͑E , t͒exp͑−bt͒, Eq. ͑26͒ is reduced to

͑29͒
In the case ͉͉ Ӷ E 0 , which means that the width of the oscillation energy distribution is small compared to the current average oscillation energy, one gets a simple onedimensional diffusion equation From these considerations it can be seen that the solution of Eq. ͑26͒ for the initial condition f͑E ,0͒ = ␦͑E − E 0 ͒ can be approximated as

f͑E,t͒
In the case of the small damping ͑t Ӷ T s Q͒, equation ͑31͒ takes the simple form For the considered phenomenological model with the constant dispersion of the noise ͗ 2 ͘, we used Eq. ͑19͒ to express the diffusion coefficient D and the parameter a in terms of the parameters of the model 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 as a function of time ͑see Fig. 11͒. To fit the energy distribution function, we estimated numerical factor ␣ to be ␣ = 0.75 in Eqs. ͑25͒, ͑27͒, and ͑33͒. As is seen in Fig. 11, solution ͑31͒ for this value of ␣ is in good agreement with the results of the numerical calculations.
In the stationary operation mode, the oscillation energy drift can be determined as the rate of energy relaxation to the stationary value E s . As it was shown above by the MD simulations and calculations with the phenomenological model, if the oscillation energy E is given a small perturbation, energy oscillations occur with the period T osc ϰ T s ͱ Q, which is considerably greater than the oscillation period T s ͓see Fig.  8͑c͔͒. Nevertheless, since the period of the energy oscillations T osc is much shorter than the characteristic time rel required for the reversion to the stationary operation mode, we believe that one can neglect the energy oscillations and approximate the drift term in the Fokker-Planck Eq. ͑22͒ by the relaxation rate of the absolute value of the oscillation energy deviation from the stationary value ͉E − E s ͉ averaged over the time T osc . The time of relaxation rel to the stationary operation mode was shown to linearly depend on the Q-factor rel Ϸ T s Q ͓see Fig. 8͑a͔͒. Thus, the drift term can be approximated as So, the Fokker-Planck equation for the energy distribution function in the stationary operation mode takes the form Therefore, the distribution function in the stationary operation mode, which is established within the time of about rel is given by where C is a constant determined by the normalization of f͑E͒. For ͉E − E s ͉ Ӷ E s , one gets Substituting expressions ͑27͒ for a ͑with ␣ = 0.75͒ and b, one gets It can be seen that the width of the oscillation energy distribution depends on the ratio of the stationary oscillation energy E s to the thermal kinetic energy k B T and is independent of the Q-factor.
In terms of the parameters of the considered phenomenological model with the constant noise dispersion, the oscillation energy distribution function in the stationary operation mode is obtained by substituting Eq. ͑27͒ for b and Eq. ͑33͒ for a into Eq. ͑37͒ Using the phenomenological model, we numerically calculated the oscillation energy distribution function in the stationary operation mode ͑see Fig. 12͒ averaged over 50 simulations for each set of parameters: Q-factor Q, amplitude of the control force F 0 , and the ratio ␦ of the noise to the friction force. It can be seen in Fig. 12 that there is a critical value of oscillation energy E c which separates two operation modes of the gigahertz oscillator: the damping oscillation for E Ͻ E c and the stationary operation mode where operation can be controlled at E Ͼ E c . The critical value E c can be found as the point where the derivative of the oscillation energy distribution function has a jump. The results of the numerical calculations are in reasonable agreement with Eq.
͑39͒ for E Ͼ E c ͑see Fig. 12͒. So neglecting the energy oscillations in the estimation of the drift term in the Fokker-Planck equation is adequate for the qualitative description of the system behavior.
The probability that E takes some critical value E c , which results in the breakdown of the stationary oscillation, is proportional to f͑E c ͒. Therefore, the lifetime of the stationary oscillation is determined by where the pre-exponential factor is chosen so as to provide s Ϸ T s Q for F 0 / F 0c =1. The critical value E c can be estimated from the following considerations. For the work of the control force to compensate the deviation of the oscillation energy E − E s and thus to stabilize the oscillation, at least N Ϸ͑E − E s ͒ / ͓2͑F 0 − F 0c ͒s͔ periods of the oscillation are required. However, within this time, a phase shift N 0 ⌬T s between the control force and the oscillation is accumulated, where ⌬T s = T s ͑E − E s ͒ / ͑2E s ͒ is the change in the oscillation period due to the fluctuation of the oscillation amplitude, and the actual work of the control force is less than 2N͑F 0 − F 0c ͒s. If the accumulated phase shift is significant N 0 ⌬T s Ն c , the actual work of the control force within the considered N periods is too small to compensate the dissipation. In this case, the damping of the oscillation should occur. Substituting expressions for N and ⌬T s into the equation N 0 ⌬T s Ϸ c , one gets that the critical value E c of the oscillation energy is approximately determined by Substituting Eq. ͑41͒ into Eq. ͑40͒, the lifetime of the stationary oscillation is roughly given by For F 0 / F 0c ӷ 1, the phase shift 0 between the control force and the velocity of the movable wall in the stationary operation mode tends to / 2; therefore, the critical phase shift c tends to . The obtained expression ͑42͒ for the lifetime of the stationary operation mode is in qualitative agreement with the calculated dependence of the lifetime on the control force amplitude ͑see Fig. 10͒.
On the basis of Eq. ͑42͒ and Fig. 10, let us analyze the possibility to increase the lifetime of the stationary operation mode s of the oscillator at a given oscillation amplitude and frequency using the external control parameters. In accordance with Eq. ͑42͒ and Fig. 10, the lifetime of the stationary operation mode s exponentially depends on the amplitude of the control force F 0 and, therefore, can be increased significantly by increasing the amplitude of the control force F 0 . In the case where the fluctuation-dissipation theorem is valid and at the given oscillation amplitude and frequency, ␦ depends only on temperature ͓see Eq. ͑21͔͒. However, due to the Q-factor being almost inversely proportional to temperature ͑see Table I͒, ␦ can only slightly be changed with temperature. As a result, the lifetime of the stationary operation mode s also weakly depends on temperature.
The critical value c of frequency corresponding to the critical value E c of the oscillation energy is determined by the relation where = / ͑2͒ is the frequency of the control force. From Eqs. ͑41͒ and ͑42͒ it follows that ⌬ c / ϰ ͑F 0 / F 0c ͒ 1/2 for F 0 / F 0c ӷ 1. Such a dependence of the critical frequency on the control force amplitude correlates with the squarerootlike dependence for the lower limit of the ratio / 0 on F 0 / F 0c presented in Figs. 6͑a͒ and 6͑b͒. This is because, in the both cases, we consider the maximum relative 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 dT s / dE ϳ T s / E. Therefore, the lifetime of the stationary operation mode should be finite for any strongly anharmonic oscillator.
Up to now atomistic simulation of nanotube-based oscillators was restricted to sizes of hundreds of nanometers and simulation times of tens of nanoseconds. [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]39 The proposed phenomenological model taking into account thermodynamic fluctuations makes it possible to extend the simulation time up to 1 ms. Let us discuss the possibility to use this model not only for long simulation times but also for sizes of the oscillator greater than 3 nm.
The proposed model based on motion equation ͑11͒ has no size restrictions and can be applied to any size of the oscillator in the case if: ͑1͒ there are no resonances between the telescopic oscillation and other vibrational modes of the system and ͑2͒ the phenomenological parameters for the ran- dom noise force and friction force are available. Let us analyze the restrictions imposed on the nanotube length by condition ͑1͒. Among the low-frequency nanotube vibrational modes, there are long-wave acoustic vibrational modes, 57-60 squash modes, 61 and relative vibrational modes of the walls. 62,63 According to our previous calculations, 24 for the ͑5,5͒@͑10,10͒ nanotube of 3-nm length, the fundamental frequency of the doubly degenerate transverse acoustic modes is 1.5 THz, the fundamental frequency of the squash mode is about 1 THz and the frequencies of the nonaxial translational relative vibrations of the inner wall inside the outer wall are 0.25-0.33 THz. For this nanotube, there is no resonance between the telescopic oscillation and any other nanotube vibrational modes. With increasing the nanotube length, the fundamental frequencies of the longitudinal and torsional acoustic modes decrease as ϰ1 / L ͑Refs. 57, 58, and 60͒ ͑where L is the nanotube length͒ while the frequencies of the relative vibrations of the nanotube walls and of the squash modes weakly depend on the nanotube length. Since the frequency of the telescopic oscillation is inversely proportional to the nanotube length ͓see Eq. ͑14͒ for the case s ϰ L͔, these modes do not become resonant with the telescopic oscillation. However, from the Euler-Bernoulli beam theory 64 it follows that the fundamental frequency of the transverse acoustic ͑flexural͒ vibrational modes 57 should decrease with increasing the nanotube length faster than the frequency of the telescopic oscillation where Y Ϸ 1 TPa is the nanotube Young's modulus, 65-67 Ϸ 2.2ϫ 10 3 kg/ m 3 is the density of carbon atoms, I Ϸ R 4 is the nanotube areal moment of inertia, A Ϸ R 2 is the nanotube cross-sectional area and R is the outer wall radius. From the equality of the fundamental frequency of the transverse acoustic modes ͓see Eq. ͑44͔͒ and the frequency of the telescopic oscillation ͓see Eq. ͑14͒ for the case s = L / 3͔, 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 L max Ϸ 140 nm. For longer nanotubes, coupling of the telescopic oscillation with the transverse acoustic modes should be taken into account.
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, 24 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 fluctuationdissipation 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 fluctuation-dissipation 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 in ␦ 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 Q-factor strongly decreases͒ with an increase in the oscillation amplitude due to the excitation of low-frequency vibrational modes. 26,35 Thus, for an oscillator with a certain length the minimum level of fluctuations ␦ 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 24 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 L Ӷ L max , see Eq. ͑45͔͒. In the same case, the oscillation energy is proportional to the oscillator length. Hence, one gets the minimum value of ␦ 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. For example, based on the data of Table I for the 3.1 nm ͑5,5͒@͑10,10͒ nanotube-based oscillator, one can estimate that, for the oscillator longer than 35 nm, ␦ should be below 0.4. For F 0 / F 0c = 2, the lifetime of the stationary operation mode would already exceed 1 s. Thus, by the example of the gigahertz oscillator, we demonstrated that thermodynamic fluctuations can impose crucial restrictions on sizes of NEMS for which control of the NEMS operation is possible.

V. CONCLUSION
In summary, we performed 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 ͑amplitude F 0 , angular frequency , and initial phase shift between the control force and the velocity of the movable wall͒ at which the stationary operation mode is possible were determined. In particular, it was shown that the ranges of and which correspond to the stationary operation mode become wider with increasing the amplitude of the control force.
Significant thermodynamic fluctuations in the nanotubebased gigahertz oscillator were observed by the MD simulations of damping oscillations. The multiscale simulations ͑both MD simulations and simulations with the phenomenological model͒ revealed that the fluctuations cause the breakdown of the stationary operation mode. The average lifetime of the stationary operation mode was found to decrease with increasing the level of fluctuations or decreasing the amplitude of the control force. These dependences were explained through the analysis of the oscillation energy distribution function calculated on the basis of the Fokker-Planck equation. The investigations performed showed that the thermodynamic fluctuations impose restrictions on the sizes of NEMS for which the control of the NEMS operation is feasible.