Computational Analysis of Strain ‐ Induced Effects on the Dynamic Properties of C 60 in Fullerite

: A hybrid discrete ‐ continuous physical and mathematical model is used to study what deformation characteristics cause the rolling effect of C 60 fullerene in a fullerite crystal. The interac ‐ tion of fullerene atoms with surrounding molecules is described using a centrally symmetric inter ‐ action potential, in which the surrounding molecules are considered as a spherical surface of uni ‐ formly distributed carbon atoms. The rotational motion of fullerene is described by the Euler dy ‐ namic equations. The results of a numerical study of the influence of the rate, magnitude, and di ‐ rection of strain on the dynamic characteristics of the rotational and translational motion of C 60 full ‐ erene in a crystalline fragment are presented.


Introduction
Interest in fullerene-containing materials is due to a wide variety of properties that are widely used in biology, medicine, chemistry, electronics, and materials science [1][2][3][4][5][6][7]. Fullerenes are a family of molecules containing from 20 to 980 carbon atoms [8,9]. Fullerenes represent a dense closed structure in the form of a convex polyhedron. Currently, the most common fullerene is the C60 molecule, which is an icosahedral hollow structure consisting of 60 sp2 hybridized carbon atoms lying on the surface of a sphere [9][10][11]. The C60 molecules have high mechanical rigidity, stability, and strength but a fullerite crystal consisting of the C60 molecules is a fairly soft material [12][13][14][15]. Fullerite in the face-centered cubic (fcc) phase has the lowest density (1.6 g/cm 3 ) compared to graphite and diamond. Due to intermolecular forces caused by van der Waals forces, fullerenes at nodes of the fullerite crystal lattice at room temperature perform the rotational motions in the gigahertz frequency range [16][17][18][19][20]. The rotation frequency of fullerenes can be increased using an external magnetic field [21] and high-power pulsed laser radiation [22]. The addition of even a small amount of fullerenes can significantly change the properties of existing materials [23]. Fullerite is also subjected to high temperatures, laser radiation, high pressure, and shear deformation to produce materials with new properties [24][25][26]. The resulting materials have a wide range of physicochemical and mechanical properties [27,28]. However, the synthesized samples are mechanical mixtures of different phases of the C60 polymers. Fullerenes are also used in tribological compositions to reduce aging and wear of parts in frictional interaction [29][30][31][32]. In this regard, it is of great interest to study the properties of fullerenes under conditions of force action on the material.
In this article, we study the effect of the simplest deformation of a fullerite crystalline fragment on the dynamic characteristics of the C60 molecule using the methods of classical mechanics [33][34][35][36][37]. A computational analysis of the behavior of the fullerene molecule depending on the speed, direction, and magnitude of the indentation deformation has been carried out.

Physical Statement of the Problem
A fullerite fragment consisting of 63 fullerene molecules is considered (see Figure 1). At room temperature, the C60 molecules form a face-centered cubic lattice with a lattice parameter of d = 1.417 nm [38]. A fragment of the fcc crystal structure is a cube, the opposite faces of which are at a distance of 2d. The central fullerene molecule is composed of 60 carbon atoms (C60, 720 am). The dynamic state of the central molecule is determined by the nature of the interaction with all surrounding fullerenes. Therefore, it plays a special role in the presented model. The remaining 62 molecules are modeled as spherical particles. At the initial moment of time, the distance between the centers of the two nearest molecules is around 1.002 nm [38,39]. The interaction of the central molecule with smoothed C60 is described using interatomic interaction potentials based on the Lennard-Jones potential [40,41]. Thus, in this article, the interaction of any carbon atom of the central fullerene with any smoothed fullerene depends only on their mutual position but does not depend on the position of any other particles, as in many-particle potentials [42][43][44]. At a certain point in time, a fullerite fragment is deformed by moving a group of molecules towards the central fullerene at a given speed. For simplicity, we do not take into account the influence of intramolecular vibrational modes, i.e., the central molecule is a set of material points (atoms) rigidly interconnected. With this assumption in mind, the central molecule has three degrees of freedom of translational motion and three degrees of freedom of rotational motion.

Mathematical Statement of the Problem
Let us place the global Cartesian coordinate system (xyz) in the initial symmetrical position of the center of mass of the central fullerene. The local Cartesian coordinate system (ξηζ) is fixed on the rotating central fullerene and its oscillating center of mass.
The interaction potential between a carbon atom of the central C60 and a smoothed spherical fullerene has the following form [33,34,40,41] (1) Here S is the area per carbon atom; a is the radius of the fullerene; ε, σ are the Lennard-Jones parameters [33]; rik is the distance from the center of the kth carbon atom of the central C60 molecule to the center of the ith smoothed fullerene, , n = 60 is the number of carbon atoms in the C60 molecule; N = 62 is the number of smoothed fullerenes.
Potential (4), based on the Lennard-Jones potential, assumes a uniform distribution of potential energy over the surface area of the sphere. We use potential (4) because it allows us to simplify the analysis of the behavior of fullerenes due to the exclusion of highfrequency oscillations from the solution.
The Euler dynamic equations are used to describe the rotational motion of the central fullerene molecule around its center of mass [33,34]: where aij ( ) are the components of the matrix of direction cosines, which connects the systems of Cartesian coordinates xyz and ξηζ.
The projections of the moment of external forces on the x-, y-, z-axes are determined from the following formulas: Here Mik is the moment of external force Fik, taken relative to the center of mass of the molecule; Xik, Yik, Ζik are the projections of force Fik on the x-, y-, z-axes; rkc = (xkc, ykc, zkc) is the radius vector of the atom relative to the center of mass of the molecule; i, j, k are the unit vectors of the fixed body coordinate system.
The force effect from the side of the ith smoothed fullerene on the kth atom of the rotating C60 molecule is determined by the following formulas: The angular velocity components can be written in terms of the Euler angles φ, ψ, θ in the following form: where the dot above the function name corresponds to the time derivative. The dynamic state of the central fullerene is determined by the nature of the interaction of sixty carbon atoms with all smoothed fullerenes. Consequently, the movement of the center of mass of the central fullerene obeys the following law: where Mf is the mass of the fullerene; VC is the velocity of the center of mass of the central fullerene; C r is the position vector of the point of the center of mass of the central fullerene. The initial conditions for solving the system of differential Equations (2-4), (10-13) are given in the following form t = 0: , , , , where the subscript i refers the ith smoothed fullerene, the superscript 0 refers to initial value. Thus, we use a hybrid discrete-continuous mathematical model in which the central fullerene has icosahedral symmetry, and the surrounding fullerenes are described as hollow rigid spheres. The system of first-order evolution Equations (2-4), (10-13) with initial conditions (14) and (15) is solved using the Runge-Kutta fourth-order method [45,46], which makes it possible to determine the characteristics of the translational and rotational motion of the central fullerene surrounded by smoothed fullerenes.

Results and Discussion
We will carry out the following computational experiment to study the influence of deformation on the behavior of fullerenes in a molecular crystal. We consider the central fullerene C60, which is at rest at the initial moment of time. In reality, the deformation force will most likely have non-zero magnitudes of projections along the ξ-, η-, ζ-axes. However, for the convenience of analyzing the force action, we introduce some simplifications to identify causal effects. We place this C60 fullerene in such a way that the vector of the angular velocity of rotation caused by the action of the surrounding fullerenes in the absence of deformation is directed along the ξ-and x-axes (see Figure 2), which coincide at the initial moment of time. In this position, the center of mass of the central fullerene in the absence of deformation (Δx = Δy = 0) must retain its original position. Figure 2 shows that ωξ(t) is a periodic and sign-changing function. This suggests that the fullerene rotates in one direction and then in the opposite direction. The angular velocity components ωη and ωζ are close to zero. The average absolute value of the angular velocity is 408 rad/ns. This result qualitatively agrees with the experimental data according to which the frequency of rotation of the C60 molecule in fullerite at room temperature is about 10 11 Hz [16][17][18][19][20].

Strain Magnitude and Direction Effect on the Central C60
We consider the deformation of a fragment of a molecular crystal along the positive x-axis (the red arrow in Figure 1). Let the extreme row (x = −d), consisting of 13 smoothed molecules, move towards the central fullerene with the absolute strain rate vdef = v*, where v* = 0.5d/(tz − t0) is specific speed, d = 1.417 nm is lattice parameter, t0 = 0, tz = 10 ps are start and end time of simulation. To analyze the external force acting on the central fullerene, we determine the position to which the central fullerene passes by finding the average value of the coordinate functions of time (xavg, yavg, zavg) and the maximum amplitude of oscillations around this position Amax. The calculations of these parameters are carried out from the moment of time at which the deformation stops to the final time of the simulation. Figure 3 shows the curves corresponding to the displacement of these molecules by a distance Δx in the range from 0 to 0.25d. Changes in other coordinates of the center of mass are not given because in this case, there is a rectilinear motion of the central molecule. As can be seen from Figure 3, the stronger the indentation deformation of the crystalline fragment along the x-axis (x-deformation), the stronger the average deviation of the central molecule from its initial position and the greater the amplitude of oscillations relative to the average position.
Thus, the central fullerene settles into average positions xavg = −0.15, −0.85, −2.34 pm with a maximum amplitudes Amax = 0.12, 0.31, 0.55 pm at Δx = 0.05d, 0.15d, 0.25d, respectively. A negative value of x(t) indicates that the molecule is moving against vdef. This is due to the fact that the movement of the most distant smoothed molecules leads to the appearance of an attractive force-directed, opposite to the direction of deformation. If the deformations involve two adjacent rows, including 25 molecules, then the nature of the movement of the central molecule will change somewhat (see Figure 4). The central molecule will move in the direction of vdef. In this case, an increase in the magnitude of the deformation leads to a decrease in the amplitude, an increase in the displacement of the molecule, and an increase in the oscillation frequency. The average positions of the molecule xavg = 35.7, 71.2, 106.3 pm, and the amplitude Amax = 20.66, 9.89, 8.76 pm for Δx = 0.05d, 0.1d, 0.15d, respectively. It can be seen from Figure 4 that the positions of the central C60 is 1-2 orders of magnitude larger than in the case of 13 molecules. With all the differences in the results presented in Figures 3 and 4, the movement of the central molecule retains its straightness in both cases. A number of differences are due to the direction and magnitude of the force that arises during deformation along the x-axis and acts on the molecule under consideration. Therefore, for simplicity of analysis, we will further consider the case with 13 molecules, which (taking into account some features and up to certain values of the force impact) qualitatively reflects the effect of deformation on the dynamics of the movement of the central fullerene. The problem considered above, when the direction of deformation is collinear to the vector of the angular velocity of rotation of the central fullerene, occurs extremely rarely in practice. Therefore, we also consider the problem when vdef and ω are non-collinear vectors.
Let the deformation be directed along the y-axis (Δy > 0, Δx = Δz = 0) (y-deformation), i.e., vdef (purple arrow in Figure 1) perpendicular to the angular velocity of rotation. Then the displacement of the central fullerene is described in Figure 5. It can be seen that the change in the y-coordinate in Figure 5 is similar to the change in the x-coordinate in Figure  3. However, as can be seen from Figure 5, the y-deformation causes a smaller displacement of the central fullerene than the x-deformation since the y-deformation leads to a qualitative change in the dynamics of the movement of the central fullerene. The central fullerene settles into yavg = −0.12, −0.69, −1.94 pm with Amax = 0.07, 0.11, 0.2 pm at Δy = 0.05d, 0.15d, 0.25d, respectively. The average positions along the deformation axis decreased by 16.9-18.6%, and the maximum amplitude Amax decreased by 40.8-67.2% compared to the xstrain (Figure 3). This result is explained by the fact that the motion of the central fullerene during ydeformation is not rectilinear. The y-deformation causes the central fullerene to move in the yz-plane, i.e., its movement is perpendicular to the vector of the angular velocity of rotation (see Figure 6). In this case, the amount of displacement along the z-axis increases with an increase in the amount of y-deformation. The central fullerene settles into a positions zavg = −0.09, −0.47, −1.28 pm with Amax = 0.75, 1.27, 2.55 pm at Δy = 0.05d, 0.15d, 0.25d, respectively. Note that the y-deformation, which affects 25 smoothed molecules, also forces the central molecule to move in the yz-plane. However, in this case, the indentation deformation can strongly affect the angular velocity of the central fullerene thus the analysis of the numerical solution becomes very difficult.
The effects of rolling and force reduction are shown schematically in Figure 7. Figure  7a describes the rectilinear motion that occurs when the angular velocity ω and the force F lie on the same straight line passing through the center of mass of the C60 molecule or when ω = 0. Figure 7b schematically shows the motion that occurs when the angular velocity ω and the force F are perpendicular vectors. We can see the effects of the angular velocity direction because we used the hybrid discrete-continuous mathematical model. The calculation results show that the models of smoothed molecules [40,41,47,48] applied to all fullerenes, including the central fullerene, will give a qualitatively incorrect result when the axis of rotation does not coincide with the direction of deformation (see Figure  7). Because the rotation of fullerenes is not present in these models (ω = 0). Thus, there will be no movement of fullerenes along the z-axis under the conditions of the computational experiment presented in this article (see Figure 7). At the same time, these models are useful for quick calculations and for isolating the various physical effects from highfrequency oscillations.  It should be noted that the described movement of the central fullerene in the z-direction (see Figure 7b) is due to the gyroscopic effect [33,37,[49][50][51], which occurs when the projection of an external force on the plane of motion of carbon atoms is non-zero. Since this feature of fullerene motion is determined by its geometry, it can be used in other molecular environments to reduce the force effect on the material. We assume that these calculation results can partly provide an additional explanation for the sliding or rolling effect [29], according to which the friction coefficient decreases with an increase in external loads or sliding speed. In this regard, we will also conduct a numerical study to test the influence of the strain rate on the gyrodynamics of the central fullerene to try to explain why the sliding speed can reduce the coefficient of friction (the force reduction effect).

Strain Rate Effect on Central C60
Let the outer row of molecules move in the positive x-direction for a distance Δx = 0.25d with a velocity vdef that varies from 1v* to 5v*. As can be seen from Figure 8, the greater the strain rate, the greater the amplitude of the oscillation. The calculation of the mean value of the x(t) function shows that at Δx = 0.25d the center of the central fullerene shifts to the position xavg = 2.33 ± 0.01 pm for any value of vdef. However, the maximum amplitudes Amax are 0.54, 0.7, 1.28 pm for strain rates vdef = 1v*, 3v*, 5v*, respectively. This suggests that higher values of the strain rate increase the effect on the central fullerene and, consequently, on the fullerite fragment.  Figure 9 shows that at Δy = 0.25d the fullerene assumes the position yavg = 1.937 ± 0.006 nm. As in the previous case (Figure 8), the strain rate has almost no effect on the average position that the fullerene assumes due to deformation. However, Figure 8 shows that xavg is 20.1% larger than yavg (Figure 9). This is due to the fact that there is a non-zero displacement of the central fullerene in the direction of the z-axis, characterized by the average positions zavg = −0.127, −0.134, −0.142 pm and the maximum amplitudes Amax = 0.255, 0.327, 0.423 pm for strain rates vdef = 1v*, 3v*, 5v*, respectively (see Figure 10). Therefore, the greater the strain rate, the greater the displacement of the fullerene in the direction perpendicular to the angular velocity and strain direction. In other words, the rolling effect [29] increases with increasing strain rate. This result suggests that the fullerene gyrodynamics are also sensitive to changes in the strain rate.

Conclusions
In this article, classical molecular dynamics calculations have been carried out to study the influence of nanoindentation on the dynamics of a rotating C60 fullerene in a crystalline fullerite using a hybrid discrete-continuous mathematical model. The results of the calculations showed that the behavior of the fullerene depends on the mutual direction of indentation deformation and the angular velocity. Their non-collinearity leads to displacement of the fullerene in the direction perpendicular to the direction of defor-mation (rolling effect) and a significant decrease in the directed force effect on the C60 fullerene. This displacement is carried out in the plane of motion of the carbon atoms of the fullerene and increases with an increase in the rate and magnitude of deformation.