Md-Simulation of Fullerene Rotations in Molecular Crystal Fullerite

Abstract

The present paper describes rotations of C₆₀ fullerene molecules in the solid phase of a fullerite. The conducted studies show that these relatively large molecules rotate according to the same laws as macroscopic bodies, i.e., according to the laws of classical mechanics. The performed calculations confirm that fullerene rotations do not cause friction. We suggest a method for a strong increase in the internal energy of the material that does not lead to its destruction. It is theoretically shown that in standard fullerite, in the absence of electric and magnetic fields, fullerene rotations occur with an average angular frequency of 0.34·× 10¹² rad·s⁻¹, which is consistent with the experimental data obtained using nuclear magnetic resonance. By means of calculations, we found that alternating magnetic fields of a certain configuration wind fullerenes encapsulated by iron. In this case, two temperatures arise in the fullerite crystal: a high rotational temperature and a vibrational temperature close to normal. For the purpose of determining this velocity, as well as the nature of rotations, the present paper suggests a way of integrating the dynamic Euler equations for the projections of a molecule’s angular velocity vector onto the coordinate axes associated with the fullerene. The stages of computer simulation of fullerene movements, which was carried out without using previously developed packages of molecular-dynamic modelling, are consistently described.

1. Introduction

A fullerite-molecular crystal with a high degree of crystalline order. The first researchers to observe a solid fullerite were V. Kretchmer and D. Hoffman (May 1990, the Institute of Nuclear Physics, Heidelberg, Germany). It was discovered that at room temperature a fullerite is observed in the so-called plastic phase which is characterized by rotations of fullerenes in the nodes of the face-centered cubic lattice (the FCC lattice). At T = 260 K, a structural transition occurs, and the FCC lattice is transformed into a simple cubic lattice (the SC lattice). This state of a fullerite is called the low-temperature phase. It is believed that in this phase fullerenes do not rotate. However, a number of recent studies claim that rotations of a specific type take place and they are implemented with a lower frequency, compared with the plastic phase of a fullerite. The main purpose of the work is to develop a simple and effective model of fullerene interactions in a fullerite material. We suggest equations for rotations and vibrations and give a scheme for their numerical implementation. The calculations found that fullerenes in the plastic fullerite rotate as well as vibrate. The frequency and amplitude of vibrations, as well as the average frequency and nature of the rotations are determined. It is shown that within each group of motions the kinetic energy is distributed uniformly in terms of the degrees of freedom. At the same time, there is a significant redistribution of energy among groups of motions. The specific energies in the groups depend on the temperature of the system. Under normal conditions, it is possible to observe a significant redistribution to the advantage of rotations. Let us consider some examples of rotations associated with fullerenes. Refs. [1,2] investigated translational and rotational movements of CH₄ inside a fullerene with an open cell. Ref. [3] demonstrated a single-molecule switch based on tunnel electron rotations of the tetrahedral cluster Sc₃N in an icosahedral cell of a C₈₀ fullerene. Ref. [4] considered rotational dynamics of C₆₀ using nuclear magnetic resonance of carbon isotopes. Ref. [5] analysed the same problem by measuring Raman active intermolecular modes at different temperatures Ref. [6] studied low-energy electronic states of a rotating fullerene in the framework of the continual model on the basis of experimental data on rotating C₆₀ molecules in a fullerite. Ref. [7] showed that dimers of a C₇₀ fullerene freely rotate around the short axis of the molecule. Ref. [8] studied rotations of fullerene chains in single-walled nanotube structures using low-voltage translucent high-resolution electron microscopy. Ref. [9] analysed the spin-lattice relaxation times in a polycrystalline C₆₀ fullerene. Ref. [10] reviewed pseudo-rotations of fullerene ions. Ref. [11] investigated stability of the C₂₀ @ C₈₀ nanoparticle and rotations of its inner shell. Ref. [12] identified that laser pulses cause rotations of elliptical polarization in solutions of C₇₀ fullerenes. The study conducted by Ref. [13] focused on orientational properties of a C₆₀ monomer. Dynamics of adatoms in exohedral C₆₀ fullerene complexes of rare atoms were studied by Refs. [9,14] investigated the effect of the functional group sizes on mobility of electrons in cells of fullerene derivatives. Ref. [15] calculated the infrared spectrum of the endohedral metal fullerene Li⁺ @ C₆₀. R.M. Ref. [16] carried out calculations for the main spectrum of vibrational-rotational motions in FCC (fullerene) lattices. Ref. [17] determined the thermal expansion coefficients of the pure fullerite on the basis of X-ray powder research. Ref. [18] described the theory of translational-rotational relations in terms of the microscopic Hamiltonian and the resulting Landau free energy. Refs. [19,20] investigated photo-absorption spectra of icosahedral fullerenes. Ref. [21] studied the local magnetic properties of carbon nanotubes using nuclear magnetic resonance. The authors of the works given in Refs. [22,23,24,25,26,27] perform molecular dynamics calculations using the Tersoff potential or its modifications (Brenner, REBO, AIREBO). These potentials are excellent for modelling carbon structures and are most commonly applied in calculations of this kind which are usually carried out in software packages, such as LAMMPS [27]. In Refs. [28,29,30] the authors suggest alternative methods for calculating fullerene vibrations. In Ref. [28] there is a consideration of modelling the vibrational modes of C₆₀ and C_70. For the purpose of analyzing the state of the considered molecular structures the authors use the finite element method. In Ref. [29] the method of atomistic structural mechanics is used to study vibrational characteristics of various fullerenes. In Ref. [30] the authors investigate axisymmetric vibrations in the case of spherical fullerenes. The nonlocal theory of elasticity is used. From the presented brief review, it follows that the work on the study of fullerene rotations is mainly experimental. There are already good theoretical models for determining the dynamics of carbon atoms, which determines the state of an individual C₆₀. However, there is no closed description of the dynamics of fullerenes in fullerite material.

2. Assessment of the Applicability of the Model of Classical Mechanics

It should be noted that, in assessing the impact of the quantum nature of the problem, it is necessary to determine the de Broglie wavelength of a fullerene in the fullerite which depends on the mass of the C₆₀ molecule and its thermal velocity. The average velocity of fullerene thermal motions is determined on the basis of calculated data on the fullerene rotation frequency, as well as the available experimental data. Introducing the found velocity into the ratio for the de Broglie wavelength, we obtain that the wavelength of the fullerene molecule is substantially smaller than its diameter. Thus, for the purpose of describing the phenomenon of these particle rotations, an approach applied in the framework of classical mechanics models is justified. Let us evaluate the applicability of classical Newtonian mechanics as opposed to quantum mechanics for the case under consideration. The de Broglie wavelength is defined as:

$$\lambda =\frac{h}{M\upsilon }$$

(1)

where h = 6.626 × 10⁻³⁴ J × s is the Planck’s constant; M = 1.197 × 10⁻²⁴ kg is the mass of the C₆₀ fullerene molecule; $\upsilon$ is its average linear rotational velocity which can be estimated based on experimental data on fullerene rotation frequency. Taking into account the linear relation of the angular rotation frequency and the velocity, $\upsilon$ = ωa, where a = 0.357 × 10⁻⁹ m is the radius of the fullerene sphere; also ω = 2π$\nu$, where $\nu ={\tau }^{-1}={\left(9.1\times {10}^{-12}\right)}^{-1}=1.1\times {10}^{11}{\mathrm{s}}^{-1}$ is the circular frequency of rotation [4]. Under these conditions, the average rotational velocity of the C₆₀ fullerene molecule is about 246.7 m/s, and its de Broilean wavelength is approximately 0.224 × 10⁻¹¹ m. If to assume that the characteristic velocity included in (1) is equal to the average vibration velocity of the fullerene which, according to our data, is 10² m·s⁻¹, we obtain $\lambda =0.553\times {10}^{-11}\mathrm{m}$. At the same time, the diameter of the C₆₀ molecule is 0.71 × 10⁻⁹ m. Thus, the size of the fullerene is two orders of magnitude larger than the de Broglie wave for this particle. Along with this, the obtained wavelength is substantially less than the spacing of the fullerite crystal lattice, which is equal to d = 1.002 × 10⁻⁹ m and even noticeably less than the distances between the carbon atoms in the fullerene (0.142 × 10⁻⁹ m). In this regard, the fullerene which is in dynamic equilibrium in the material of the plastic fullerite must behave like a particle, not like a wave. Therefore, in describing rotations of such particles, the approach implemented in the framework of classical mechanics is fully justified. The calculations of the state of plastic fullerite carried out in the present work using the model of classical mechanics made it possible to find the rotation frequency of fullerenes, as well as the frequency and amplitude of the vibrations.

2.1. Packings of Fullerenes in a Fullerite

C₆₀ molecules are the most common closed form of a surface carbon crystal. They differ from other fullerenes in the perfection of the form due to the fact that all atoms in these molecules are located on the spherical shell. The nature of the location corresponds to icosahedral symmetry.

Left part of picture 1 shows that a hexa- and pentagonal packing of carbon atoms is observed on a spherical surface. The straight-line segments conventionally show C-C bonds, including the double bonds which are present in the structure of the molecule. Considering the model of a smoothed fullerene developed [31], it can be noted that there is a spherically symmetric potential, and the radial component of that has a minimum which is located at a distance of 2.88 Å from the surface (Figure 1, right). This value can be considered the average distance of the potential interaction energy minimum in a real system of two fullerenes. In the first approximation, a fullerite can be considered as a set of contiguous fullerenes located in potential wells of individual molecules. The presented distribution can be obtained on the basis of the potential of atom–atom interaction in two ways. Add the energy of all cross interactions of carbon atoms belonging to two different fullerenes and average the result over the rotation angles of each fullerene, or introduce into consideration the energy of the interaction of an elementary site of a homogeneous interacting sphere and integrate over the area of each of the spheres. As a result, both methods will lead to the same result, shown in Figure 1. At the same time, taking into consideration that C₆₀ molecules are rather large and experience rotations, the material under study can be reasonably considered as a molecular mechanism with numerous degrees of freedom.

It is possible to distinguish SC nodes and FC nodes in the face-centered lattice of a plastic fullerene. Let us take an FC node as a central node of a certain minimum number of C₆₀ molecules still affecting the dynamic state of the selected central fullerene. Such a minimal fragment (see Figure 2) can always be distinguished since the Van der Waals interactions within the structure depend significantly on the distances. Moreover, these forces are such that they quickly decrease when increasing the distance between atoms.

The FCC-lattice node selected for consideration is surrounded by twelve nodes located at the distance $\frac{d\sqrt{2}}{2}$, here d is parameter of lattice. They will be called nodes of the first environment echelon. Six more nodes will be found at the distance d (nodes of the second environment echelon). Eight more nodes will be found at the distance $\frac{d\sqrt{6}}{2}$. The minimal fragment is a system of nodes including the central node under consideration and the nodes of the three environment echelons. Despite the fact that, compared with repulsion forces, the Van der Waals forces are considered long-range, other fullerenes outside the selected fragment will no longer influence the state of the central fullerene since they are quite far from it.

If an SC-node is selected as the central node of the minimum fragment, the fragment will also contain twelve nodes of the first environment echelon, six nodes of the second echelon and eight nodes of the third echelon, with the same distances as in the case of the central FC-node. Therefore, rotation of the central fullerene will be equivalent to the originally designated case.

It is most convenient to consider the crystal lattice of a fullerite as a system of material points (carbon atoms) with superimposed covalent bonds acting on the surfaces of each individual fullerene. Within the framework of such a presentation, an explanation of the reasons for fullerene rotations can be as follows. A mechanical system consisting of a large aggregate of carbon atoms assembled into fullerene cells tends to minimize cross-interactions, which results in rotations of C₆₀ molecules. For ideal systems this tendency is manifested in the nature of ongoing oscillations.

In order to calculate dynamics of a system with multiple degrees of freedom, it is necessary to try to limit the number of these degrees, while maintaining basic qualities of the system. We considered each internal fullerene of the selected fragment as a free solid body which, as known, has six degrees of freedom. Since for the eight fullerenes constituting the third echelon of the environment we use dynamic symmetry, nineteen fullerenes of the fragment remain free indeed. Therefore, the representative element of the material under consideration has 114 degrees of freedom.

In the presented mathematical model, all the fullerenes of the minimal fragment are supposed to ensure the correct effect on the central fullerene with their force fields. These force fields generate rotations and vibrations of C₆₀. Naturally, with an increase in the distance from the central fullerene, the contributions of impacts which determine its dynamic state decrease. The third environment echelon fullerenes are the last among those that still have a noticeable effect on the state of the central fullerene. However, these edge fullerenes themselves are in an incomplete environment. Therefore, due to the unilateral effects of fullerenes adjacent to them, they can acquire a different kind of rotation. To prevent this, we assume that eight fullerenes of the third environment echelon rotate in the same way as the central fullerene; moreover, they have the same vibrations of their mass centers as the central fullerene. This is exactly the geometric and dynamic symmetry of the crystal fragment.

In the present work the deterministic description proposed by Euler is used to determine movements of fullerene cells. It consists in calculating the moments of forces rotating a fullerene, as well as integrating the equations for angular velocities, rotation angles and centers of masses of fullerene molecules.

2.2. The Equation of Rotational Motion for a C₆₀ Molecule

All the turns which fullerenes make around their own centers of mass are considered as rotations. Vibrations are translational movements of fullerenes along with their centers of mass. During translational motions, all elements of motion of a solid body (a fullerene molecule in the considered case) are determined by the motion of one of its points (which in the given case is the center of mass of a C₆₀). Thus, if we have the equations which determine motions of the centers of mass, as well as the necessary initial conditions for these points, the vibrations can be completely determined. In connection with the foregoing, in the plastic phase of fullerite, fullerenes themselves both rotate and vibrate. They rotate around their own centers of mass, and vibrational oscillations occur near their average positions. Let us now write the equations of rotational motion of a fullerene, which in classical mechanics are called the dynamic Euler equations.

Let be, Oxyz will be the fixed frame of reference, and $C\xi \eta \zeta$ will be associated with the rotating C₆₀ molecule, as shown in Figure 3.

Then the values of force effects on a single atom of the C₆₀ molecule can be determined as follows:

$$X_k=-{\displaystyle \sum _{j=1}^{N_p}\frac{\partial U}{\partial x}\left(r_{kj}\right)};$$

(2)

$$Y_k=-{\displaystyle \sum _{j=1}^{N_p}\frac{\partial U}{\partial y}\left(r_{kj}\right)};$$

(3)

$${{\rm Z}}_k=-{\displaystyle \sum _{j=1}^{N_p}\frac{\partial U}{\partial z}\left(r_{kj}\right)}.$$

(4)

Here, X_k, Y_k, Ζ_k are the projections of the all the resulting forces (acting on the k-th atom of the molecule under consideration from other atoms of the crystal) on the axis of the stationary reference system; $r_{kj}=\left|r_k-r_j\right|$, r_k is the radius vector of the k-th atom of the central molecule C₆₀; r_j is the radius vector of the j-th atom of the surrounding molecules; N_p is the number of carbon atoms in the minimum fragment of the face-centered structure, with the exception of central fullerene atoms, i.e., Np = 1560. U(r_kj) is the potential of the “atom–atom” interaction, which can be found in the book by [32]. We consider atomic vibrations in a simplified form as harmonic vibrations of a certain frequency and amplitude. This approach can be justified since in the plastic phase of fullerite C₆₀ molecules rotate freely, chemical bonds are not created and are not destroyed in the process of dynamics, and atomic vibrations have practically no effect on fullerene rotations, moreover, fullerenes are characterised by a very low friction [33,34,35]. Naturally, models in which the bonded potential is taken into account [25,26,27,28,29,30] are more rigorous; however, in connection with the foregoing, consideration of atomic vibrations with the use of bonded potentials does not provide a significant correction. Moreover, the suggested model is significantly simpler and more economical. Nevertheless, even this model, with good optimization of the calculation program, requires up to six hours of estimated time on a stationary PC with an intel core 7/9 processor of the 6th generation or higher, or AMD ryzen 5/7 with a RAM frequency of at least 3400 MHz with a step of 10⁻⁷ nm and a calculation time of 0.3 ns. It should be noted that the interaction potential of a specific bound carbon atom in a fullerene, in terms of another atom in a neighboring fullerene, is not centrally symmetric. However, in the case under consideration, three planes of symmetry can be distinguished in the force field of each particular atom; therefore, the real interaction can be replaced by some effective but centrally symmetric potential. Moreover, the problem of fullerene rotations is perfectly suitable for finding the effective potential form of a bound carbon atom.

The moment of forces external to the extracted fullerene is the sum of the moment of forces resulting from cross-interactions of the corresponding atoms M and the moment of magnetic forces M^(m). Let us consider the first of these moments:

$$M_{}^{}={\displaystyle \sum _{k=1}^S{M}_k^{},}{M}_k^{}=\left[r_{kC},F_k\right]=\left|\begin{array}{ccc}i& j& k\\ x_{kC}& y_{kC}& z_{kC}\\ X_k^{}& Y_k^{}& {{\rm Z}}_k^{}\end{array}\right|$$

(5)

Here, $M_k^{}$ is the moment of the resultant of all the forces acting on the k-th atom of the molecule taken with respect to the center of mass of the molecule; r_kC is the radius vector of the atom relative to the center of mass of the molecule; i,j,k are the units of the fixed coordinate system; S = 60 is the number of carbon atoms in the C₆₀ molecule.

The projections of the principal moment of forces acting from the structure on the rotating fullerene will be composed of the projections of the moments of forces acting on each atom of the molecule:

$$M_x^{}={\displaystyle \sum _{k=1}^S\left(y_{kC}{{\rm Z}}_k-z_{kC}{Y}_k\right)};$$

(6)

$$M_y^{}={\displaystyle \sum _{k=1}^S\left(z_{kC}{X}_k-x_{kC}{{\rm Z}}_k\right)};$$

(7)

$$M_z^{}={\displaystyle \sum _{k=1}^S\left(x_{kC}{Y}_k-y_{kC}{X}_k\right)}.$$

(8)

Here, $x_{kC}=x_k-x_C,y_{kC}=y_k-y_C,z_{kC}=z_k-z_C$. However, to use the dynamic Euler equations written in terms of the body-fixed axes, we need the projections of the moments of interaction between the C₆₀ molecule and the structure on the main central axes rigidly connected with the moving molecule (body-fixed axes).

As known from analytic geometry, the Cartesian coordinates of the vectors in both systems are related by the equations, presented by [36]:

$$M_{\xi }^{}=c_1{M}_x^{}+c_4{M}_y^{}+c_7{M}_z^{},$$

(9)

$$M_{\eta }^{}=c_2{M}_x^{}+c_5{M}_y^{}+c_8{M}_z^{},$$

(10)

$$M_{\zeta }^{}=c_3{M}_x^{}+c_6{M}_y^{}+c_9{M}_z^{}.$$

(11)

Here, $M_{\xi }^{},M_{\eta }^{},M_{\zeta }^{}$ are the projections of the moments of forces on the body-fixed coordinate axes. The projections of the moments of magnetic forces on the axis of the moving base are defined as follows:

$$M_{\xi }^{(m)}={\mu }_{\eta }{B}_{\zeta }-{\mu }_{\zeta }{B}_{\eta },$$

(12)

$$M_{\eta }^{(m)}={\mu }_{\zeta }{B}_{\xi }-{\mu }_{\xi }{B}_{\zeta },$$

(13)

$$M_{\zeta }^{(m)}={\mu }_{\xi }{B}_{\eta }-{\mu }_{\eta }{B}_{\xi }.$$

(14)

Here, ${\mu }_{\xi },{\mu }_{\eta },{\mu }_{\zeta }$ are the projections of the intrinsic magnetic moment of the fullerene on the axis of the moving coordinate system; $B_{\xi },B_{\eta },B_{\zeta }$ - are the projections of the magnetic induction vector which are found from the known values $B_x,B_y,B_z$ using direction cosines (19).

$$B_{\xi }=c_1{B}_x+c_4{B}_y+c_7{B}_z,$$

(15)

$$B_{\eta }=c_2{B}_x+c_5{B}_y+c_8{B}_z,$$

(16)

$$B_{\zeta }=c_3{B}_x+c_6{B}_y+c_9{B}_z.$$

(17)

External magnetic field is defined as:

$$B_x=B\cos ft,B_y=B\sin ft,B_z=0,$$

(18)

here $f$ is the frequency of generated external magnetic field.

In this case, the directional cosine c_i ($i=\overline{1,9}$) have the following form:

$$\begin{array}{c}{c}_1^{}=\cos \psi \cos \phi -\sin \psi \sin \phi \cos \theta ,c_2^{}=-\cos \psi \sin \phi -\sin \psi \cos \phi \cos \theta ,c_3^{}=\sin \psi \sin \theta ,\\ c_4^{}=\sin \psi \cos \phi +\cos \psi \sin \phi \cos \theta ,c_5^{}=-\sin \psi \sin \phi +\cos \psi \cos \phi \cos \theta ,c_6^{}=-\cos \psi \sin \theta ,\\ c_7^{}=\sin \phi \sin \theta ,c_8^{}=\cos \phi \sin \theta ,c_9^{}=\cos \theta .\end{array}$$

(19)

Then, the dynamic Euler equations can be written following [37]:

$$A\frac{dp}{dt}+\left(C-B\right)qr=M_{\xi }^{}+M_{\xi }^{(m)};$$

(20)

$$B\frac{dq}{dt}+\left(A-C\right)pr=M_{\eta }^{}+M_{\eta }^{(m)};$$

(21)

$$C\frac{dr}{dt}+\left(B-A\right)pq=M_{\zeta }^{}+M_{\zeta }^m.$$

(22)

Here, p,q,r are the angular velocity projections on the axis of the moving frame; A,B,C are the main moments of inertia of the molecule for its center of mass.

The latter equations are the equations of rotational motion of a molecule around its center of mass. These equations are closed by the kinematic relations connecting the projections of the angular velocity vector with the Euler angles and their derivatives:

$$p=\dot{\psi }\sin \theta \sin \phi +\dot{\theta }\cos \phi ;$$

(23)

$$q=\dot{\psi }\sin \theta \cos \phi -\dot{\theta }\sin \phi ;$$

(24)

$$r=\dot{\psi }\cos \theta +\dot{\phi }.$$

(25)

Solving these equations with respect to the time derivatives of the Euler angles, we obtain:

$$\dot{\psi }=\frac{p\sin \phi +q\cos \phi }{\sin \theta };$$

(26)

$$\dot{\theta }=p\cos \phi -q\sin \phi ;$$

(27)

$$\dot{\phi }=r-\left(p\sin \phi +q\cos \phi \right)\frac{\cos \theta }{\sin \theta }.$$

(28)

It should be noted that, in the calculation process, it is necessary to control the angle θ which can pass through the points 0, ±π, ±2π, etc. If such a transition is realized, accuracy of calculations can be lost. The presented system of differential Equations (20) to (28) requires specification of initial conditions. In the general case, they can be written as:

$$t=0:\psi ={\psi }_0,\theta ={\theta }_0,\phi ={\phi }_0,p=p_0,q=q_0,r=r_0$$

(29)

The problem presented in (20) to (29) is a Cauchy problem for six first-order differential equations. For each of the free fullerene of the minimal fragments it is formulated so that the final number of differential equations is one hundred and sixteen.

2.3. Motions of the Centers of Mass of Molecules

The motion equation of the centers of mass of polyatomic molecules is a necessary element in formulating the problem of fullerene rotations since the cross-forces of interatomic interactions included in the torques depend both on the orientation of fullerenes in space and their convergence or divergence, i.e., on oscillations of their centers relative to equilibrium positions:

$$M\frac{d{\mathit{\upsilon }}_{\mathit{C}}}{dt}=-{\displaystyle \sum _{j=1}^{N_p}{\displaystyle \sum _{k=1}^S\nabla U\left(r_{jk}\right)}}.$$

(30)

Here M is the mass of the molecule, v_s. is the velocity of the center of mass of the moving molecule; S is the number of atoms in the molecule, N_P – is the number of atoms that make up the environment in the selected fragment of the crystal; $r_{jk}=\sqrt{{\left(x_j-x_k\right)}^2+{\left(y_j-y_k\right)}^2+{\left(z_j-z_k\right)}^2}$, U(r) is the atom–atom interaction potential which depends on the type of atoms (in the considered case, carbon–carbon interactions), $\nabla$ is the gradient operator.

As can be seen from (30), motions of the center of mass of a polyatomic molecule are determined by the effect of all surrounding atoms on each atom of the molecule.

According to the Euler approach, the position of the body (in the present case, a C₆₀ molecule), along with the three coordinates of the center of mass, is also determined by the three Euler angles. Consequently, in the general case, along with Equation (30), it is necessary to use three more scalar equations for the angular velocity vector projections onto the moving molecules of the coordinate axis associated with the framework structure—the Euler dynamic Equations (20) to (22).

However, the equation of the center of mass motion, within the framework of each time step, can be integrated, to a certain extent, independently of the equations of angular velocities. In any case, it must be supplemented by the kinematic vector relation:

$$\frac{d{r}_C}{dt}={\mathit{\upsilon }}_{\mathit{C}}$$

(31)

Then the system of scalar equations defined by vector Equation (30) is closed, and it can be integrated numerically by stepwise methods of high accuracy.

The system of vector Equations (30), (31) must be integrated with the following elementary initial conditions:

$$t=0;{\mathit{\upsilon }}_{\mathit{C}}=0;{\mathit{r}}_c=0.$$

(32)

Here ${\mathit{r}}_C^{}$ is the radius vector of the center of mass of the molecule; ${\mathit{\upsilon }}_{\mathit{C}}^{}$ is the velocity of the center of mass.

Here we also deal with a Cauchy problem defined by the equations and relations given in (30) to (32). There are twenty-one problems of this kind in terms of the number of free fullerenes. In the process of solving these problems, one hundred and fourteen first-order differential equations are integrated.

2.4. Determining the Moments of Inertia of a Fullerene

The moment of inertia of a fullerene relative to some axis passing through its center of mass can be found on the basis of the following formula:

$$I=m_c{\displaystyle \sum _{i=1}^{60}{r}_i{}^2,}$$

(33)

where m_c is the mass of the carbon atom; r_i is the distance from this atom to the selected central axis. It is clear that for each specific axis direction this formula will result in a certain value, which will be different from the previous instantaneous value each time. However, due to the fact that there are a lot of nodes in the crystal lattice of the molecule and they are more or less evenly distributed over the spherical shell, the values of the moments of inertia will be close to each other. In this connection, it is convenient to find some average value of these moments of inertia, assuming that the mass of the molecule is uniformly distributed over the mentioned spherical shell (a smoothed fullerene).

The moment of inertia of a homogeneous spherical shell is known:

$$I=\frac{2}{3}M{a}^2.$$

(34)

Here $M$ is the mass of the fullerene; $a$ is its radius.

Formula (25) gives an approximate value of axial moments, namely:

$$A=B=C=I.$$

(35)

However, due to the fact that in the course of calculations the position of the instantaneous axis of rotation of the fullerene is known, it is not difficult to directly apply formula (33). In this case, the values will differ by 2% to 3%, with these differences being fluctuations around the mean value defined by formula (34).

2.5. The Mechanical Energy Balance of the Central Fullerene

In the absence of a magnetic field for each fullerene inside the material the energy integral is valid:

$$\frac{M{v}_C{}^2}{2}+\frac{1}{2}\left(A{p}^2+B{q}^2+C{r}^2\right)+{\displaystyle \sum _{j=1}^{N_p}{\displaystyle \sum _{k=1}^{S}U\left(r_{j,k}\right)}}=0.$$

(36)

Here A, B, C are the main central moments of inertia of the rotating fullerene. This relation is not necessary in the mathematical formulation of the problem; however, using Equation (36) allows maintaining accuracy of calculations.

3. Results

In order to calculate the inertial motion of a single internal fullerene in the fullerite material, it is necessary to correctly calculate or determine the dynamic state of the nearest environment of the selected crystal lattice site. In our work, this problem is solved by a minimal fragment of the material in which the state of fullerenes of the first and second environment echelons is determined by the solutions of the Cauchy problems formulated in exactly the same way as in the case of the central fullerene. Additionally, the extreme fullerenes of the considered particle of the material are in asymmetric conditions if their own environment is taken into account. In this regard, we set the same character of their motion as that of the central fullerene, thereby increasing the degree of symmetry of the dynamic state of the material. Along with that, this method allows solving the problem of integrity of the minimum fragment. When using the symmetry conditions in this form, the average time positions of the centers of mass of all calculated fullerenes remain unchanged. In this way, the eight fullerenes of the third echelon, which can be considered extreme in the minimum fragment, follow the symmetry conditions implying that deviations from the average positions of their mass centers are the same as those of the central fullerene and, moreover, they have the same Euler angles. As a result, only nineteen fullerenes are relatively free. In this case, carbon atoms in all fullerenes perform thermal vibrations, regardless of the nature of fullerene convergence. Simulation of free fullerene rotations, according to the mathematical model presented here, led to the following calculation frequency: $0.34\times {10}^{12}\mathrm{rad}\cdot {\mathrm{s}}^{-1}$, which is consistent with experimental data obtained using nuclear magnetic resonance [4].

At the initial moment of time, all fullerenes of the minimal fragment were at rest while being in their middle positions. At the same time, carbon atoms had harmonic vibrations around their own equilibrium positions with a frequency of $3\times {10}^{13}{\mathrm{s}}^{-1}$ and an amplitude of $\left(0.05\right)\times \left(0.142\right)\mathrm{nm}$. The calculations show that rotations of C₆₀ molecules are generated regardless of the presence of vibrations of carbon atoms. However, vibrations of fullerenes are completely determined by these vibrations. Calculations show that, if the plastic phase existed at temperatures below critical (T = 260 K), rotations would still take place. However, in reality, as the temperature decreases, the very crystal structure of the fullerite material experiences changes (rearrangement of the first kind). All calculations were performed using the standard high-order scheme with an unchanged time step of $\Delta t={10}^{-7}$ ns. Below are the results of calculations for rotation angles of the central fullerene with respect to the axes of the fixed reference frame, as well as the trajectory of a separately selected carbon atom on the visualizing sphere.

Figure 4, Figure 5 and Figure 6 show the values of Euler angles determining the instantaneous orientation of the fullerene relative to the absolute basis. According to the figures, there are areas of approximately constant values of Euler angles. This suggests local stabilization of the instantaneous axis of rotation of the fullerene around which angular oscillations of the fullerene take place. After that, the direction of the rotation axis changes, i.e., there occurs reorientation of the fullerene. In Figure 4, Figure 5 and Figure 6 these are the areas of significant changes in the angles. The length of such sections corresponds to the reorientation time of C_60s. Figure 7 also shows the trajectory of one carbon atom of the central fullerene. This trajectory shows the character of rotational motion of the fullerene. From the distributions, we can conclude that the instantaneous axis of rotation is constantly changing its direction. At the same time, in the process of motion, it is possible to select the states when the positions of the instantaneous axis stabilize, and the fullerene makes one or several rotations around it. Generally, the reasoning performed does not reveal any define regularity. Therefore, the process of rotating C₆₀ in fullerite is taken as stochastic.

In the presence of a magnetic field, there occur consistent patterns in the rotational motions of fullerenes. The angular vibration frequency behaves in a regular way (Figure 8). The atomic trajectories contain stable turns (Figure 9). Additionally, the vibration frequency of fullerenes has an average value (Figure 10) of 112 m·s⁻¹, which corresponds to the temperature of 355 K.

Instantaneous distribution of the angular frequency of rotations ω shown in Figure 8 is determined by the formula $\omega =\sqrt{p^2+q^2+r^2}.$ Averaging the values given above, in the selected time interval, we obtain the value ${\omega }_{}=3.4\times {10}^{11}\mathrm{rad}\cdot {\mathrm{s}}^{-1}$, which corresponds to the equilibrium temperature T = 283 K.

The vibration amplitude of C₆₀ depends on the amplitude and the vibrational frequency of carbon atoms in the molecule and weakly depends on fullerene rotations. Thus, if a fullerite crystal is affected by a variable magnetic field of a certain configuration, it is possible to obtain dynamic equilibrium of the system of fullerenes constituting the material under consideration so that the vibrational temperature is equal to 355 K, and the rotational temperature is 22080 K with a rotating magnetic induction vector with an amplitude of 6.7 T and a rotation frequency of 100 GHz. It is known that fullerenes in plastic fullerite rotate freely, which essentially means that the system in this case does not exhibit friction [33,34,35]. For bodies in a vacuum, the situation when the kinetic energy of translational motion differs significantly from the energy of rotations is typical (pulsars, neutron stars). Therefore, when the fullerite crystal is in vacuum and there are no adsorbed molecules or atoms in it, it will not be a big problem to obtain a high temperature of rotation while maintaining the integrity of the material. This means that solid fullerite can be considered as a material which is capable of storing energy at forced degrees of freedom.

4. Discussion

When assessing the influence of the quantum nature of the problem, we will determine the de Broglie wavelength of the fullerene in the fullerite, which depends on the mass and the velocity of the C₆₀ molecule. The average velocity of motion of the fullerene is determined from the calculated data on vibrations frequency of the fullerene, as well as available experimental data on fullerene rotation. Introducing the found velocity in the ratio for the de Broglie wavelength, we obtain that the wavelength of the fullerene molecule is substantially less than C₆₀ diameter. Thus, to describe the phenomenon of rotation for these particles, the approach implemented in the framework of classical mechanics models is fully justified.

The desire to carry out calculations in the most economical way led to the number of fullerenes in the considered area of fullerite. For this purpose the minimal element of a fullerite, which in the plastic phase contains twenty-seven fullerenes, was selected. Since the representative particle of the material belongs to a structure with a high degree of the element arrangement order, it is necessary to use periodicity conditions in the calculations. These conditions consist in the fact that the eight fullerenes of the third echelon of the environment (the edging fullerenes of the minimal fragment) have the same changes in positions of mass centers and Euler angles as those of the central fullerene of the selected fragment. In the presented model of non-deformable carbon cells, where vibrations of carbon atoms are still present, only interactions of atoms belonging to different buckminsterfullerenes are important. These cells have a random orientation relative to each other. Therefore, the effective interaction potential should be simple and, consequently, symmetrical. We used the Lennard–Jones potential, selecting the “carbon–carbon” interaction energy so as to obtain agreement with the experimental data on fullerene rotations in the crystal under consideration.

The calculations carried out within the framework of the classical approach suggest that to make the nature of rotation of the central fullerene not dependent on the fullerenes additionally included in the calculation procedure, the minimum fragment of the plastic phase of fullerite must have a symmetrical shape and contain at least twenty-seven fullerenes. Calculations show that in order for the rotation of fullerenes to have the observed frequency in the experiment, the minimum fragment of the plastic phase of fullerite must have a symmetrical structure and contain at least twenty-seven fullerenes. Thus, each fullerene can be considered as a molecular pendulum, which does not have one material element in the function of its weight but sixty elements—carbon atoms of C₆₀. These sixty atoms generate the forces of Van der Waals cross-effects so that it is these forces that cause oscillations of such pendulums. If at the initial moment of time the system is squinted, the gravitational pendulum begins to move from the maximum potential energy gradually reducing the value of this energy. In the case of a molecular pendulum, the system of a single fullerene and its immediate environment also tends to minimize the potential energy of cross-effects. Therefore, within the framework of this analogy, the root cause of rotations is the presence of an initial amplitude turn of the fullerene cell in relation to the positions of atoms belonging to other fullerenes of the nearest environment. The calculations confirmed the independence of the fullerene rotations on their own vibrations, as well as on vibrations of carbon atoms inside C₆₀ molecules. This makes it possible to spin fullerenes and accumulate a significant amount of energy inside the material. If fullerenes are encapsulated by nickel, cobalt or iron, they have a magnetic moment. Carrying out formal calculations for fullerenes having ferromagnetic properties, we obtain that in an alternating high-frequency magnetic field the nodes of the crystal lattice acquire oriented rotations with an angular velocity which is three times higher than the inertial velocity in ordinary fullerite. In addition, the internal energy and the rotational temperature increase proportionally with respect to the square of the angular velocity. Moreover, all the considered vibrations remain about the same as in the case of an unexcited crystal.

On the basis of the developed model, the task of [38] on rotations of fullerene chains in the structures of single-walled nanotubes, the problem formulated by [11] on rotations of the inner shell of the C₂₀@C₆₀ nanoparticle, the problem stated by [19] on nested fullerene cells, and many other tasks can be considered.

5. Conclusions

An economical and affordable model of molecular interaction inside fullerite was developed. Within the framework of the model, it is shown that using variable electromagnetic fields it is possible to spin the nodes of the crystal under consideration up to certain speeds at which the vibrational temperature is not too high. This makes it possible to increase both the hardness of the crystal and its strength with respect to effects of the impact type. The field configuration can determine the nature of fullerene rotations; for example, it can orient the axis of fullerene rotations throughout the whole fullerite in view of the unique structure of the material, while other properties of the material, such as electrical conductivity, thermal conductivity and others, are changed as well. The mathematical model which is quite universal and is limited only by the applicability of the classical mechanics approach is described in detail. It requires that the potential correctly describes the structure under consideration, and its parameters are well adapted to the elements of the structure. Moreover, it is quite possible to change one potential to another or combine them, depending on the conditions of the problem.