Statistical Analysis and Molecular Dynamics Simulations of the Thermal Conductivity of Lennard–Jones Solids Including Their Pressure and Temperature Dependencies

Aspects of the thermal conductivity, λ, of a Lennard–Jones (LJ) solid along an isotherm and the sublimation line are studied using equilibrium molecular dynamics (MD) simulations. A reformulation of the Green–Kubo time correlation function expression for λ in the form of a probability distribution function (PDF) of single trajectory contributions (STC) exhibits the same characteristic statistical trends as found previously for liquids, even at high pressures and low temperatures. The analysis reveals that for short periods of time the thermal conductivity can be negative. This feature is evident along the sublimation line isobar and a low‐temperature isotherm going to high densities. Along the isobar and isotherm lines, λ is to a good approximation a power law in temperature and density, respectively. This behavior is used in a more general thermodynamics‐based analysis description of the state point dependence of the thermal conductivity. The heat flux autocorrelation function increasingly develops a damped oscillatory appearance as pressure increases or temperature decreases, consistent with the phonon formulation of thermal conductivity.

aspects of the pressure and temperature dependence of the thermal conductivity of the LJ solid phase using MD, are compared with experimental thermal conductivity data of solid argon and previous MD simulations. This preliminary investigation builds on pioneering MD studies of LJ solid thermal conductivity as a function of pressure and density in refs. [6][7][8][9][10][11]. At the temperatures of concern here, energy waves associated with the lattice dynamics or phonons can be important in the thermal transport, which favors an atomistic model similar to MD rather than a continuum approach. [6,[8][9][10][11] The focus of interest here is, λðP, TÞ, and to express and interpret the P and T dependencies within the framework of classical thermodynamic relations developed in the literature many decades ago, see for example those given in refs. [12,13]. It is timely to take another look at these treatments as a number of the quantities which are difficult to obtain experimentally can now be computed directly by MD simulation. Therefore these equations can be more fully evaluated and hence exploited better than was possible decades ago.
In Section 2, the theory and LJ MD results relating to formal aspects of the PDF route to the thermal conductivity of solids are discussed. In Section 3, the MD study is extended to consider the effects of pressure at constant temperature. Section 4 covers the dependence of the thermal conductivity on temperature along the sublimation (P ≃ 0) line. Comparisons with the temperature dependence of the thermal conductivity of solid argon which can be represented by the LJ potential are made. A summary of the conclusions of this work is made in Section 5.

PDF Route to the Thermal Conductivity
In this section, some fundamental statistical mechanical aspects underpinning the thermal conductivity based on a PDF description are discussed. The new emphasis here is its application to the solid phase, whereas previous applications of this theory were confined to liquids. [4,5]

Theory and MD Results
Refs. [4,5] discuss the background to the STC methodology and here we will concentrate on its application to the thermal conductivity. The GK method for calculating the thermal conductivity of an equilibrium system involves the heat flux vector, J q , [2,14] where N is the number of molecules in volume V (i.e., the volume of the simulation cell here). The relative velocity between molecules i and j is v ij ¼ v i À v j where v i is the velocity of particle i, and r ij is the pair separation vector between the two molecules. The first derivative of the pair potential, ϕðrÞ, with respect to r is denoted by ϕ 0 . Also is the total energy of a molecule i in the fluid. The thermal conductivity, λ, by GK is, [14,15] where k B is Boltzmann's constant, and where "x" is a dummy variable in this and the subsequent equations. The angular brackets, h : : : i represent an average over time origins, and C T ðtÞ is the heat flux time autocorrelation function (TACF). The quantity, λ GK ðtÞ is known as the time-dependent thermal conductivity, whose limit at long times is the thermal conductivity. In practice, for any given simulation length, there is a limit on how large the upper limit of time (t) can be before statistical noise dominates and there is no benefit in carrying out the integration of C T to longer time. Despite many publications having been written on this problem, [16][17][18] it still is an unavoidable feature of the implementation of the GK method which limits the precision of the final value for the transport coefficient. The STC expression for the thermal conductivity, λ ST is The thermal conductivity STC is denoted by λ u . Each λ u may be considered to constitute a single collective dynamical "event" in the heat transmission process. In the last line of Equation (4), the quantity, P, is the PDF of λ u , which note is implicitly a function of t. The PDF of λ u for long times has a significant negative λ u region. The first moment of P also defined in the last line of Equation (4) is the thermal conductivity, which must be positive as there must be a net production of heat arising from the imposition of a temperature gradient in an experimental system or nonequilibrium molecular dynamics (NEMD). [4] The STC decomposition analysis throws light on the distribution of dynamical events that combine to yield the thermal conductivity. It also enables, in principle, PDF analytic tools to be exploited in the context of transport coefficients. Its computation adds little computational time to the simulation as the STC are already calculated as an intermediate step in the GK procedure.
Ref. [5] discusses the possible PDF definitions which are relevant to the STC quantity. The PDF, PðxÞ, in Equation (4) is for the case in which the quantity, x, is in real units (i.e., which here is LJ, principally), where À∞ < x < ∞, and PðxÞ is normalized to unity in this complete argument range. As PðxÞ is asymmetric on the negative and positive sides, it is also convenient to define the PDFs, P À ðxÞ for the negative argument range, À∞ < x < 0 and P þ ðxÞ for 0 < x < ∞. The PDF on each side is normalized separately so that the two integrals are unity. The standard deviation of the two PDFs, σ À and σ þ are defined through, σ 2 À ¼ ∫ 0 À∞ x 2 P À ðxÞ dx and σ 2 þ ¼ ∫ ∞ 0 x 2 P þ ðxÞ dx. It was also found, [5] that a further simplification can usefully be made, in which the argument of the PDF is normalized by the appropriate www.advancedsciencenews.com www.pss-b.com standard deviation. This new PDF definition, P s,AE , is defined through the following relationships It was found in ref. [5] that within statistics and for not too small a value of t that, P s,À ðÀx=σ À Þ ¼ P s,þ ðx=σ þ Þ ≡ P 0 ðyÞ for x > 0 and where y ¼ Àx=σ À or x=σ þ , as appropriate. The nondimensionalized PDFs, P s,À and P s,þ are therefore observed to be symmetric on the negative and positive argument sides.
The LJ pair potential, ϕ LJ ðrÞ ¼ 4ϵ½ðσ=rÞ 12 À ðσ=rÞ 6 , where r is the center-to-center separation of the two particles, was used in the MD simulations. The simulated quantities are mostly reported in units of ϵ and σ, and the mass of the molecule, m. The MD time step was, Δt ¼ 0.004= ffiffiffi ffi T p , and the interaction truncation distance, r c , was 3.5 (see ref. [3]). The number of particles in the simulation cell, N, was typically 1372. Some of the computed quantities reported here are expressed in real units using the usual LJ potential parameters for argon, with conversion factors between reduced and real units for the relevant quantities in this study given in Table 1. The parameters used for argon were, ϵ=k B ¼ 119.8 K and σ ¼ 0.3405 nm.
The computations were conducted typically for 6 Â 10 5 to 1 million time steps. Most of the simulations carried out were for the solid state at low temperatures below the triple-point temperature. For comparison, some results for fluid-state points are reported, notably at the reduced number density, ρ ¼ 0.8442 and T ¼ 0.722, which is close to the triple point. This state point was first introduced by Levesque, Verlet and coworkers, [19,20] and has been used as a standard or reference state point in many simulation studies since. Simulations were carried out mainly using constant temperature or NVT dynamics where N, V, and T are the number of particles, volume, and temperature, respectively. The Nosé-Hoover thermostat [21,22] with a time constant of 3 LJ time units was employed. Figure 1a shows the LJ phase diagram on a ρ À T plane. The solid region of interest is in the bottom right-hand side of the diagram. The isotherm and isobar (sublimation) lines investigated here are highlighted in the diagram. The solid region of the LJ phase diagram focused on in this study is shown in Figure 1b, plotted on a T À ρ plane, which is a typical form of presentation for such investigations. The melting line above the triple-point temperature ≃ 0.69 shown in the figure is taken from various MD and Monte Carlo simulation sources. [23][24][25][26][27] The experimental argon melting line data, [28] shown in the figure (see the figure caption for further details) agrees within their mutual statistical uncertainties with the LJ data, giving support for the adequacy of the LJ potential in representing argon, at least in the present context. It is worth noting that the sublimation line (the solid-vapor coexistence boundary below the triple-point temperature) is not an extrapolation of the melting line to lower temperatures (which is also shown in the figure) but takes a different course on the ρ, T surface. There is, in fact, an increase in the coexisting density of the solid phase with decreasing temperature below the triple-point temperature, [29,30] rather that a decrease as would be predicted by extrapolating the melting line data to 0 K. A recent accurate analytic expression for the sublimation line, [31] is also shown in the figure as a solid black line. This is compared with points used in these simulations for the sublimation line based on the P ¼ 0 approximate condition. The two sets of values agree very well. The sublimation line is to a very good aproximation a P ¼ 0 isobar. Figure 2 shows the PDF, Pðλ u Þ, for liquid and solid-state points. The STC quantity, λ u , is in LJ units. The PDFs on the negative and negative STC sides are shown. In both cases, the negative side is larger for jλ u j close to zero, but decays more rapidly than the positive-side PDF for larger argument values (i.e., there is a cross-over). Figure 3 presents these data where the λ u are normalized by the standard deviation, calculated separately for each side. Therefore, there are two standard deviations of relevance in the P s,AE description. The thermal conductivity PDF, P s,AE ðλ u =σ AE Þ, defined in Equation (5) for the liquid and solid-state points are shown in Figure 3. In Figure 3a,b, the arguments of the PDFs are λ u =σ AE and jλ u =σ AE j, respectively. The figure shows that the PDFs expressed in this way are statistically the same for the liquid and solid states. In ref. [5], it was shown that these reduced quantity PDFs can be represented well by a sum of exponentials (three were found sufficient), which for the thermal conductivity can be written using the following formula

MD Simulations
where a i and b i are constants fitted to the simulation data. The constants are, a 1 ¼ 0.92726, a 2 ¼ 1.45707, and a 3 ¼ 0.50139, and b 1 ¼ 2.57893, b 2 ¼ 27.86215, and b 3 ¼ 0.87189. Figure 4 shows the time-dependent thermal conductivity, λðtÞ, obtained by the GK formula of Equation (3) and the STC procedure of Equation (4), for the liquid and solid states in Figure 4a,b, respectively. The separate negative and positive-side contributions of the SCT method are presented. One notable feature of Figure 4a is that the positive and negative contributions to the thermal conductivity increase continuously with time without reaching a plateau, even though the combined value reaches a limiting plateau value at comparatively short times. The STC and GK results are statistically indistinguishable, as they should be, as they are both constructed from exactly the same numerically acquired quantities. The same trends were noted in ref. [5] for the shear viscosity. The qualitative difference in behavior of the component and total value of the time-dependent STC was attributed in ref. [5] to random fluctuations in the STC ("viscuit") which increase with t. The same trends are evident in λ u ðtÞ. Statistically, these fluctuations have equal probability of occurring for the negative and positive-side PDFs, and for  The magenta and green horizontal lines mark out the boundary between different regions. The sublimation (solid-vapor coexistence region) is below the green line. Frame (b), the solid part of the phase diagram of the LJ system and argon which is of relevance to this work, given in LJ reduced units. The melting line above the triple-point temperature (≃0.69) derived from molecular simulation. [23][24][25][26][27] using the LJ potential is shown (open circles). The lower temperature sublimation line ("SUBL.") was mapped out by MD simulations carried out here. The density for a given T is determined when P ¼ 0 in the MD simulation (shown as filled-in circles). Experimental data for argon, [28] (solid squares) where ρ ¼ 4=a 3 and the lattice parameter, aðTÞ ¼ A þ BT þ CT 2 , where in LJ units A, B and C ¼ 1.666109, À0.09506882, and 0.01187766. This line obtained in ref. [28] is also shown on the figure as a blue solid line with positive slope. From ref. [31], an analytic fit to the sublimation line is, ρ ¼ 1.09151 À 0.14081T À 0.04152T 2 þ 0.01828T 3 À 0.18547T 4 þ0.31686T 5 À 0.24139T 6 , which is also shown on the figure as a solid black line with negative slope.
www.advancedsciencenews.com www.pss-b.com practical purposes, their effects (almost) cancel out in the total transport coefficient. The data in Figure 4a is plotted as a function of ffi ffi t p . The negative and positive PDF side contributions to the thermal conductivity are linear with ffi ffi t p in the long time limit. This is the expected dependence for stochastic processes, [16][17][18] which therefore supports this hypothesis on the origin of the continual growth in the negative and positive STC contributions. The problem of establishing a long-time plateau value with increasing correlation time due to random fluctuations at long time is also manifested in the GK and its formally equivalent transformation, the Einstein-Kubo-Helfand (EKH) method. [32][33][34][35] Whichever formulation is used to compute the transport coefficient, one cannot escape the need to accommodate the long-time growth in statistical uncertainty. It should not be concluded, however, that the PDFs are entirely noise. At short times, of order, the correlation time of the dynamical processes leading to the transport coefficient, the negative and positive sides of the PDF, Pðλ u Þ have a different functional form and do not cancel out. This difference contains the information needed to determine the thermal conductivity. Figure 5 shows the PDF ratios, R ¼ P þ =P À plotted on a log scale for the solid, where the abscissa is jλ u ðtÞj. The lnðRÞ data are linear with jλ u j (for not too large values)   with the slope decreasing as t increases. The linearity of the data and dependence on t suggests the analytic form, P þ ðλ u ðtÞ=P À ðÀjλ u ðtÞjÞ ¼ BðtÞ expðAðtÞjλ u ðtÞtjÞ, where AðtÞ and BðtÞ are time-dependent constants. A simple approximate model derived in ref. [5] predicted that It was shown in Figure 11 of ref. [5] that σ þ ðtÞ and σ À ðtÞ converge with increasing t. The trends in R seen in Figure 5 were also noted in ref. [5] for the shear viscosity and thermal conductivity of the liquid phase.
To summarize this section, the analysis of the MD data for solid LJ systems reveals that qualitatively and for certain properties (e.g., the P s,AE ) quantitatively, the same behavior is exhibited in both the liquid and solid states for the thermal conductivity and its STC decomposition. The constancy of P s,AE across the fluid and solid-phase diagrams is a useful simplification, but this does not mean we can dispense with simulations to evaluate the thermal conductivity at any specific state point, as the STC standard deviations σ À and σ þ still need to be computed by MD at that state point to give the value of the thermal conductivity. Note also that the STC or "viscuit" is not a trivial quantity similar to the mean of a system property over a time t, but contains information about individual property correlations over time (the STC involves the product of two heat fluxes at different times).
The present investigation is extended to consider the density and pressure dependence of the thermal conductivity. The MD results and theoretical analysis are presented in the next section. Table 1 presents the conversion factors between LJ reduced and SI units for the properties computed in this study. Table 2 gives the thermal conductivity computed by MD for some representative regions of the fluid (including liquid) and solidphase diagrams. State point E which was studied in ref. [6] is far into the solid phase at a very high density (ρ ¼ 1.414). The value for λ obtained here is in agreement with the value obtained in that work, within the mutual statistical uncertainties. State point E simulations were conducted with N ¼ 256 particles, as in the previous work, and a further simulation was carried out with N ¼ 1372, which showed a significantly lower value of the thermal conductivity. Although it has been found that the MD λ have a very weak N-dependence, [8] down to a molar volume, V m , value of 22 μm 3 mol À1 , the E state point has a much smaller V m of 17 μm 3 mol À1 than those for which the previous N-dependent analysis was made. It is not surprising that for this state point there is a noticeable N-dependence in the value of the transport coefficient. As discussed earlier, the main difficulty can be in assigning the plateau value in the GK integrand which can be difficult to identify within a certain range of possible values. Table 3 lists thermodynamic quantities and λ as a function of density and pressure, given in LJ reduced units. The contribution of the repulsive and attractive parts of the LJ potential to the total energy per particle, u, from the simulation are also given in Table 3 for future reference. Figure 6 shows the radial distribution function, gðrÞ, of the LJ solid at three densities along the T ¼ 0.578 isotherm (70 K for argon). The pair separation is scaled by ρ 1=3 to account for the homogeneous part of the change in gðrÞ due to density, which enables density-induced structural changes to be more readily identified. The peaks become sharper with increasing density (pressure). At the highest density on the figure, ρ ¼ 1.5 additional peaks appear for %rρ 1=3 > 2.5, which are not present for lower densities (pressures). These new features indicate a  www.advancedsciencenews.com www.pss-b.com subtle structural change in the long range order of the lattice at high (%GPa) pressures, which may have an influence on the phonon spectrum. Figure 7a presents the pressure dependence of the density at constant temperature. Two literature formulas can be fitted to the MD data very well. These are the Tait, [36,37] Table 3. Calculated properties along the T ¼ 0.578 isotherm in the solid phase. The thermal conductivity was calculated using the GK method formula given in Equation (3). The quantities, u r , u a , and u are the r À12 , r À6 and total contributions to the average potential energy per particle, respectively. P is the total pressure. The numbers in brackets for λ are the uncertainties in the final digits. The quantities are given in LJ units.    www.advancedsciencenews.com www.pss-b.com

Pressure Dependence of the Thermal Conductivity
and Bair, [37] ρðPÞ ¼ ρð0Þ formulas. The parameters, ρð0Þ, A T and B T , and A B and B B are given in the figure caption. Figure 7b shows that these expressions can also be used to fit ln λðρÞ quite well. It is intuitively reasonable that the thermal conductivity will follow the pressure dependence of the isothermal compressibility, κ T ¼ ð ∂ ln ρ= ∂PÞ T , which is the slope of the curve in Figure 7a. It is generally the case that the density dependence of λ or any other transport coefficient is not easy to measure in experiment, in contrast to what can be achieved in MD simulation where density is the natural input parameter. The pressure is usually the independent variable in tribology and geology, and other fields where the effects of large loads are of importance. Figure 8a presents the density dependence of ln λ along the T ¼ 0.578 isotherm, and Figure 8b shows ln λ as a function of ln ρ. Both plots are statistically close to being linear in both ρ and ln λ, and the regression constants are given in the figure caption. There is no systematic difference between the thermal conductivity values obtained with N ¼ 864 and 1372 even at the highest densities considered. Interpretation of these trends is aided using the following thermodynamic relationship, [12] ∂ ln λ Experimental measurements of the parameter, g, defined in Equation (10), give values typically between 6-10 for crystals, and the value obtained by linear regression in Figure 8b is 8.1 AE 0.1, which is in the middle of this range. This indicates that for this temperature and wide density range, λ ≃ Aρ 8 , where A is a constant valid up to a pressure of ≃4 GPa for argon. This simple density dependence was also discovered by Tretiakov and Scandolo, [7] who carried out similar simulations using an exponential-6 (E6) potential for solids at higher temperatures, and obtained an exponent of 6. This is consistent with the present data, as the LJ potential is steeper at short range than E6. There is evidence that the E6 potential is a better representation of argon for pressures in excess of %1 GPa. [7,38,39] On the basis of the present results, combined with those of Tretiakov and Scandolo, it is clear that the simple density scaling of the thermal conductivity covers a wide range of the solid-phase diagram and is relatively independent of the potential form (for argon at least).
The pressure dependence of the thermal conductivity therefore from Equation (10) is seen to track well that of the isothermal compressibility. A comparison of Figure 8a,b indicates that ln λ is quite linear when plotted against both ρ and ln ρ, with the latter performing slightly better in this respect. An important conclusion is that within the simulation statistics, g is independent of ρ and therefore P, along the isotherm which could be exploited in further analytical analysis. In contrast, for hard spheres, g is a strong function of P in the solid phase. [40] Various attempts have been made to predict the value of the parameter, g, by approximate model treatments. For example, [12] a simple lattice approximation for the thermal conductivity is, λ ≃ ffiffi ffi 2 p νC v r À1 , where C v is the isochoric heat capacity per unit volume, ν, is a characteristic average lattice vibrational frequency, and r is the nearest-neighbor distance in the lattice. For an FCC (face-centered cubic) lattice, r ¼ 2 1=6 =ρ 1=3 , and as C v per molecule changes by only a few percent along the isotherm (at T ¼ 0.578, it only increases from 2.71 to 2.95 between the www.advancedsciencenews.com www.pss-b.com densities ρ ¼ 0.949 and 1.52), the following approximation can be derived using this approximate formula for the thermal conductivity is the Grüneisen parameter which is a measure of the anharmonic contribution to dynamical properties. The frequency, ν, can be approximated by where ω E is the Einstein frequency, which within a cell lattice or Lindemann approximation, [41][42][43][44] is given by where m is the mass of the particle, and N nn is the number of nearest neighbors for a molecule (i.e., 12 here for an FCC lattice), and ϕ 0 and ϕ 00 are the first and second derivatives of the potential with respect to r (here, the nearest neighbor distance). This approximation using Equation (12) and (13) leads to There is an even simpler formula in the literature which only retains the highest order derivative in the numerator and denominator of Equation (14), [45] i.e., γ ¼ Àrϕ 000 =6ϕ 00 . Equation (14) gives a value of ≃1 between ρ ¼ 0.9 À 1.5, and the simplified version of Ramani and Ghodgaonkar, [45] decays monotonically from 4.3 to 2.6 in the same density range, rather than the expected value ≃8 obtained directly from the MD data. Therefore, approximate lattice expressions for the Grüneisen parameter and hence g give values which are quite sensitive to the truncations and approximations used in their derivation. Despite this drawback, it is undeniable that from the MD results, ln λ and ln ρ are to a good approximation linearly related over a significant pressure range, which is equivalent to the thermal conductivity being proportional to a high power of ρ up to vary large pressures. Figure 9a presents the density dependence of the normalized TACF, defined in Equation (3) of five state points along the T ¼ 0.578 isotherm of the LJ solid. The figure shows that with increasing density this function progressively becomes more oscillatory, which indicates that a collective excitation ("phonon") mechanism of heat transfer becomes more important in that limit. Figure 9b shows the standard deviation normalized PDFs defined in Equation (5) for four of the state points represented on Figure 9a. The absolute value of λ u is used to show both sides of the PDF on the figure which facilitates their comparison. The figure demonstrates that within simulation statistics the negative and positive sides are symmetric and the same for these four state points. This is consistent with Figure 3 and the fluid-state behavior found in ref. [5]. Therefore, the PDFs in reduced form are independent of pressure and statistically the same as for a liquid composed of the same types of particle.
Temperature is the other intensive thermodynamic variable of much importance in experimental situations. The next section considers the temperature dependence of λ, along the sublimation or P ≃ 0 line.   figure. Frame (b) gives the corresponding P s,AE for STC averaging times: (i) t ¼ 2.63, (ii) t ¼ 10.5, (iii) t ¼ 9.9, and (iv) t ¼ 9.2 for (on the figure), ρ ¼ 0.999, 1.0125, 1.200, and 1.500, respectively. The absolute value of λ u is plotted, and both negative and positive sides of each PDF are shown and seen to be coincident within the simulation statistics.

Temperature Dependence of the Thermal Conductivity
In this section, the temperature dependence of the thermal conductivity along the sublimation line is considered. Again λ is calculated by MD using the GK and STC methods. The temperature range of the solid is below the triple-point temperature.
The triple-point parameters of the LJ system are, T ¼ 0.695ð1Þ, P ¼ 0.001ð1Þ, and for the liquid phase, ρ ¼ 0.845ð1Þ, and the density of the coexisting solid is 0.961 (1). [31,46,47] The vapor-solid or sublimation line has been the focus of a number of molecular simulation investigations. [31,38,[48][49][50] A number of experimental measurements of the temperature dependence of the thermal conductivity of solid argon along this line have been made. [51,52] The TðρÞ boundary has a discontinuity in slope at the triple-point temperature. [29,30] This line was computed here by carrying out simulations in which the density was adjusted iteratively to give an average total pressure of zero at a given temperature. A crystal cannot be thermodynamically stable at negative pressure and one can assume the pressure of the vapor in this temperature range is relatively small so that the P ¼ 0 isobar is a good approximation to the sublimation line. These densities and the corresponding thermal conductivities obtained by MD are listed in Table 4.
The thermal conductivity in crystalline van der Waals solids has a strong temperature dependence. At low temperatures below about 10 K, it increases as %T 3 due to the heat capacity dependence predicted within the Einstein and Debye models, an effect which cannot be captured by classical MD. The thermal conductivity peaks at about T ¼ 5 K. Above about this temperature, the thermal conductivity decreases with increasing T, and is affected by phonon-phonon scattering due to the effects of anharmonicity. The scattering mean free path has %1=T dependence. This effect is most pronounced in the vicinity of the low temperature peak in λðTÞ. The phonon-description of the thermal conductivity of noble gas solids based on a pair potential description has been demonstrated to give good agreement with experimental data. [53] Figure 10a shows the temperature dependence of the thermal conductivity along the sublimation line calculated in this study. The MD simulation data from different sources are also shown in the figure. They agree very well with those performed in this study. A transition between phonon-dominated and molecular diffusive mechanisms of energy transfer occurs as temperature decreases between the low and high-temperature limits, respectively. [9] Figure 10a confirms a scaling of the %1=T b form where the exponent b is %1.5 which is the same value as that obtained in ref. [9] (see the figure caption for further details), giving some support for the phonon-phonon scattering mechanism and the damped oscillatory behavior of C T ðtÞ evident in Figure 9a. The simulation data tends to be slightly lower than the experimental λ data of ref. [52]. This is a difficult quantity to measure experimentally and some systematic differences between different samples used in that study are evident in the figures of that work. This discrepancy could originate in the experimental values, but this would need to be verified by new measurements. Nevertheless, the three simulation studies whose λ values are given on Figure 10a are, within the statistical uncertainties, in agreement.
Following the macroscopic thermodynamic treatment given in ref. [12], assuming again that λ is a state function, using standard thermodynamic relations [12] ∂ ln λ where α ¼ Àð ∂ ln ρ= ∂TÞ P in Equation (15) is the isobaric expansivity (it is almost always a positive quantity). As discovered earlier, it is a good approximation to take g to be a constant, independent of ρ (at least along an isotherm). Figure 10b is a plot of ln ρ against T, which is seen to be nearly linear in T for small T, indicating in that temperature range, α ≃ A expðÀBTÞ where A and B are positive constants (given in the figure caption). The plot of ln ρ against ln T (not shown) is even more curved downward for large T than Figure 10b. It is tempting to take ð ∂ ln λ= ∂TÞ ρ to be a constant as has been done in the literature. [12] Although a reasonable assumption, as isochoric quantities tend to be less temperature dependent than isobaric ones, this proposal would require further practical investigation to confirm to what extent this is the case. Table 4. Calculated properties along the P ≃ 0 isobar or sublimation line in the solid phase (apart from the last three rows). The thermal conductivity was calculated using the GK method formula given in Equation (3). The numbers in brackets for λ are the uncertainties in assigning the thermal conductivity from the plateau value, in the last digits. The quantities are given in LJ units. Key: all data is for N ¼ 1372. www.advancedsciencenews.com www.pss-b.com Table 4 shows that for the temperature range covered and P ≃ 0, the density decreases with increasing temperature which indicates a positive isobaric expansivity over the entire temperature range studied. A positive α requires a positive Grüneisen parameter, γ, as α ¼ γρC v =κ T , [54] and the three remaining quantities are all positive. This is ensured in the approximate formula in Equation (14) because ϕ 000 << 0 at the nearest-neighbor distance. The third derivative term outweighs the contributions from the first and second derivative terms in that formula using the LJ potential. Figure 11a shows the isobaric dependence of ln λ against ln T, which complements the lin-lin plot of Figure 10a. The power law decay in T of λ is evident in Figure 11a.    Figure 10. Frame (a), temperature dependence of the thermal conductivity along the sublimation line. "EXPT." refers to experimental data from ref. [52]. Three sets of MD data are shown on the figure, where "TW" is MD data from this work in Table 4, "MK" is MD data from ref. [11] and "KLYK" is MD data from ref. [9]. The solid line is a fit to the MD data with the functional form, λ ≃ A=T b , where A ¼ 118ð14Þ and b ¼ 1.46ð6Þ in LJ reduced units, with the statistical uncertainties given in brackets. Frame (b), the temperature dependence of the P ¼ 0 sublimation isobar. The black line is a linear regression fit of the MD data up to a temperature of T ¼ 0.2 in LJ reduced units, where the intercept and slopes are 0.0884(3) and À0.140(3), respectively. The red line is an analytic fit to sublimation line data from ref. [31].
www.advancedsciencenews.com www.pss-b.com An alternative approach is to follow the spirit of Equation (10) and use ∂ ln λ As with g, the quantity h is not easy to access in experiment but can be evaluated in MD in a more routine way. Figure 11b shows the isobaric dependence of ln λ against ln ρ. The isobaric temperature dependence in Figure 11a is almost linear in a narrow ln ρ range 0.05-0.07, which reflects the expansivity behavior in Figure 10b. Clearly, the density dependence of λ along the P ¼ 0 isobar is more complicated than its temperature dependence. This dependence may ultimately provide a link between the isotherm and isobar trends in the thermal conductivity.
The temperature-dependent normalized TACF, C T ðtÞ, defined in Equation (3) for five temperatures along the P ¼ 0 isobar are shown in Figure 12a. The figure shows that with decreasing temperature C T ðtÞ develops a more pronounced damped oscillatory decay appearance. Just as for the analogous isotherm TACF in Figure 9b which shows the effects of increasing pressure, this suggests that heat transfer becomes more phonon-like in origin as the intensive variable becomes more "extreme" (here where there is a decrease in temperature). Figure 12b shows the standard deviation normalized PDFs defined in Equation (5) for the same four state points. The absolute value of λ u is used to show both sides of the PDF and the extent to which they collapse on a single curve. The figure demonstrates that the negative and positive sides are symmetric and independent of temperature for the higher three temperatures. The simulation carried out at T ¼ 10 K produced a PDF, which is symmetric but slightly systematically lower than the other three. This may be a consequence of the dominance of phonons at very low temperatures. The phonon collision mean free path length scales as %1=T, so the small departure from scaling at this temperature may be a finite size effect.

Conclusions
This work is an extension and a new application of the so-called "viscuit" or "single trajectory" (ST) decomposition of the GK formulae for obtaining transport coefficients. [4,5] For the last 50 years, the GK procedure has been used in MD computer simulations as a standard tool to determine the transport coefficients of liquids. The simulation is broken up into a sequence of trajectories each of which acts as the start of a ST contribution to the total transport coefficient value. These starting points in time are referred to as "time origins". The normal GK practice is to form an average correlation function from all of these individual components and then integrate this with time. The plateau value of the integral is proportional to the transport coefficient of interest. In the viscuit approach, each time origin is integrated, and then the statistics of these individual contributions to the transport coefficient (here, the thermal conductivity, λ) are analyzed to provide further insights into the distribution of dynamical "events" which lead to a specific value of λ. In principle, the statistical tools of PDF could be taken advantage of to provide further insights into the molecular origins of the processes that lead directly to the transport coefficient. This new method brings the evaluation of the transport coefficient within the framework of PDF theory. The ST approach makes use of information which has not been exploited in past applications of GK, even though the necessary ST contributions are already available as a step in the GK computational procedure.  The new aspect of the ST treatment is to explore whether the PDF scaling features for liquids discovered in ref. [5] also apply to crystalline solids. Simulations of the LJ FCC lattice below the triple-point temperature carried out in this study reveal that the universal features in the PDF found for liquids in ref. [5] also apply to solids. In this case, the thermal conductivity ST PDF expressed in terms of the standard deviations of the two sides treated separately is symmetric, and collapse onto those of the liquid. It is worth noting that many of the ST trajectories contribute a negative amount to the total value of λ, but the positive contributions outweigh those on the negative side of the PDF, leading to a total thermal conductivity that is positive. This means that the thermal conductivity can be expressed in terms of standard deviations of the ST on the two sides (further details are given in ref. [5]).
This work shows therefore from the thermal conductivity viscuit PDF of the LJ solid that for short periods of time the thermal conductivity of the solid is negative. Therefore, especially for nanoscale atomic solid systems, the thermal conductivity of the system can be negative, which is an effect that may be important to consider and could even be exploited in practical applications. This behavior is present over the range of temperatures and pressures considered in this study.
Another new feature of this work is an investigation of the pressure and temperature dependence of the thermal conductivity, λ, of the LJ solid. These were computed principally along an isotherm, and an isobar close to the sublimation line on the phase diagram, respectively. It is shown that along the isotherm, λ, is to a good approximation proportional to a high power of the density. The pressure dependence of λ follows a similar analytic form to that of the pressure dependence of the density. The temperature dependence of λ along the P ¼ 0 isobar is a power law in T, confirming previous MD investigations. The TACFs increasingly developed a damped oscillatory appearance as pressure increases or temperature decreases, suggesting an increasing contribution of phonon heat waves to the heat transmission mechanism in both of these cases.
In order for the GK integral to reach a plateau requires the heat flux time correlation functions to be computed for comparatively long times (for typically %50-100 reduced LJ time units) when the pressure is high and/or the temperature is low. The trends in the pressure and temperature dependence of the thermal conductivity are interpreted and formulated using thermodynamic relationships developed in pioneering work on thermal conductivity of liquids and solids carried out in the 1960s and 1970s.
The GK and viscuit adaptation explored in this study are based on Fourier's linear relationship between the heat flux and an applied temperature gradient. [55] Shear flow can also induce an anisotropic heat flux, which has recently been explored using nonequilibrium MD simulation. [56] There are NEMD methods using homogeneous synthetic equations of motion which enable the thermal conductivity in the zero temperature gradient limit to be computed. [57] The reformulated GK method used in this study is much more useful in the present context, however, as this method lends itself to showing up the negative contributions to the thermal conductivity which occur naturally as the system evolves through time and space. Future work could investigate the extent to which the intermolecular potential affects the value of the g-parameter (which characterizes the density dependence of the transport coefficient), in particular its density (pressure) dependence. Also further details on the similarities and differences between the pressure and temperature dependences of the thermal conductivity could be obtained. Is there an equivalence or mapping between the P and T dependencies of the thermal conductivity? Such interchangeability is known in other fields, such as the time-temperature superposition found in the polymer rheology literature. [58] This would provide greater understanding and potential applications of the simulation results in practice. In regard to the isobaric trends, are there differences in qualitative behavior of the thermal conductivity as a function of temperature for different pressure values? Similarly, for isothermal state points to what extent do the trends change as a function of temperature?
With a wider perspective, recent research on high-speed dislocation movement and plastic relaxation processes in crystals, [59] demonstrates the importance of collective motion underpinning many physical effects. Coupling of these processes with atomistic treatments of heat waves in crystals would appear to be a natural extension of this work (i.e., heat conduction in crystals with defects). Also, similar MD studies of heat conduction in glassy systems using the same analytical tools would be another possible direction of this work.