# Semi-quantum molecular dynamics simulation of thermal properties and heat transport in low-dimensional nanostructures

## Abstract

We present a detailed description of semi-quantum molecular dynamics simulation of stochastic dynamics of a system of interacting particles. Within this approach, the dynamics of the system is described with the use of classical Newtonian equations of motion in which the effects of phonon quantum statistics are introduced through random Langevin-like forces with a specific power spectral density (the color noise). The color noise describes the interaction of the molecular system with the thermostat. We apply this technique to the simulation of thermal properties and heat transport in different low-dimensional nanostructures. We describe the determination of temperature in quantum lattice systems, to which the equipartition limit is not applied. We show that one can determine the temperature of such system from the measured power spectrum and temperature- and relaxation-rate-independent density of vibrational (phonon) states. We simulate the specific heat and heat transport in carbon nanotubes, as well as the heat transport in molecular nanoribbons with perfect (atomically smooth) and rough (porous) edges, and in nanoribbons with strongly anharmonic periodic interatomic potentials. We show that the effects of quantum statistics of phonons are essential for the carbon nanotube in the whole temperature range K, in which the values of the specific heat and thermal conductivity of the nanotube are considerably less than that obtained within the description based on classical statistics of phonons. This conclusion is also applicable to other carbon-based materials and systems with high Debye temperature like graphene, graphene nanoribbons, fullerene, diamond, diamond nanowires etc. We show that the existence of rough edges and quantum statistics of phonons change drastically the low temperature thermal conductivity of the nanoribbon in comparison with that of the nanoribbon with perfect edges and classical phonon dynamics and statistics. The semi-quantum molecular dynamics approach allows us to model the transition in the rough-edge nanoribbons from the thermal-insulator-like behavior at high temperature, when the thermal conductivity decreases with the conductor length, to the ballistic-conductor-like behavior at low temperature, when the thermal conductivity increases with the conductor length. We also show how the combination of strong nonlinearity of periodic interatomic potentials with the quantum statistics of phonons changes completely the low-temperature thermal conductivity of the system.

## I Introduction

Molecular dynamics (MD) is a method of numerical modeling of molecular systems based on classical Newtonian mechanics. It does not allow for the description of pure quantum effects such as freezing out of high-frequency oscillations at low temperatures and the related decrease to zero of heat capacity for . In classical molecular dynamics (CMD), each dynamical degree of freedom possesses the same energy , where is Boltzmann constant. Therefore, in classical statistics the specific heat of a solid almost does not depend on temperature when only relatively small changes, caused the anharmonicity of the potential, can be taken into account (1). On the other hand, because of its complexity, a pure quantum-mechanical description does not allow in general the detailed modeling of the dynamics of many-body systems. To overcome these obstacles, different semiclassical methods, which allow to take into an account quantum effects in the dynamics of molecular systems, have been proposed (2); (3); (4); (5); (6); (7); (8). The most convenient for the numerical modeling is the use of the Langevin equations with color-noise random forces (5); (7). In this approximation, the dynamics of the system is described with the use of classical Newtonian equations of motion, while the quantum effects are introduced through random Langevin-like forces with a specific power spectral density (the color noise), which describe the interaction of the molecular system with the thermostat. Below we give a detailed description of this semi-quantum molecular dynamics (SQMD) approach in application to the simulation of specific heat and heat transport in different low-dimensional nanostructures.

The paper is organized as follows. In Section II we describe the temperature-dependent Langevin dynamics of the system under the action of random forces. If the random forces are delta-correlated in time domain, this corresponds to the white noise with a flat power spectral density. This situation corresponds to high enough temperatures, when is larger that the quantum of the highest-frequency mode in the system, . But for low enough temperature, , the stochastic dynamics of the system should be described with the use of random Langevin-like forces with a non-flat power spectral density, which corresponds to the system with color noise. In Section III we describe the determination of temperature in quantum lattice systems, to which the equipartition limit is not applied. We show how one can determine the temperature of such system from the measured power spectrum and temperature- and relaxation-rate-independent density of vibrational (phonon) states. In Section IV we describe a method for the generation of color noise with the power spectrum, which is consistent with the quantum fluctuation-dissipation theorem (1); (9). In Sections V and VII we apply such semi-quantum molecular dynamics approach to the simulation of the specific heat and phonon heat transport in carbon nanotubes. In Section VI we apply this method to the simulation of heat transport in a molecular nanoribbon with rough edges. Previously, we have predicted analytically (10) and confirmed by classical molecular dynamics simulations (11) that rough edges of molecular nanoribbon (or nanowire) cause strong suppression of phonon thermal conductivity due to strong momentum-nonconserving phonon scattering. Here we show that the rough edges and quantum statistics of phonons change drastically the thermal conductivity of the rough-edge nanoribbon in comparison with that of the nanoribbon with perfect (atomically smooth) edges. The semi-quantum molecular dynamics approach allows us to model the transition in the rough-edge nanoribbons from the thermal-insulator-like behavior at high temperature, when the thermal conductivity decreases with the conductor length, see Ref. (11), to the ballistic-conductor-like behavior at low temperature, when the thermal conductivity increases with the conductor length. In Section VIII we apply this technique to the simulation of thermal transport in a nanoribbon with strongly anharmonic periodic interatomic potentials. We show how the combination of strong nonlinearity of the interatomic potentials with quantum statistics of phonons changes completely the low-temperature thermal conductivity of the nanoribbon. In Section IX we provide a summary and discussions of all the main results of the paper.

## Ii Langevin equations with color noise

In the presence of a Langevin thermostat, equations of motion of coupled particles have the following form:

(1) |

where the three-dimensional radius-vector gives the position of the -th particle at the time instant (, being the total number of particles), is the particle mass, is the Hamiltonian of the system and gives the force applied to the -th particle caused by the interaction with the other particles, is the friction coefficient ( being the relaxation time due to the interaction with the thermostat), and are random forces with the Gaussian distribution.

In the semi-quantum approach the random forces do not represent in general white noise. The power spectral density of the random forces in that description should be given by the quantum fluctuation-dissipation theorem (1); (9):

(2) |

where , are Cartesian components, and are Kronecker delta-symbols. Here the positive and even in frequency function

(3) | |||

gives the dimensionless power spectral density at temperature of an oscillator with frequency . The first term in Eq. (3) gives the contribution of the zero-point oscillations to the power spectral density of random forces. With the use of Eqs. (2) and (3), we get the correlation function of the random forces as

(4) |

Let be the maximal vibrational eigenfrequency of the molecular system. At high temperatures, , one has and the molecular system will perceive the random forces as a delta-correlated white noise:

(5) |

Therefore for high enough temperatures we deal with molecular dynamics with Langevin uncorrelated random forces. But for low temperatures, the random forces are correlated according to the function given by Eqs. (3) and (4). In other words, for low temperatures the random forces in the system represent the color noise with the dimensionless power spectral density given by Eq. (3).

## Iii Determination of temperature in quantum lattice systems

In connection with the consideration of the thermodynamics and dynamics of non-classical lattices, in this Section we discuss possible determination of the temperature in the systems, to which the equipartition limit is not applied. We start from Eqs. (1), rewritten as the equations of motion for the lattice vibrational mode with eigenfrequency :

(6) |

where are random forces with the spectral density given by Eqs. (2) and (3).

and power spectrum of lattice vibrations with nonzero eigenmodes:

(8) |

We can use and the power spectrum , Eqs. (7) and (8), to compute the average energy of the system. Taking into account that the above quantities are even functions of frequency , for positive the average energy of the system can be computed as:

(9) |

In the weakly dissipative limit, with for all , which is realized in all the considered system, for positive and one has, in the sense of generalized functions:

(10) |

We can introduce in this limit the density of vibrational (phonon) states (for positive frequencies) as:

(11) |

Then the power spectrum (for positive ) consists of an array of the delta-functions at the system eigenfrequencies, weighted by temperature,

(12) | |||||

and the average energy of the system (9) has the following form:

(13) |

Here the first term in the sum gives the temperature-independent contribution of zero-point oscillations to the energy of the system. As one can see from Eqs. (11)-(13), in the weakly dissipative limit the density of vibrational states, power spectrum and average energy of the system do not depend explicitly on the relaxation rate .

From Eq. (12) we get the ”quantum definition” of temperature through the power spectrum of lattice vibrations at the given frequency:

(14) |

where

(15) |

In the case when one omits the zero-point contribution to the spectrum of the random forces, namely considers

(16) |

in Eqs. (2) and (6), the temperature can be determined from the following equation:

(17) |

where in the definition of the function in Eq. (15) the function , in contrast to function in Eq. (8), does not have now the zero-point contribution,

(18) |

and, correspondingly, one has

(19) |

In the equipartition limit, which is realized for , one has , , and both Eqs. (14) and (17) turn into the identity .

Therefore one can determine the temperature of the quantum lattice system from the measured power spectrum and temperature- and relaxation-rate-independent density of vibrational (phonon) states . In Fig. 6 we show the power spectrum, computed within the semi-quantum approach with , ps, line 1, and in the limit, when the spectrum is given by an array of the smoothed delta-functions, line 3. As one can see in this figure, the particular choice of the (long enough) relaxation time indeed does not affect the presented results. The determination of temperature, given by Eqs. (14) and (17), is valid for weakly anharmonic systems, in which the power spectrum of lattice vibrations is close to the one, given by Eq. (8). Carbon-based materials like carbon nanotubes, graphene and graphene nanoribbons, and crystal structures with stiff valence bonds belong to such systems. It is worth mentioning that the temperature, given by Eq. (14) or Eq. (17), does not depend on frequency only in the case of equilibrium harmonic lattice excitations, driven by random forces with power spectrum given by the Bose-Einstein distribution, Eqs. (2) and (3). The anharmonicity of the lattice can change the power spectrum of the excitations, see Fig. 12 (b) below. In such case the ”phonon temperature” starts to depend on frequency (similar to the frequency- and direction-dependent temperature of photons in non-equilibrium radiation, see Ref. (1)), and ”hot phonons” with non-equilibrium high frequencies appear in the system, see Section VI below.

## Iv Color noise generation

For the implementation of the semi-quantum approach in molecular dynamics, random forces with the power spectral density given by must be generated. Existing numerical techniques allow the generation of a random quantity from a prescribed correlation function (12). This approximation was used in (7) in the modeling of thermodynamic properties of liquid He above the point. But universal methods require in general the usage of high numerical resources. For instance, the utilization of the method proposed in (12) requires the execution of the Fast Fourier Transform for every random force value. We propose a technique which takes into account the specific form of the dimensionless power spectral density , where is a dimensionless frequency.

In terms of , the dimensionless power spectral density of random forces looks as

(20) |

The random force with dimensionless power spectral density, given by Eq. (20), can be presented as a sum of two independent functions, , where the first function, , has a linear power spectral density, while the second one is .

To generate the color noise with a given dimensionless power spectral density , we will construct the dimensionless random vector functions of the dimensionless time , which will give the power spectral density of the random forces in Eqs. (1) as

(21) |

such that

(22) |

and

(23) |

Here, the uncorrelated random functions and will generate, in a finite frequency interval, the power spectra and , respectively.

For a molecular system with maximal vibrational eigenfrequency , we will construct the dimensionless random vector function , whose power spectral density generates in the frequency interval (). Each component of this vector can be written as a sum of dimensionless random functions

(24) |

The functions which generate the color noise, are the solutions of the linear equations:

(25) |

where is the derivative with respect to . Here are dimensionless white-noise random forces, with the correlation functions

(26) |

By solving Eqs. (25) in the frequency domain, can be substituted in Eq. (24) and the power spectra of can be obtained as

(27) |

The dimensionless parameters and can be found by minimizing the integral

(28) |

with respect to the parameters , where , . The coefficients , , , obtained within this procedure, are shown in Table 1. For these parameters, integral (28) reaches its lower value, equal to .

coefficient | value | coefficient | value |
---|---|---|---|

1.763817 | 1.8315 | ||

0.394613 | 0.3429 | ||

0.103506 | 2.7189 | ||

0.015873 | 1.2223 | ||

1.043576 | 5.0142 | ||

0.177222 | 3.2974 | ||

0.050319 | |||

0.010241 |

As one can see from Fig. 1, the function given by Eq. (27) approximates with high accuracy the linear function in the frequency interval .

On the other hand, the random function , which will generate the power spectral density , can be approximated by a sum of two random functions with relatively narrow frequency spectra:

(29) |

In this sum the dimensionless random functions , , satisfy the equations of motion as

(30) |

where, as before, are delta-correlated white noise functions:

(31) |

We can solve, as previously, Eq. (30) in the frequency domain and insert into Eq. (29) to obtain the power spectrum of :

(32) |

The dimensionless parameters can be found by minimizing the integral

(33) |

with respect to the parameters . The obtained parameters are shown in Table 1. For these parameters, integral (33) reaches its lower value, equal to . Since the integral (33) is determined by the rapid-descending functions, the infinite integration limit can be replaced by a finite one . Numerical integration of the integral (33) shows that its minimal value practically does not depend on the upper integration limit for .

As one can see from Fig. 2, the function, given by Eq. (32), approximates with high accuracy the function in the frequency interval .

Therefore in the semi-quantum molecular dynamics approach one has to solve the Langevin equations (1) with random forces , whose power spectral density is determined as:

(34) |

From Eq. (22) we get the relation between the dimension and dimensionless random forces: .

In Fig. 3 we show the comparison of the power spectral density for the dimensionless frequency , given by the sum of Eqs. (27) and (32), with the function , Eq. (20). We see a very good coincidence in the frequency interval for .

To describe the dynamics of molecular system with an account for quantum statistics of molecular vibrations but without an account for zero-point oscillations, one has to keep only the last two terms, with , in Eq. (34) and to solve numerically Eqs. (30) in order to simulate random forces in Eq. (1). Equations (30) play the role of the filter, which filters out the dimensionless high-frequency (”non-quantum”) component of the white noise at low temperature.

It is worth mentioning in this connection that the considered semi-quantum approach permits one to obtain the correct value of the energy of zero-point oscillations. A semi-quantum account for zero-point oscillations in the modeling of the properties of liquid He above the point has allowed to describe correctly the liquid state of the helium (7), while the classical molecular dynamics description (without an account for zero-point oscillations) predicts the solid state of the helium at the same low temperature.

## V Temperature dependence of specific heat of carbon nanotube

For the illustration of the implementation of the semi-quantum approach in molecular dynamics, in this Section we find the temperature dependence of specific heat of a single-walled carbon nanotube. We consider the nanotube with the (5,5) chirality index (the zig-zag structure), which consists of carbon atoms, see Fig. 4. In the following we will use the molecular-dynamics model of the carbon nanotube, which was discussed in details in Ref. (13).

To model nanotube dynamics at temperature in the classical approximation, we will solve the Langevin equation (1) with a white noise, where is a Hamiltonian of the nanotube and all , being a mass of a carbon atom. The correlation functions of the random forces , which describe the interaction of an -th atom with the thermostat, are given by Eq. (5).

We take the initial conditions, which correspond to the equilibrium stationary state of the nanotube at , and integrate a system of equations of motion (1) for the time , during which the system will reach the equilibrium state at finite temperature. The integration beyond this time will allow us to find the average thermal energy of the system . Then the specific heat of the molecular system can be found from the temperature dependence of the average thermal energy . As one can see in Fig. 5, the average thermal energy of the carbon nanotube, placed in classical Langevin thermostat, is strictly proportional to temperature, , line 1, and, correspondingly, nanotube specific heat almost does not depend on temperature, line 3. This result shows both the correctness of the modeling and that carbon nanotubes are stiff structures which possess weak anharmonicity of the dynamics of the constituting atoms. One can also see in Fig. 5 that the average thermal energy per mode without an account for zero-point oscillations has the feature that in general and for temperature much less than the Debye one. For instance, in the carbon nanotube is almost six times less than at room temperature =300 K.

To obtain power spectra of the thermal oscillations of the atoms in the carbon nanotube, one has to switch off the interaction with the thermostat after the establishment of the thermal equilibrium in the system. In other words, one has to solve numerically the frictionless equations of motion, Eqs. (1) with and , with the initial conditions for and , which correspond to the thermalized state of the molecular system. By performing the Fast Fourier Transform of the time-dependent particle velocities , one can get the power spectrum , Eq. (8). To increase the accuracy of the measurement, it is necessary to perform the averaging of the obtained results on different initial thermalized states of the system. Power spectral density of the thermal vibrations of the carbon nanotube atoms is shown in Fig. 6. As one can see from this figure, all vibrational modes of the system are excited in the classical Langevin thermostat.

Carbon nanotube is a rigid structure in which the anharmonicity of atomic dynamics at room temperature is weakly pronounced: as one can see from line 2 in Fig. 6, the characteristic Debye temperature of the carbon nanotube is very high, K, where cm is the maximal phonon frequency in the system. Therefore we can obtain the temperature dependence of the carbon nanotube specific heat with the use of the spectrum of the harmonic (non-interacting) vibrational eigenmodes: , where the average energy of the system is given by Eq. (13). To obtain the eigenfrequencies , we have to find all the eigenvalues of the square symmetric matrix of the size: if is one of the eigenvalues, the eigenfrequency is determined as . For atoms in the nanotube, we have eigenmodes, from which the first 6 modes, which describe rigid translations and rotations of the nanotube, have zero frequency, , while the other eigenfrequencies are nonzero and are distributed in the interval 25.1 cm cm, . We can obtain the temperature dependence of the nanotube specific heat with the use of Eq. (13):

(35) |

As one can see in Fig. 5, the nanotube specific heat in the quantum description monotonously goes to zero for , and the effects of quantum statistics of phonons are essential for the carbon nanotube in the whole temperature range K, in which the specific heat of the nanotube is considerably less than the classical-statistics value . This conclusion also applies to other carbon-based materials and systems with high Debye temperature like graphene, graphene nanoribbons, fullerene, diamond, diamond nanowires etc.

To model the nanotube stochastic dynamics in the semi-quantum approach, we will use the Langevin equations of motion (1) with random forces with the power spectral density, given by , see Eq. (34). Here it is explicitly taken into account that zero-point oscillations do not contribute to the specific heat of the system, and therefore one needs to account only for the random functions, given by Eq. (32).

Therefore the semi-quantum approach in this case amounts to the integration of Eqs. (1) and (30), to find the average energy and the specific heat of the system. As one can see from Fig. 5, the temperature dependence of the carbon nanotube specific heat, obtained by studying the nanotube stochastic dynamics with color noise, coincides completely with the dependence, obtained from the spectrum of the harmonic vibrational eigenmodes in the system. Power spectral density of carbon atoms vibrations, see Fig. 6, shows that the semi-quantum approach describes the thermalization of only the low-frequency eigenmodes, with . In such a way, the approach models the quantum freezing of the high-frequency eigenmodes at low temperature . For the carbon nanotube, the classical description is definitely not valid at room temperature K, when the characteristic frequency cm is much less than the maximal one in the system, cm, see Fig. 6.

It is worth mentioning in this connection that the obtained results for the temperature dependence of the specific heat of the carbon nanotube C almost do not depend on the relaxation time . All the five values of , , 0.2, 0.4, 0.8, 1.6 ps, give the same average energy and specific heat of the nanotube. One cannot use either too short , when the system dynamics becomes the forced one, or too long , which increases the computation time. We find that the optimal relaxation time for the carbon nanotube is ps.

Above we have introduced and discussed the lattice temperature of the system embedded in the Langevin thermostat with color noise determined by the quantum fluctuation-dissipation theorem, see Eqs. (14) and (17). But to compute the thermal conductivity, one also needs to determine the local temperature of the part of the molecular system, which is placed between two thermostats with different temperatures. In the classical-statistics approximation, the lattice temperature can be determined as the average double kinetic energy per mode (or per degree of freedom):

(36) |

where is a number of atoms and is assumed that . As we have shown in Chapter III, such determination of temperature does not work in the semi-quantum approximation when one needs to deal with the density and population of vibrational (phonon) states. In this case the determination of temperature is given by Eq. (14) or by Eq. (17). But in the MD simulations, one can also use the determination of temperature, which is equivalent to the integral presentation of Eq. (12).

Having in mind further simulation of thermal conductivity, we omit in Eqs. (2) and (6) the zero-point contribution to the spectrum of the random forces. With the use of Eqs. (9), (12), (18) and (19), we find the average energy per mode (per degree of freedom) as

(37) |

where is density of vibrational states, cf. Eq. (11).

From the MD simulations, we can also measure the spectral density of the average double kinetic energy per mode as

(38) |

In the case of correct ”quantum” level occupation, the both Eqs. (37) and (38) should correspond to the same value of temperature, therefore we obtain the following equation for the (numerical) determination of the lattice temperature:

(39) |

Since function monotonously increases with , Eq. (39) has a unique solution for the temperature. Similar equation, but with an account for zero-point oscillations, for the rescaling of MD temperature for an effective one was used in Ref. (14) for quantum corrections of MD results. As we have mentioned above in discussion of Fig. 5, the neglect of zero-point oscillations in the computation of the average thermal energy leads to the feature that in general and for temperature much less than the Debye one, see also Fig. 11 (b) and (c) below.

Essentially one can use Eq. (39) for the determination of the lattice temperature only for the quantum-statistics population of phonon states. But the situation can emerge when high-frequency phonons have an excess excitation, in comparison with the Bose-Einstein distribution for the given temperature, see section VII below. In this case the lattice temperature can be determined with the use of only part of the phonon spectrum:

(40) |

where is the frequency interval with the Bose-Einstein distribution of the average energy per mode . For , and , Eq. (40) gives the classical definition of temperature, . Clearly Eq. (39), which is more convenient for the numerical modeling, gives the correct value of lattice temperature only for and .

If the molecular system is completely embedded in the Langevin thermostat with color-noise random forces, their spectrum, given by Eq. (3), provides the correct ”quantum” population of all the phonon states. As one can see in Fig. 6, at K the spectral density of the average energy per mode in the carbon nanotube (red line 1) coincides exactly with the smoothed function (black line 3). In this case Eq. (39) determines the temperature, which is exactly equal to the one of the thermostat.

## Vi Thermal conductivity of nanoribbon with rough edges

First we apply the semi-quantum approach for the molecular dynamics simulation of thermal conductivity in the nanoribbon with rough edges and harmonic interparticle potential. In such system the vibrational eigenmodes do not interact and the considered semi-quantum approach turns to be the exact one.

It was predicted analytically and confirmed by classical molecular dynamics simulations that rough edges of molecular nanoribbon (or nanowire) cause strong suppression of phonon thermal conductivity due to strong momentum-nonconserving phonon scattering (10); (11). In the case of nonlinear interatomic potential, the ribbon has a finite, length-independent thermal conductivity, while in the case of harmonic interatomic potential, the thermal conductivity decreases with the ribbon length (and the ribbon behaves as a thermal insulator). Essentially for both interatomic potentials, the thermal conductivity increases with the length of the nanoribbon with perfect (atomically smooth) edges. In the case of harmonic interatomic potential, phonons in the nanoribbon with rough edges experience effective localization due to strong antiresonance multichannel reflection from side atomic oscillators in the rough edge layers (the Anderson-Fano-like localization due to interference effects in phonon backscattering), see (10); (11). Apparently such strong backscattering can localize phonons with the wavelength shorter than the ribbon length. Within the classical description, such effect does not depend on temperature. But this picture will be changed if we take into account the quantum effects. The latter result in thermalization and correspondingly in participation in low temperature thermal transport of long wave phonons only. The long wave phonons are much less affected by surface roughness than the short wave phonons and therefore the effect of thermal conductivity suppression should decrease with the decrease of temperature. Below we demonstrate this effect within the proposed semi-quantum approach.

We consider the system which consists of parallel molecular chains in one plane (11). Let be the chain number, be the molecular number in the chain, then the equilibrium position of the -th atom in the -th chain will be , , where and are, respectively, the intra- and interchain spatial periods, see Fig. 7.

Only the longitudinal displacements are taken into account in a simplest scalar model of two-dimensional crystal: , , where is a dimensionless longitudinal displacement of the -th molecule from its equilibrium position. Hamiltonian of the system has the form as

(41) | |||||

where is a number of molecules in each chain, is a mass of molecule, is a characteristic interaction energy, is the dimensionless potential of the nearest-neighbor intrachain interactions, describes relative displacements of the nearest-neighbor atoms, describes the interchain interaction.

Dimensionless energy of the nearest-neighbor interchain interactions for odd is

(42) |

and for even is

(43) |

where is the dimensionless potential of the nearest-neighbor interchain interaction.

We will consider a finite rectangular with free edges, and harmonic intra- and interchain interaction potentials:

(44) |

For the convenience of the modeling, we introduce the dimensionless quantities: temperature , time , and energy , where , , and . Then the dimensionless Hamiltonian will have the form as

(45) |

where the dot denotes the derivative with respect to the dimensionless time .

For the dimensionless quantities, equations of motion will take the form: