{'fraction_non_alphanumeric': 0.025887539276837217, 'fraction_numerical': 0.010598883834357267, 'mean_word_length': 4.741518578352181, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Calculation of Thermal Conductivity from Einstein Relation in Equilibrium MD

In equilibrium molecular dynamics, Einstein relation can be used to calculate the thermal conductivity. This method is equivalent to Green-Kubo relation and it does not require a derivation of an analytical form for the heat current. However, it is not commonly used as Green-Kubo relationship. Its wide use is hindered by the lack of a proper definition for integrated heat current (energy moment) under periodic boundary conditions. In this paper, we developed an appropriate definition for integrated heat current to calculate thermal conductivity of solids under periodic conditions. We applied this method to solid argon and silicon based systems; compared and contrasted with the Green-Kubo approach.

In equilibrium molecular dynamics, Einstein relation can be used to calculate the thermal conductivity. This method is equivalent to Green-Kubo relation and it does not require a derivation of an analytical form for the heat current. However, it is not commonly used as Green-Kubo relationship. Its wide use is hindered by the lack of a proper definition for integrated heat current (energy moment) under periodic boundary conditions. In this paper, we developed an appropriate definition for integrated heat current to calculate thermal conductivity of solids under periodic conditions. We applied this method to solid argon and silicon based systems; compared and contrasted with the Green-Kubo approach.

I. INTRODUCTION
Accurate determination of thermal conductivity of materials is of fundamental importance especially for technological applications where the control of thermal transport is critical.
For instance, microelectronics, where the device failures are mostly due to overheating, require a strict thermal management. The development of efficient thermoelectrics require materials and materials systems with low thermal conductivity. A similar requirement arise in the development of thermal barrier coatings. Experimental evaluation of thermal conductivity is not as straightforward as electrical conductivity. Measurements generally suffer from heat losses, non-uniform heating and errors introduced due to approximations that account for sample size and structure 1 . As the dimensions of the samples get smaller, these effects become more pronounced, increase the inaccuracy of the measurements. Due to their wide use in electronics, thermal transport in silicon and silicon based alloys has been studied extensively. Within the last decade, silicon based materials have been considered for thermoelectric and thermal barrier applications where a small thermal conductivity is desired 2 . The studies in this area generally focus on the influence of defects, impurities and nano-structure (superlattices 3 , nanowires 4 , quantum dots 5 , etc.) on thermal transport.
As an alternative to experiments, computational approaches offer controllable virtual measurements to determine thermal conductivity, with the added potential of probing the processes and mechanisms that influence the thermal transport. The two most frequently used computational methods that include atomistic level details in studying thermal transport are based on Boltzmann transport equation (BTE) and molecular dynamics (MD) simulations. The first one, BTE, is a semiclassical approach. The atomistic level granularity in BTE, is provided by the determination of phonon frequencies and group velocities from the lattice dynamics of an N-body system. However, the detailed physics of the heat transfer is incorporated into the formulation through phonon relaxation times [6][7][8] . This requires either experimental results or a detailed analysis of the fundamental scattering processes to provide reliable theoretical estimates of relaxation times. When such data is not available due to the complexity of the system, the relaxation times can be probed with MD. Details of such calculations are described in earlier studies where MD was used to obtain the phonon relaxation times in fcc crystals 9 and silicon 10,11 using the decay of phonon mode amplitudes or energies.
In contrast, the full atomistic detail of the physics of heat transport is inherently present in MD through the definition of interactions between the atoms that make up the system. The evolution of the state of the system is given by the equations of motion based on interaction potentials that are usually anharmonic. Hence, MD is a powerful tool, as it includes normal, Umklapp and higher order phonon-phonon scattering events through the dynamics of the thermodynamic state of the system. Of course due to the classical nature of MD, vibrational modes are equally excited (i.e. have the same energy k B T ), 12,13 and the phonon scattering processes do not include electronic contributions. However, some approaches 14,15 have also implemented quantum corrections. There are three MD based approaches for calculating thermal conductivity: the steady state non-equilibrium (direct) method, synthetic hamiltonian non-equilibrium method and the equilibrium (Green-Kubo) method. The direct method mimics experiment by imposing a temperature gradient on the system and determines the thermal conductivity from the phenomenological Fourier law [16][17][18][19][20][21] . Temperature gradients can be created by scaling the particle velocities, changing the ensemble averages or switching particle velocities between the hot and cold ends. There can be some complications with this method due to the nonlinear response of a system associated with very large temperature gradients. The second technique, proposed by Evans 22 applies a homogenous synthetic field to create a heat flux, as opposed to applying a temperature gradient. Finally, Green and Kubo 23,24 showed that through the fluctuation-dissipation theorem, transport properties i.e. the diffusion coefficient (D), shear viscosity (η µν ) and thermal conductivity (κ) can be calculated by integration of autocorrelation functions of particle velocity (v i ), pressure tensor (P µν ) and heat current (J ).
where µ and ν are the vector components, V is the system volume, k B is the Boltzmann constant, T is temperature and N is the total number of particles. The same transport coefficients may also been obtained via 'Einstein' relation; which are related to doubly integrated form of these expressions.
In the Eqs. 4, 5 and 6, r i , p iµ and R are the position and momentum of i th particle and integrated heat current respectively. The integrated heat current will be referred to as the energy moment in this paper and is defined for any N-body system as: where r i is the position and ǫ i is the energy content of the i th particle.
Detailed derivations on the time correlations functions and associated Einstein relations for transport properties can be found in several sources [25][26][27] . In the case of thermal conductivity, the heat current autocorrelation function (HCACF) approach has been employed for a number of materials such as bulk diamond 14 , bulk silicon 28 , metal organic frameworks 29 , nanotubes 30 , nanoribbons 31 and nanoparticles 32 . However, unlike the wide use of the Einstein relation in determining mass transport/diffusion, the use of the Einstein relation for thermal conductivity is rather rare. This situation has arisen from a lack of proper implementation in the calculation of integrated heat current/energy moment in simulations.
In the application of the Einstein relation based thermal conductivity calculation in MD, discrepancies arise regarding the calculation of the energy moment and energy current. The techniques used for calculating R give an expression which produces a bounded behavior and results in zero thermal conductivity in non-diffusive solid systems, while J yields a reliable result. Similar problems have been reported for the calculation of shear viscosity from the Einstein relation using MD with periodic boundary conditions 33,34 . Here, we propose an expression for R that gives a proper thermal conductivity and can be reduced numerically to give the exact value of J . Further, for solid systems, the resulting energy moment expression can be generalized to be independent of the potential model. In this report, we formulate the new method and demonstrate its use through applications to a two-body potential (Lennard-Jones) and then an N-body potential (Tersoff) 35 .
The paper is organized as follows: In the next section we present the background and reformulated form of the energy moment based method of calculating thermal conductivity.
In section III, we present the application of the method, to Argon, Silicon, and Porous Silicon and silicon Germanium nanodots in Silicon matrix, providing a detailed comparison of energy moment based method with HCACF based methods. Where we observe that HCACF of heterogeneous systems (Si-Ge) or systems with extended defects (Si with pores) are showing large fluctuations that one needs to filter out or large variance in resulting thermal conductivity values. In contrast the Einstein relation for all systems has a much better behavior. In section IV, we summarize our results and present our conclusions.

II. BACKGROUND AND THEORY
The fundamental variable for calculating thermal conductivity from Eq. 6 is the energy moment, R is defined above in Eq. 7. The energy content of particle i, ǫ i , is defined as: Here, the first and the second terms are the total kinetic and potential energy of the particle i. The velocity and the mass of the same particle are represented by v i and m i . The heat current, J , is the time derivative of R.
As the general definition of R and J are given, we will describe the method of determining thermal conductivity with two-body potential and generalize to n-body potential. For twobody interactions, the total potential energy content of a particle is defined as u i = 1 2 j =i u ij where u ij is the potential between particle i and j. Then, the corresponding heat current can be calculated through the following expression: where r ji is defined as the nearest image distance and f i ij the interatomic force on i due to j.
A straightforward implementation of Eq. 6 through Eq. 8 we call formulation 1, by considering a bulk non-convective solid in a periodic box, as we will show produces no thermal conductivity. Moreover, if r ji in Eq. 10 is the nearest image distance, then Eq. 7 cannot be reduced to Eq. 10. The issue of defining R in periodic systems has been pointed out by Donadio and Galli 36 , who attributed the inconsistency to the ill-definition of particle positions in a periodic system. However, we believe this is a consequence of the omission of what we call cross-boundary interactions, interactions or transfer of energy between the N atoms in the simulation region and the image atoms across the periodic boundary. These interactions are vital to accurately defining the energy dispersion of the system and must be included into the definition of R. With the inclusion of the cross-boundary interactions, energy moment can be separated into potential (R p ) and kinetic (R k ) contributions. The combination of R p + R k will be called formulation 2 for the rest of the report.
The new form for the potential energy portion of R, R p , is where u ij is the pair potential, i represents atoms in the box and j represents the (real or image) neighbors of i. To avoid double counting, j > i constraint is imposed. Though the form is very similar to Eq. 7, if an atom in the simulation box interacts with an image atom, the actual atomic position of the image is used, instead of contributing the image interaction to its corresponding particle in the box. Although a simple equipartition of potential energy is the reasonable choice in the case of two-body interactions, n-body interactions may require other treatments and a rule of energy distribution may not preferable to another. In these circumstances, one needs to derive a separate expression for the heat current from R for each choice of potential splitting. However, this derivation is avoided when using R p directly.
The kinetic portion of the energy moment requires a similar approach. In short, the transfer of kinetic energy from atoms in the box to image atoms must be included. This transfer is mediated by power term, f · v. The kinetic portion, R k , can then be written as: where f i ij represents the force on atom i due to j. This equation represents the transfer of kinetic energy from atom i to all particles that interact with it. Again, if the neighbor is an image atom, the image position is used. This discussion has been carried out for the two body interaction, but it can also be formatted for n-body potential. In fact, the kinetic part can be put in a form that does not depend explicitly on the potential form by simply rearranging the summation as: where l is a summation over particles (real or image) that make up unique n-body interaction groups. Though similar to the previous expression, Eq. 13 is well defined for any interatomic potential as f l is simply the force on particle l due to all unique interaction groups that include l.
The computations can be restricted to the kinetic portions of the energy moment when dealing with the non-convective systems. In solids, there is very little convection, and to see diffusion in a simulation is highly improbable. This, along with the fact that the potential contribution to the atomic energy is bound, leads to the conclusion that R p does not contribute to the thermal conductivity in perfect solid systems. Then, R k must be main component driving the increase in the Einstein relation. On the other hand, for highly diffusive systems such as liquids, R p is expected to contribute significantly.

A. Model Systems
In order to further illustrate the formulations given in the theory, we studied archetypical two-body (argon) and n-body (silicon) systems. Moreover, we discuss the behavior of Einstein relation when impurities (germanium) and pores are present in silicon.

Two-Body Potential
We calculate the lattice thermal conductivity of face centered cubic argon. We have chosen the same system size (256 atom) and equilibration procedure as in McGaughey and Kaviany 37 in order to calculate HCACF. After equilibration, the system is allowed to run in NVE ensemble at 50 K for 2 ns with a time step of 0.1 femtosecond. We are aware that this time step is small but we chose it to be consistent with the previous literature 37 and to enable us to calculate more accurate numerical derivatives of (R). We performed five independent simulations to average HCACF's. The normalized HCACF is plotted in Fig.1 which is similar to what was obtained in McGaughey and Kaviany's work. In Fig. 2 we have plotted the time evolution of R. It is clear that there is a quantitative difference between formulation 1 and formulation 2. The former, as expected, fluctuates around a constant value. Whereas the latter displays an increasing tendency.
In principle, the analytical and numerical derivatives of R should produce the same heat current (J ). Fig. 3 compares the behavior of the heat current obtained from the analytical and numerical derivative of R. The inconsistency between the formulation 1 and analytical J is obvious. On the other hand, the numerical derivatives of formulation 2 perfectly matches the analytical form.
For the calculation of thermal conductivity, the system size is enlarged to 4000 atoms and the total simulation time is extended to 20 ns. Again, five different initial conditions are evaluated and the results are averaged. Ensemble average for the autocorrelation is obtained by selecting 6 millions time origins 3 fs apart from each other. In calculation of κ from HCACF, the data was fitted to a double exponential, Eq. 14, and the infinite integration was carried out over this function, an approach adopted by many previous studies 14, 28,30,37 . f (t) = a 1 e −t/τ 1 + a 2 e −t/τ 2 In this equation a 1 , a 2 , τ 1 and τ 2 are parameters to be fitted and t is the correlation time.
Both the autocorrelation function and corresponding double exponential function are plotted in Fig. 4. The thermal conductivity is calculated as 0.403±0.01 W/mK which is very close to the value (0.417 W/mK) given in Ref. 37.
The Einstein's relations corresponding to formulation 1 is presented in Fig. 5. The formulation in Eq. 6 shows that in order to obtain a thermal conductivity, there should  In the previous section we pointed out that bmR is composed of a kinetic (R k ) and a potential (R p ) portion and further stated that R p should not contribute to the thermal conductivity of the non-convective solid system. For a diagonal element, Einstein relation can be separated and written as We calculated the first and second terms on the right hand side of the Eq. 15 and plotted in Fig. 7. It can be seen that Einstein relation for R k is equal to the Einstein relation for total R. The second term on the right hand side of the Eq. 15 is just a flat line and do not add to the total. Accordingly, the coupled term in Eq. 15, on the average, does not contribute to the thermal conductivity.

N-Body Potential
We selected silicon system represented by the Tersoff potential for many body simulations.
We applied the same procedure to calculate Einstein relation for 10×10×10 supercell of Si, see Fig. 8. By also simulating 6×6×6 and 8×8×8 we have established that the calculated conductivity is well converged within 5%. The calculated thermal conductivity of silicon is 160.5±10.0 W/mK which is approximately 10% different than the experimental value of 142 -148 W/mK 38,39 at 300 K. Again, the Einstein relation for R k is almost identical to the one for total R. This confirms that the thermal conductivity does not depend on how the potential energy is split between interacting atoms and only the phonons are relevant to calculations. Ladd et al. 40 draw a similar conclusion that the heat current should be independent of the localization of the potential energy in a system at or near equilibrium.
In their study, Schelling et al. 28 re-derived the heat current for different ways of energy distribution for a 3-body potential and calculated the thermal conductivity of silicon by using HCACF. Their results also did not differ for the reasons stated here. However, splitting may be important for the system where convective motion is dominant.

B. Porosity and Ge-Clusters in Silicon
Heterogeneous systems contain a mixture of elements having different masses, elastic con-  The minimum thermal conductivity in bulk silicon was estimated as 1.0-1.3 W/mK [41][42][43] when the structure is amorphous. These studies suggest that for amorphous silicon the phonon mean free path is on the order of average nearest neighbor distance. Later experimental studies claimed that even lower thermal conductivities than the ones in amorphous systems are possible in doped crystalline silicon with nanopores 39,44,45 . This modification, to lower extend (i.e. an order of magnitude decrease from bulk thermal conductivity) was also observed in microporous silicon by Song and Chen 46 . This result is not surprising because in addition to decrease in phonon mean free paths and relaxation times due to internal surface scattering, there are less conduction channels due to missing silicon atoms. We have tested these ideas on porosity with the corrected R and Einstein relation. The porous silicon structures and the corresponding thermal conductivities are given in Fig. 10 and Table I respectively. In Table I, a* represents the same octahedron in Fig. 10 without germanium atoms inside. Porosity is defined as the ratio of number of missing atoms to number of atoms in a perfect system. We have also presented the ratio of 2 bonded atoms to 4 bonded atoms (r2) and the ratio of 3 bonded atoms to 4 bonded atoms (r3) as an indication of the internal surface area. The tubular pores in Fig. 10 In general, as the porosity and internal surface area increase, the thermal conductivity decreases as seen in Table I    It overcomes the difficulties inherent in taking the analytical derivative of R to obtain heat current. This also removes any ambiguity that arises from fitting the HCACF to a function in order to calculate thermal conductivity. Partition of R into potential and kinetic terms, moreover, demonstrates clearly the null contribution of potential energy term for any nondiffusive system. However, this is not case for systems with diffusion (directly contributing mainly as the convective heat flow). We also have demonstrated the utility of the method for bulk Si structures. With the studies performed on porous Si (cylindrical pores, and holes) and Ge-nanodots embedded into Si matrix, we have explored specific features of R-based or HCACF based thermal conductivity calculations. In these systems, we observe increased noise on HCACF when compared to pure Si simulations, in contrast R-based calculation for these complex systems produced smoother behavior leading to well defined thermal conductivity values. In our earlier studies, we have extensively taken advantage of this formulation in determining thermal conductivity of Si, Si-Ge nanostructures, 49