# Elementary vibrational model for thermal conductivity of Lennard-Jones
fluids: Applicability domain and accuracy level

S. A. Khrapak<EMAIL_ADDRESS>Joint Institute for High Temperatures,
Russian Academy of Sciences, 125412 Moscow, Russia  A. G. Khrapak  Joint
Institute for High Temperatures, Russian Academy of Sciences, 125412 Moscow,
Russia

###### Abstract

Exact mechanisms of thermal conductivity in liquids are not well understood,
despite rich research history. A vibrational model of energy transfer in dense
simple liquids with soft pairwise interactions seems adequate to partially
fill this gap. The purpose of the present paper is to define its applicability
domain and to demonstrate how well it works within the identified
applicability domain in the important case of the Lennard-Jones model system.
The existing results from molecular dynamics simulations are used for this
purpose. Additionally, we show that a freezing density scaling approach
represents a very powerful tool to estimate the thermal conductivity
coefficient across essentially the entire gas-liquid region of the phase
diagram, including metastable regions. A simple practical expression serving
this purpose is proposed.

##  I Introduction

The thermal conductivity coefficient is an important characteristic of a
material. It can vary greatly depending on specific substance as well on its
phase state. While heat transport in gases, plasmas and solids is relatively
well understood, mechanisms of heat transfer in liquids have remained elusive.

A vibrational model of thermal conductivity in simple liquids with soft
pairwise interactions has been recently proposed and discussed Khrapak
(2021a). In this model it is assumed that atoms in liquids exhibit solid-like
oscillations about temporary equilibrium positions corresponding to a local
minimum on the system’s potential energy surface Frenkel (1955); Stillinger
and Weber (1982); Zwanzig (1983). The equilibrium positions are not fixed like
in solids, but are allowed to diffuse (and this is why liquids can flow).
However, this diffusion occurs on long enough time scales, which are
irrelevant for the process of energy transfer. Furthermore, liquid is
approximated by a quasi-layered structure with layers perpendicular to the
temperature gradient and separated by the distance $\Delta=\rho^{-1/3}$, where
$\rho$ is the atomic number density. The average interatomic separation in
each quasi-layer is also $\Delta$. Then elementary consideration based on the
assumption that the energy difference between the atoms in neighbouring layers
is transferred at an average frequency of their solid-like vibrations leads
immediately to a simple expression Khrapak (2021a)

$\lambda\simeq c_{\rm v}\frac{\langle\omega\rangle}{2\pi\Delta}.$ (1)

Here $\lambda$ is the thermal conductivity coefficient (which does not include
the Boltzmann constant and hence is measured in cm-1s-1), $c_{\rm v}$ is the
specific heat at constant volume, and $\langle\omega\rangle$ is the average
vibrational frequency that approximates the energy transfer rate. Actually,
specific heat at constant volume $c_{\rm v}$ appears in Eq. (1) under the
assumption of a small thermal expansion coefficient. Dense liquids with soft
interactions near the liquid-solid phase transition are essentially
incompressible and this represents a good approximation.

We note in passing that the vibrational paradigm sketched above is useful not
only in relation to the problem of thermal conductivity of dense fluids. It
also allows to relate the self-diffusion and the shear viscosity coefficients
in the form of Stokes-Einstein relation without the hydrodynamic diameter
Zwanzig (1983); Costigliola  _et al._ (2019); Khrapak and Khrapak (2021a);
Khrapak (2021b). The present paper is, however, mostly focused on the thermal
conductivity mechanism.

Equation (1) emphasises the relation between the thermal conductivity and
collective modes properties of dense liquids. To evaluate the average
vibrational frequency one would require to know the liquid vibrational density
of states (VDOS). This is a formidable task, because collective properties can
greatly differ from one liquid system to another and depend considerably on
the state point in the phase diagram. Some progress has recently been reported
Zaccone and Baggioli (2021); Stamper  _et al._ (2022). Zaccone and Baggioli
have developed an analytical model for VDOS, based on overdamped Langevin
liquid dynamics Zaccone and Baggioli (2021). Distinct from the Debye
approximation, $g(\omega)\propto\omega^{2}$, for solids, the universal law for
liquids reveals a linear relationship, $g(\omega)\propto\omega$, in the low-
energy region. Stamper et al. have confirmed this universal law with
experimental VDOS measured by inelastic neutron scattering on real liquid
systems Stamper  _et al._ (2022). Nevertheless, the applicability regime and
accuracy level of this model still require clarification. In the meantime, it
is natural to employ some approximations and simplifications which can differ
depending on the type of liquid under investigation (some examples are
provided below).

Interestingly enough, Eq. (1) can be reduced to the known previous results
under special simplifying assumptions about the vibrational properties Khrapak
(2021a). In the simplest approximation all atoms are oscillating with the same
Einstein frequency $\Omega_{\rm E}$ (this approximation is known as the
Einstein model in the solid state physics). Averaging is then trivial,
$\langle\omega\rangle=\Omega_{\rm E}$. We get

$\lambda\simeq c_{\rm v}\frac{\Omega_{\rm E}}{2\pi\Delta}.$ (2)

This is very similar to to the results obtained by Rao Rao (1941) and later by
Horrocks and McLaughlin Horrocks and McLaughlin (1960, 1963). This form is
particularly suitable for extremely soft interactions such as Coulomb and
screened Coulomb potentials, relevant in the plasma-related context Khrapak
(2021c, d, 2023a). In dense liquids near freezing conditions we can set
$c_{\rm v}\sim 3$, according to Dulong-Petit law. The Einstein frequency at
the liquid-solid phase transition can be estimated using the Lindeman’s
melting criterion in its simplest formulation. Then, quasi-universality of the
thermal conductivity coefficient at freezing conditions emerges in the form
Rao (1941); da C. Andrade (1952)

$\lambda\simeq{\rm const}\sqrt{\frac{T}{m\Delta^{4}}},$ (3)

where $T$ is the temperature in energy units ($=k_{\rm B}T$) and $m$ is the
atomic mass. The value of the constant is about $\simeq 10$ for usual simple
model systems such as Lennard-Jones, Coulomb, screened Coulomb, as well as
monatomic liquids (e.g. liquefied noble gases) Khrapak and Khrapak (2022a, b).
It is somewhat larger for the hard-sphere fluid at freezing and reaches values
$\simeq 15-20$ for liquids with more complex molecular structure Khrapak
(2023b).

Alternatively, we can relate the average frequency to the sound velocity.
Namely, consider an acoustic dispersion of the form $\omega=kc_{\rm s}$, where
$k$ is the wave-number and $c_{\rm s}$ is the sound velocity. The
characteristic wave-number for energy transfer between nearest neighbours is
$k\sim 2\pi/\Delta$. This yields $\langle\omega\rangle\sim 2\pi c_{\rm
s}/\Delta$ and hence

$\lambda\sim c_{\rm v}\frac{c_{\rm s}}{\Delta^{2}}.$ (4)

This resembles Bridgman’s expression for the thermal conductivity coefficient
Bridgman (1923). Actually, Bridgman postulated that the energy is transferred
at the speed of sound, thus leading to a linear correlation between the
thermal conductivity coefficient and the sound velocity. He also used a
constant coefficient of two, that is $\lambda\simeq 2c_{\rm s}/\Delta^{2}$.
This numerical coefficient remained somewhat controversial, values between 2
and 3 were used in the literature Zhao  _et al._ (2021); Xi  _et al._ (2020);
Khrapak (2021a); Bird  _et al._ (2002). In a recent extensive study of
correlations between the thermal conductivity coefficient and the sound
velocity of various liquids Khrapak (2023b); Khrapak and Khrapak (2023), it
has been demonstrated that linear correlations are well reproduced for model
liquids as well as real monatomic and diatomic liquids. However, they are less
convincing in polyatomic molecular liquids. The actual coefficient of
proportionality in Bridgman’s formula is not fixed. It is about unity for
monatomic liquids and generally increases with molecular complexity.

Just like in solids, dense liquids support one longitudinal and two transverse
collective modes Hansen and McDonald (2006); Balucani and Zoppi (1994); Ohta
and Hamaguchi (2000); Hosokawa  _et al._ (2015); Bryk  _et al._ (2017). A
Debye-like approximation for the vibrational spectra of dense liquid is thus
not completely irrelevant and has demonstrated reasonable success in deriving
the Stokes-Einstein coefficient without the hydrodynamic diameter Zwanzig
(1983) and estimating thermodynamic properties of dense liquids Bolmatov  _et
al._ (2012); Khrapak and Yurchenko (2021). In this approach the dispersion
relations of one longitudinal and two transverse collective modes are
approximated by their corresponding acoustic asymptotes $\omega_{\rm
l}(k)=kc_{\rm l}$ and $\omega_{\rm t}(k)=kc_{\rm t}$, terminating at
respective cut-off wavelengths $k_{\rm cut}$. Application to the problem of
thermal conductivity results in the expression Khrapak (2021a)

$\lambda\simeq\frac{1}{4}\left(\frac{3}{4\pi}\right)^{1/3}c_{\rm
v}\frac{c_{l}+2c_{t}}{\Delta^{2}}.$ (5)

Using $c_{\rm v}\simeq 3$ near freezing we get a formula similar to that of
minimal thermal conductivity model proposed by Cahill and Pohl Cahill and Pohl
(1989); Cahill  _et al._ (1992) for amorphous solids.

Equation (5) can be applicable to liquids with conventional acoustic
dispersion relation of the longitudinal collective mode, but should not be
used when dispersion is non-acoustic (e.g. in the case of Coulomb one-
component plasma and weakly screened Yukawa systems). In such cases actual
dispersion relations can be used to perform averaging Khrapak (2021c, b)).
Another complication arises due to the presence of the so called $k$-gap
(zero-frequency domain at sufficiently long wavelengths) in the dispersion
relation of the liquid transverse mode, which has received considerable
attention in recent years Goree  _et al._ (2012); Yang  _et al._ (2017);
Trachenko and Brazhkin (2015); Khrapak  _et al._ (2019); Kryuchkov  _et al._
(2019); Kryuchkov and Yurchenko (2021). The effect of $k$-gap in the
transverse mode has not yet been discussed in detail in the context of the
vibrational paradigm of transport properties in liquids and will not be
addressed here.

Formula (5) provides a good compromise between simplicity and accuracy.
Notably, it does not contain free parameters and adjustable coefficients.

One of us demonstrated previously that Eq. (5) applies reasonably well to
dense Lennard-Jones (LJ) fluids Khrapak (2021a). However, only a single
exemplary slightly supercritical isotherm was considered in that work. The
purpose of this paper is to provide a detailed and extensive verification of
its applicability domain and accuracy level across the LJ fluid phase diagram.
We demonstrate that for the LJ fluid, Eq. (5) allows to express the thermal
conductivity coefficient in terms of thermodynamic properties such as specific
heat, excess internal energy, and excess pressure. An appropriate equation of
state is then used to evaluate the thermal conductivity coefficient and to
quantify the accuracy of the approach. This is a remarkable example of a
direct relation between transport and thermodynamics. The applicability domain
is defined with reference to a gas-to-liquid dynamical crossover in
supercritical fluids (Frenkel line in the phase diagram), which represents an
important current research topic. Finally, we reiterate our recent freezing
density scaling approach to transport properties Khrapak and Khrapak (2021b,
2022c, 2022b), and put forward an ad hoc expression for the thermal
conductivity coefficient. Overall, the results reported represent an important
step toward better understanding main mechanisms and peculiarities of
transport phenomena in the liquid state.

##  II Methods

The LJ model is one of the most popular and best studied systems in condensed
matter research. It combines relative simplicity with adequate representation
(at least at the qualitative level) of interatomic interactions in real
substances, exhibiting steep short-range repulsion and softer long-range
attraction. The LJ potential is

$\phi(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(6)

where $\epsilon$ and $\sigma$ are the energy and length scales (or LJ units).
The reduced density and temperature expressed in LJ units are
$\rho^{*}=\rho\sigma^{3}$, $T^{*}=T/\epsilon$.

Transport properties of the LJ system have been extensively studied in the
literature. Recent overview of available simulation data can be found in Refs.
Bell  _et al._ (2019); Harris (2020); Allers  _et al._ (2020). Particularly
extensive and useful datasets have been published by Meier et al. Meier
(2002); Meier  _et al._ (2004a, b) and by Baidakov et al. Baidakov  _et al._
(2011, 2012); Baidakov and Protsenko (2014). Transport data have been
tabulated along different isotherms in a wide regions of the LJ system phase
diagram. Good agreement between the two datasets for overlapping regimes has
been reported Harris (2020). Previously, we used the numerical results of
Meier Meier (2002) to demonstrate the adequacy of Eq. (5). Using the data for
the thermal conductivity coefficient, specific heat, reduced energy and
pressure along a single close-critical isotherm $T^{*}=1.35$ tabulated in Ref.
Meier  _et al._ (2004a) we documented a very good accuracy of Eq. (5) in the
dense liquid regime. Here we first make use of the thermal conductivity data
from Ref. Baidakov and Protsenko (2014), which cover a rather extended area on
the LJ system phase diagram.

Figure 1: (Color online) Phase diagram of the LJ system on the density-
temperature plane ($\rho^{*}$, $T^{*}$). The squares correspond to the fluid-
solid coexistence boundaries as tabulated in Refs. Sousa  _et al._ (2012), the
fluid-solid coexistence curves are simple fits by Eq. (7). The liquid-vapour
coexistence boundary is plotted using the formulas provided in Ref. Heyes  _et
al._ (2019). The reduced triple point and critical temperatures are
$T^{*}_{\rm tr}\simeq 0.69$ Sousa  _et al._ (2012) and $T^{*}_{\rm c}\simeq
1.33$ Heyes  _et al._ (2019), respectively. The crosses correspond to the
state points investigated in Ref. Baidakov and Protsenko (2014). Only data
points with the temperature above the triple point are considered. The dashed
line corresponding to a fixed reduced density $\rho/\rho_{\rm fr}\simeq 0.6$
marks the onset of the vibrational model applicability (for further details
see the text).

The phase diagram of LJ system is shown in Fig. 1. Here the solid-liquid
coexistence data are taken from Ref. Sousa  _et al._ (2012). For convenience,
they are fitted by simple expressions

$\displaystyle T^{*}_{\rm fr}$ $\displaystyle\simeq
2.111(\rho^{*})^{4}-0.615(\rho^{*})^{2},$ (7) $\displaystyle T^{*}_{\rm m}$
$\displaystyle\simeq 1.988(\rho^{*})^{4}-1.019(\rho^{*})^{2},$

where subscripts “${\rm fr}$” and “${\rm m}$” denote freezing and melting.
respectively. The functional form of these expressions arises in various
simple approaches to fluid-solid coexistence in the LJ model Rosenfeld (1976a,
b); Khrapak  _et al._ (2010); Khrapak and Morfill (2011); Pedersen  _et al._
(2016); Costigliola  _et al._ (2016). The liquid-vapour coexistence boundary
is plotted using the formulas provided in Ref. Heyes  _et al._ (2019). The
crosses are the state points for which the thermal conductivity coefficients
was numerically evaluated by means of equilibrium molecular dynamics with the
use of the Green-Kubo formalism in Ref. Baidakov and Protsenko (2014). Some of
investigated state points correspond to metastable regions: superheated and
supercooled liquids, supersaturated vapour. Here we only consider state points
with the temperature above the triple point temperature (lower temperatures
are not relevant because they are beyond the applicability limits of the
vibrational model). An additional dashed line marks the onset of validity of
the vibrational model. This was identified in Ref. Khrapak and Khrapak (2021a)
from the analysis of the Stokes-Einstein (SE) product $D\eta(\Delta/T)$, where
$D$ is the self-diffusion coefficient and $\eta$ is the shear viscosity
coefficient. At low densities the SE product scales as $(\rho^{*})^{-4/3}$ as
it should in the gaseous regime. However, it approaches a constant asymptotic
value of $\simeq 0.15$ at sufficiently high densities, not too far from that
at freezing. The two asymptotes intersect at $\rho/\rho_{\rm fr}\simeq 0.35$
and this can be considered as an indication of the gas-to-fluid dynamical
crossover Khrapak (2022b). At about twice this density, $\rho/\rho_{\rm
fr}\gtrsim 0.6$, the SE product becomes practically constant and this can be
identified as a lower boundary of the validity of the SE relation and hence of
the vibrational picture of atomic dynamics Khrapak and Khrapak (2021a). Note,
that the analysis of various model systems (such as Lennard-Jones, Coulomb,
Yukawa) indicated that the onset of validity of vibrational dynamics
corresponds to nearly the same value of excess entropy, $s_{\rm ex}\simeq-2$
Khrapak and Khrapak (2021a), and this can be considered as a more general
condition of the validity of the vibrational model.

The numerical results for the thermal conductivity coefficient are
traditionally expressed in LJ units as

$\lambda^{*}=\frac{\lambda\sigma^{2}}{k_{\rm
B}}\left(\frac{m}{\epsilon}\right)^{1/2}.$ (8)

Since we measure temperature in energy units, the Boltzmann constant $k_{\rm
B}$ should disappear. Also, we find it more convenient to work with
macroscopically reduced units (sometimes referred to as Rosenfeld’s
normalization Rosenfeld (1999)). This is particularly advantageous when
comparing transport properties of different systems Rosenfeld (1977, 1999);
Khrapak (2023, submitted). The macroscopically reduced thermal conductivity
coefficient reads

$\lambda_{\rm R}=\lambda\rho^{-2/3}\left(\frac{m}{T}\right)^{1/2}.$ (9)

From Eqs. (8) and (9) a trivial relation emerges

$\lambda_{\rm R}=\frac{\lambda^{*}}{(T^{*})^{1/2}(\rho^{*})^{2/3}}.$ (10)

This relation has been used to evaluate $\lambda_{\rm R}$ from the available
numerical results.

To compare with the vibrational mechanism of heat transfer, the longitudinal
and transverse instantaneous sound velocities have to be evaluated. For the LJ
system these can be expressed using the excess energy and pressure of the
system as demonstrated by Zwanzig and Mountain Zwanzig and Mountain (1965)
(see also Ref. Khrapak (2020) for further details):

$\frac{c_{l}^{2}}{v_{\rm T}^{2}}=3-\frac{72}{5}u_{\rm ex}+11p_{\rm ex},$ (11)
$\frac{c_{t}^{2}}{v_{\rm T}^{2}}=1-\frac{24}{5}u_{\rm ex}+3p_{\rm ex}.$ (12)

Here $u_{\rm ex}=U/NT-3/2$ is the excess internal energy, $p_{\rm ex}=P/\rho
T-1$ is the excess pressure, $U$, $P$, and $N$ being the internal energy,
pressure, and the number of atoms, respectively. The sound velocities are
expressed in units of the thermal velocity $v_{\rm T}=\sqrt{T/m}$. It is a
simple exercise to verify that Cauchy relation is satisfied in this approach:

$\left(\frac{c_{l}}{v_{\rm T}}\right)^{2}-3\left(\frac{c_{t}}{v_{\rm
T}}\right)^{2}=2p_{\rm ex}.$ (13)

In order to evaluate the sound velocities, an appropriate equation of state
(EoS) is required. A large number of different EOS of have been proposed in
the literature to describe thermodynamic properties of LJ fluids (see e.g.
Refs. Stephan  _et al._ (2019, 2020) for review). Here we use the equation of
state developed by Thol et al. Thol  _et al._ (2016). This is an empirical
equation of state, formulated in terms of the Helmholtz free energy and based
on a large molecular simulation data set and thermal virial coefficients. Its
applicability range $0.7<T^{*}<9$, $\rho^{*}<1.08$ is sufficient for the
present purpose. This EoS demonstrates reasonable level of accuracy when
compared to others Stephan  _et al._ (2020), is relatively simple and
convenient in implementation. We use Thol et al. EoS to calculate the specific
heat $c_{\rm v}$, excess pressure $p_{\rm ex}$ and energy $u_{\rm ex}$. From
the latter thermodynamic quantities the sound velocities $c_{l}$ and $c_{t}$
are calculated. This is all input needed to calculate the thermal conductivity
coefficient from Eq. (5). Thus, application of the vibrational model to LJ
fluids represents a remarkable example when a transport property (thermal
conductivity in the considered case) is determined explicitly by thermodynamic
quantities and knowledge of EoS is sufficient to calculate it. In Section III
we will see how well theory compares with results from numerical simulations.

Figure 2: (Color online) Reduced thermal conductivity coefficient
$\lambda_{\rm R}$ of the Lennard-Jones fluid as a function of the reduced
density $\rho^{*}$. Results for eight isotherms are shown: $T^{*}=0.7$ (a),
0.85 (b), 1.0(c), 1.15 (d), 1.2 (e), 1.3 (f), 1.5 (g), and 2.0 (h). The
circles denote the numerical results from Ref. Baidakov and Protsenko (2014),
the solid curves are calculated by means of vibrational model using Eq. (5),
as discussed in the text. The shaded regions in each figure mark the
conditions of applicability of the theory.

##  III Results

Figure 2 demonstrates the comparison between the vibrational model, as
calculated using Eq. (5) complemented by the thermodynamic quantities from
Thol et al. EoS Thol  _et al._ (2016), with the MD results from Baidakov et
al. Baidakov and Protsenko (2014). The circles correspond to the results from
MD simulation. We have selected only the isotherms with $T^{*}>T^{*}_{\rm
tp}$, where vibrational model makes sense and Thol et al. EoS is applicable.
We have also considered supercritical densities for which the vibrational
model might be adequate. For this reason we, have omitted the two isotherms
$T^{*}=1.25$ and $T^{*}=1.35$ for which only irrelevant low-density data are
available. The remaining data for eight isotherms are depicted in Fig. 2
(a)-(h). The solid curves in each figure correspond to our theoretical
calculation.

The shaded regions in each figure correspond to the regime of applicability of
the present model. From the side of high densities the applicability is
limited by the fluid boundary of the fluid-solid coexistence. From the side of
low densities the onset of the applicability is either the condition
$\rho/\rho_{\rm fr}\simeq 0.6$ where vibrational picture is adequate or the
liquid boundary of the gas-liquid coexistence. Although metastable state
points are present in Fig. 2 they should not be used to judge the
applicability of the vibrational model. Neither the model itself was designed
to deal with coexisting phases, nor the EoS from Ref. Thol  _et al._ (2016)
can be expected to be reliable there.

Figure 3: (Color online) Relative disagreement between theory [Eq. (5)] and MD
simulations of Ref. Baidakov and Protsenko (2014). The quantity $(\lambda_{\rm
Theory}-\lambda_{\rm MD})/\lambda_{\rm Theory}$ is plotted as a function of
the reduced density $\rho^{*}$ for eight isotherms. Only the state points
within the applicability regime of the theoretical model are considered.

Careful examination of eight panels in Fig. 2 demonstrates a very good
agreement between the theoretical calculation and MD simulation in the regime
where theory is applicable. This is further quantified in Fig. 3. Relative
deviations are typically limited by few percent, increasing to $\simeq 5\%$
towards the boundaries of the applicability domain. Some simulation data for
the thermal conductivity reported previously have considerably larger standard
deviations compared to the observed deviations (see e.g. Table II in Ref.
Nasrabad  _et al._ (2006)). For this reason the agreement can be considered as
excellent, especially taking into account that the theoretical model does not
contain any free adjustable parameters.

To reinforce this conclusion and demonstrate that the documented agreement is
not the result of some fortunate coincidence, we have performed an additional
comparison using the set of simulation data reported by Galliero and Boned
Galliero and Boned (2009). The purpose is two-fold. First, in contrast to Ref.
Baidakov and Protsenko (2014), Galliero and Boned employed non-equilibrium
molecular dynamics simulations to obtain thermal conductivity of LJ fluids. In
this way, we can potentially quantify discrepancies between different methods,
if such exist. Second, the data set from Ref. Galliero and Boned (2009) covers
considerably higher range of temperatures and hence we can verify whether the
vibrational model remains meaningful as the temperature further increases.
Figure 4 presents the comparison between the theory and non-equilibrium
simulations. We observe that for $\rho^{*}\gtrsim 0.6$ the deviation between
theory and numerical experiment is mostly within numerical data uncertainty.
For a common temperature $T^{*}=2$ the theoretical model compares equally well
with equilibrium and non-equilibrium simulations, demonstrating no major
discrepancies between the two latter. For higher temperatures there is no any
sign that the theory may behave inappropriately.

Figure 4: (Color online) Reduced thermal conductivity coefficient
$\lambda_{\rm R}$ of the Lennard-Jones fluid as a function of the reduced
density $\rho^{*}$. The symbols denote the numerical results from the non-
equilibrium simulations in Ref. Galliero and Boned (2009). The solid curves
are calculated using Eq. (5).

Thus, where applicable the vibrational model of heat transport agrees very
well with the results from different numerical simulations. Yet, there is
another practical approach – the freezing density scaling (FDS) – which allows
to predict the transport properties, including the thermal conductivity
coefficient, in an even wider parameter regime. We would like to use this
opportunity and reiterate the application of FDS scaling to the thermal
conductivity coefficient of the LJ fluid. Section IV serves this purpose.

##  IV Freezing density scaling

It has been recently demonstrated that macroscopically reduced self-diffusion,
shear viscosity, and thermal conductivity coefficients of LJ fluids along
isotherms exhibit quasi-universal scaling on the density divided by its value
at the freezing point, $\rho_{\rm fr}$ Khrapak and Khrapak (2021b). Originally
considered as a useful empirical observation Khrapak and Khrapak (2021b), it
has been later discussed in the context of quasi-universal excess entropy
scaling and isomorph theory Khrapak and Khrapak (2022c, b); Heyes  _et al._
(2023). FDS approach implies that the transport coefficients are functions of
a single variable $\rho/\rho_{\rm fr}$, and thus it represents a very
convenient corresponding states principle to estimate transport properties in
LJ fluids. Importantly, it has been shown that FDS holds even at quite low
densities where neither the original form of the excess entropy scaling nor
the isomorph theory are expected to work. Additionally, the functional form of
the FDS scaling is similar (although not identical) to that in the hard-sphere
fluid. Thus, it can be expected that FDS is not a special property of LJ
fluids, but applies (possibly with some modifications) to a wider class of
fluids.

Figure 5: (Color online) Reduced thermal conductivity coefficient
$\lambda_{\rm R}$ vs reduced density $\rho/\rho_{\rm fr}$ demonstrating the
success of the freezing density scaling of transport coefficients. Symbols
correspond to the numerical data from Baidakov et al Baidakov and Protsenko
(2014). along different isotherms (see the label). The solid curve corresponds
to a quasi-universal fit of Eq. (14).

Detailed analysis of the numerical dataset from Ref. Baidakov and Protsenko
(2014) has been further extended to reinforce the FDS concept. The results are
shown in Fig. 5, demonstrating that the data tend to collapse on a universal
master curve. The only region where the FDS fails is the vicinity of the
critical point, where critical enhancement of the thermal conductivity
coefficient becomes important. For the dataset produced in Ref. Baidakov and
Protsenko (2014) this concerns only a few data points, mostly on the near-
critical isotherm $T^{*}=1.3$ (see Fig. 5). Importantly, FDS approach
seemingly works well for metastable state points, which are also shown in Fig.
5. No clear distinction between thermodynamically stable and metastable states
can be seen. Also, it is remarkable that the FDS approach applies to rather
low densities. This is rather unexpected, because at such low densities it is
natural to assume that any ”memory” about the location of the fluid-solid
phase transition and the freezing density is completely lost. However, it is
necessary to remark that in the low-density gaseous regime established methods
to calculate the transport cross sections and the transport properties of the
LJ system do exist (see e.g. Refs. Hirschfelder  _et al._ (1954, 1948); Smith
and Munn (1964); Khrapak (2014a, b); Kim and Monroe (2014); Kristiansen (2020)
and references therein for some related works).

The appealing quasi-universality of the data shown in Fig. 5 calls for an
appropriate mathematical description. Among various fitting functions that we
attempted, a particularly simple but appropriate form is provided by

$\lambda_{\rm
R}=\frac{\alpha}{{\mathcal{R}}^{2/3}}+\beta+\gamma{\mathcal{R}}^{\delta}.$
(14)

Here ${\mathcal{R}}$ is the density ratio, ${\mathcal{R}}=\rho/\rho_{\rm fr}$.
The first term is chosen as to provide a correct asymptote of the thermal
conductivity coefficient in the low-density limit, $\lambda_{\rm
R}\propto(\rho^{*})^{-2/3}$ Khrapak and Khrapak (2021b). Furthermore,
$\alpha$, $\beta$, $\gamma$, and $\delta$ are the fitting coefficients. Based
on the dataset for the supercritical isotherm $T^{*}=2$, the coefficients
$\alpha\simeq 0.43$, $\beta\simeq 1.18$, $\gamma\simeq 8.39$, and
$\delta\simeq 2.30$ have been obtained. The solid curve in Fig. 5 represents
this fit. It is representative for other isotherms, except those very near the
critical temperature, because the critical enhancement is not accounted for.
Importantly, it appears also representative for all metastable state points
that were investigated in Ref. Baidakov and Protsenko (2014).

##  V Conclusion

Thermal conductivity is an important characteristic of a material. While this
property is relatively well understood in gases, solids and plasmas, the same
cannot be said about liquids. In this paper we have considered the Lennard-
Jones fluid as an important simple model system to investigate the
applicability limits and check the accuracy of theoretical approaches to the
thermal conductivity coefficient. This program has been realised with the help
of an extensive MD simulation datasets provided by Baidakov et al. Baidakov
and Protsenko (2014) and Galliero and Boned Galliero and Boned (2009).

Two theoretical methods have been highlighted in this paper. The first one –
the vibrational paradigm of heat transfer – leads to an expression [Eq. (5)],
which relates the thermal conductivity coefficient to the thermodynamic
properties of the LJ fluid. Although very simplistic and to some extent naive,
the vibrational model is able to describe the numerical data on the thermal
conductivity coefficient with a remarkable accuracy. No adjustable parameters
are involved and the theoretical calculation requires only the input from
thermodynamics. The applicability region spans from the gas-liquid coexistence
boundary to the liquid-solid coexistence boundary for subcritical temperatures
and from a constant density ratio $\rho/\rho_{\rm fr}\simeq 0.6$ to the fluid-
solid coexistence boundary for supercritical temperatures.

The second method is the corresponding state principle, based on the freezing
density scaling approach. This is more intuitive approach, although solid
relations to the excess entropy scaling and the isomorph theory do exist. The
FDS approach implies that the reduced thermal conductivity coefficients
depends quasi-universally on the density reduced by its value at the freezing
point (at a given temperature). The applicability domain of the FDS approach
covers almost the entire gas-liquid-fluid region in the LJ system phase
diagram (probably excluding very low densities and deep gas-liquid coexistence
region). It is applicable even to metastable regions of superheated and
supercooled liquids and supersaturated vapour. This can make the FDS approach
a very useful practical tool to estimate the thermal conductivity coefficient
under various conditions and a tentative fitting formula is provided for this
purpose. The only domain where FDS clearly fails is the close vicinity of the
critical point, because it does not account for critical enhancement.

Overall, the reported results shed new light on the properties and mechanisms
of heat transport in liquids and provide new accurate methods for the
estimation of the thermal conductivity coefficient. This can be of interest
for researchers in condensed matter, physics of fluids, materials science and
beyond.

The authors have no conflicts of interest to disclose.

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * Khrapak (2021a) S. A. Khrapak, “Vibrational model of thermal conduction for fluids with soft interactions,” Phys. Rev. E 103, 013207 (2021a).
  * Frenkel (1955) Y. Frenkel,  _Kinetic theory of liquids_ (Dover, New York, NY, 1955).
  * Stillinger and Weber (1982) F. H. Stillinger and T. A. Weber, “Hidden structure in liquids,” Phys. Rev. A 25, 978–989 (1982).
  * Zwanzig (1983) R. Zwanzig, “On the relation between self-diffusion and viscosity of liquids,” J. Chem. Phys. 79, 4507–4508 (1983).
  * Khrapak (2022a) S. Khrapak, “Vibrational model of heat conduction in a fluid of hard spheres,” Appl. Sci. 12, 7939 (2022a).
  * Costigliola  _et al._ (2019) L. Costigliola, D. M. Heyes, T. B. Schrøder, and J. C. Dyre, “Revisiting the Stokes-Einstein relation without a hydrodynamic diameter,” J. Chem. Phys. 150, 021101 (2019).
  * Khrapak and Khrapak (2021a) S. A. Khrapak and A. G. Khrapak, “Excess entropy and Stokes-Einstein relation in simple fluids,” Phys. Rev. E 104, 044110 (2021a).
  * Khrapak (2021b) S. A. Khrapak, “Self-diffusion in simple liquids as a random walk process,” Molecules 26, 7499 (2021b).
  * Zaccone and Baggioli (2021) A. Zaccone and M. Baggioli, “Universal law for the vibrational density of states of liquids,” Proc. Natl. Acad. Sci. 118, e2022303118 (2021).
  * Stamper  _et al._ (2022) C. Stamper, D. Cortie, Z. Yue, X. Wang, and D. Yu, “Experimental confirmation of the universal law for the vibrational density of states of liquids,” J. Phys. Chem. Lett. 13, 3105–3111 (2022).
  * Rao (1941) M. R. Rao, “Thermal conductivity of liquids,” Phys. Rev. 59, 212–212 (1941).
  * Horrocks and McLaughlin (1960) J. K. Horrocks and E. McLaughlin, “Thermal conductivity of simple molecules in the condensed state,” Trans. Faraday Soc. 56, 206 (1960).
  * Horrocks and McLaughlin (1963) J. Horrocks and E. McLaughlin, “Temperature dependence of the thermal conductivity of liquids,” Trans. Faraday Soc. 59, 1709 (1963).
  * Khrapak (2021c) S. A. Khrapak, “Thermal conduction in two-dimensional complex plasma layers,” Phys. Plasmas 28, 010704 (2021c).
  * Khrapak (2021d) S. A. Khrapak, “Thermal conductivity of strongly coupled Yukawa fluids,” Phys. Plasmas 28, 084501 (2021d).
  * Khrapak (2023a) S. A. Khrapak, “Vibrational model of heat transfer in strongly coupled Yukawa fluids (dusty plasma liquids),” Plasma Phys. Rep. 49, 15–22 (2023a).
  * da C. Andrade (1952) E. N. da C. Andrade, “Viscosity and thermal conductivity of liquid argon,” Nature 170, 794–794 (1952).
  * Khrapak and Khrapak (2022a) S. A. Khrapak and A. G. Khrapak, “Minima of shear viscosity and thermal conductivity coefficients of classical fluids,” Phys. Fluids 34, 027102 (2022a).
  * Khrapak and Khrapak (2022b) S. A. Khrapak and A. G. Khrapak, “Freezing density scaling of fluid transport properties: Application to liquefied noble gases,” J. Chem. Phys. 157, 014501 (2022b).
  * Khrapak (2023b) S.A. Khrapak, “Bridgman formula for the thermal conductivity of atomic and molecular liquids,” J. Mol. Liq. 381, 121786 (2023b).
  * Bridgman (1923) P. W. Bridgman, “The thermal conductivity of liquids under pressure,” PNAAS 59, 141 (1923).
  * Zhao  _et al._ (2021) A. Z. Zhao, M. C. Wingert, R. Chen, and J. E. Garay, “Phonon gas model for thermal conductivity of dense, strongly interacting liquids,” J. Appl. Phys. 129, 235101 (2021).
  * Xi  _et al._ (2020) Q. Xi, J. Zhong, J. He, X. Xu, T. Nakayama, Y. Wang, J. Liu, J. Zhou, and B. Li, “A ubiquitous thermal conductivity formula for liquids, polymer glass, and amorphous solids,” Chin. Phys. Lett. 37, 104401 (2020).
  * Bird  _et al._ (2002) R. B. Bird, E. N. Lightfoot, and W. E. Stewart,  _Transport Phenomena -_ (J. Wiley, New York, 2002).
  * Khrapak and Khrapak (2023) S. A. Khrapak and A. G. Khrapak, “Sound velocities in liquids near freezing: Dependence on the interaction potential and correlations with thermal conductivity,” Phys. Fluids 35, 077129 (2023).
  * Hansen and McDonald (2006) J.-P. Hansen and I. R. McDonald,  _Theory of Simple Liquids -_ (Elsevier, Amsterdam, 2006).
  * Balucani and Zoppi (1994) U. Balucani and M. Zoppi,  _Dynamics of the Liquid State_ (Clarendon Press, Oxford, 1994).
  * Ohta and Hamaguchi (2000) H. Ohta and S. Hamaguchi, “Wave dispersion relations in Yukawa fluids,” Phys. Rev. Lett. 84, 6026–6029 (2000).
  * Hosokawa  _et al._ (2015) S. Hosokawa, M. Inui, Y. Kajihara, S. Tsutsui, and A. Q. R. Baron, “Transverse excitations in liquid Fe, Cu and Zn,” J. Phys.: Condens. Matter 27, 194104 (2015).
  * Bryk  _et al._ (2017) T. Bryk, A. Huerta, V. Hordiichuk, and A. D. Trokhymchuk, “Non-hydrodynamic transverse collective excitations in hard-sphere fluids,” J. Chem. Phys. 147, 064509 (2017).
  * Bolmatov  _et al._ (2012) D. Bolmatov, V. V. Brazhkin, and K. Trachenko, “The phonon theory of liquid thermodynamics,” Sci. Rep. 2, 421 (2012).
  * Khrapak and Yurchenko (2021) S. A. Khrapak and S. O. Yurchenko, “Entropy of simple fluids with repulsive interactions near freezing,” J. Chem. Phys. 155, 134501 (2021).
  * Cahill and Pohl (1989) D. G. Cahill and R.O. Pohl, “Heat flow and lattice vibrations in glasses,” Solid State Commun. 70, 927–930 (1989).
  * Cahill  _et al._ (1992) D. G. Cahill, S. K. Watson, and R. O. Pohl, “Lower limit to the thermal conductivity of disordered crystals,” Phys. Rev. B 46, 6131–6140 (1992).
  * Goree  _et al._ (2012) J. Goree, Z. Donkó, and P. Hartmann, “Cutoff wave number for shear waves and Maxwell relaxation time in Yukawa liquids,” Phys. Rev. E 85, 066401 (2012).
  * Yang  _et al._ (2017) C. Yang, M. T. Dove, V. V. Brazhkin, and K. Trachenko, “Emergence and evolution of the k-gap in spectra of liquid and supercritical states,” Phys. Rev. Lett. 118, 215502 (2017).
  * Trachenko and Brazhkin (2015) K. Trachenko and V. V. Brazhkin, “Collective modes and thermodynamics of the liquid state,” Rep. Progr. Phys. 79, 016502 (2015).
  * Khrapak  _et al._ (2019) S. A. Khrapak, A. G. Khrapak, N. P. Kryuchkov, and S. O. Yurchenko, “Onset of transverse (shear) waves in strongly-coupled Yukawa fluids,” J. Chem. Phys. 150, 104503 (2019).
  * Kryuchkov  _et al._ (2019) N. P. Kryuchkov, L. A. Mistryukova, V. V. Brazhkin, and S. O. Yurchenko, “Excitation spectra in fluids: How to analyze them properly,” Sci. Rep. 9, 10483 (2019).
  * Kryuchkov and Yurchenko (2021) N. P. Kryuchkov and S. O. Yurchenko, “Collective excitations in active fluids: Microflows and breakdown in spectral equipartition of kinetic energy,” J. Chem. Phys. 155, 024902 (2021).
  * Khrapak and Khrapak (2021b) S. A. Khrapak and A. G. Khrapak, “Transport properties of Lennard-Jones fluids: Freezing density scaling along isotherms,” Phys. Rev. E 103, 042122 (2021b).
  * Khrapak and Khrapak (2022c) S. A. Khrapak and A. G. Khrapak, “Freezing temperature and density scaling of transport coefficients,” J. Phys. Chem. Lett. 13, 2674–2678 (2022c).
  * Bell  _et al._ (2019) I. H. Bell, R. Messerly, M. Thol, L. Costigliola, and J. C. Dyre, “Modified entropy scaling of the transport properties of the Lennard-Jones fluid,” J. Phys. Chem. B 123, 6345–6363 (2019).
  * Harris (2020) K. R. Harris, “Thermodynamic or density scaling of the thermal conductivity of liquids,” J. Chem. Phys. 153, 104504 (2020).
  * Allers  _et al._ (2020) J. P. Allers, J. A. Harvey, F. H. Garzon, and T. M. Alam, “Machine learning prediction of self-diffusion in Lennard-Jones fluids,” J. Chem. Phys. 153, 034102 (2020).
  * Meier (2002) K. Meier,  _Computer Simulation and Interpretation of the Transport Coefficients of the Lennard-Jones Model Fluid (PhD Thesis)_ (Shaker, Aachen, 2002).
  * Meier  _et al._ (2004a) K. Meier, A. Laesecke, and S. Kabelac, “Transport coefficients of the Lennard-Jones model fluid. I. Viscosity,” J. Chem. Phys. 121, 3671–3687 (2004a).
  * Meier  _et al._ (2004b) K. Meier, A. Laesecke, and S. Kabelac, “Transport coefficients of the Lennard-Jones model fluid. II. Self-diffusion,” J. Chem. Phys. 121, 9526–9535 (2004b).
  * Baidakov  _et al._ (2011) V.G. Baidakov, S.P. Protsenko, and Z.R. Kozlova, “The self-diffusion coefficient in stable and metastable states of the Lennard–Jones fluid,” Fluid Phase Equilibria 305, 106–113 (2011).
  * Baidakov  _et al._ (2012) V. G. Baidakov, S. P. Protsenko, and Z. R. Kozlova, “Metastable Lennard-Jones fluids. I. Shear viscosity,” J. Chem. Phys. 137, 164507 (2012).
  * Baidakov and Protsenko (2014) V. G. Baidakov and S. P. Protsenko, “Metastable Lennard-Jones fluids. II. Thermal conductivity,” J. Chem. Phys. 140, 214506 (2014).
  * Sousa  _et al._ (2012) J. M. G. Sousa, A. L. Ferreira, and M. A. Barroso, “Determination of the solid-fluid coexistence of the n - 6 Lennard-Jones system from free energy calculations,” J. Chem. Phys. 136, 174502 (2012).
  * Heyes  _et al._ (2019) D. M. Heyes, D. Dini, L. Costigliola, and J. C. Dyre, “Transport coefficients of the Lennard-Jones fluid close to the freezing line,” J. Chem. Phys. 151, 204502 (2019).
  * Rosenfeld (1976a) Y. Rosenfeld, “Additivity of melting curves,” Chem. Phys. Lett. 38, 591–593 (1976a).
  * Rosenfeld (1976b) Y. Rosenfeld, “Universality of melting and freezing indicators and additivity of melting curves,” Mol. Phys. 32, 963–977 (1976b).
  * Khrapak  _et al._ (2010) S. A. Khrapak, M. Chaudhuri, and G. E. Morfill, “Liquid-solid phase transition in the Lennard-Jones system,” Phys. Rev. B 82, 052101 (2010).
  * Khrapak and Morfill (2011) S. A. Khrapak and G. E. Morfill, “Accurate freezing and melting equations for the Lennard-Jones system,” J. Chem. Phys. 134, 094108 (2011).
  * Pedersen  _et al._ (2016) U. R. Pedersen, L. Costigliola, N. P. Bailey, T. B. Schrøder, and J. C. Dyre, “Thermodynamics of freezing and melting,” Nature Commun. 7, 12386 (2016).
  * Costigliola  _et al._ (2016) L. Costigliola, T. B. Schrøder, and J. C. Dyre, “Freezing and melting line invariants of the Lennard-Jones system,” Phys. Chem. Chem. Phys. 18, 14678–14690 (2016).
  * Khrapak (2022b) S. A. Khrapak, “Gas-liquid crossover in the Lennard-Jones system,” J. Chem. Phys. 156, 116101 (2022b).
  * Rosenfeld (1999) Y. Rosenfeld, “A quasi-universal scaling law for atomic transport in simple fluids,” J. Phys.: Condens. Matter 11, 5415–5427 (1999).
  * Rosenfeld (1977) Y. Rosenfeld, “Relation between the transport coefficients and the internal entropy of simple systems,” Phys. Rev. A 15, 2545–2549 (1977).
  * Khrapak (2023, submitted) S. A. Khrapak, “Elementary vibrational model for transport properties of dense fluids,” Phys. Rep. (2023, submitted).
  * Zwanzig and Mountain (1965) R. Zwanzig and R. D. Mountain, “High-frequency elastic moduli of simple fluids,” J. Chem. Phys. 43, 4464–4471 (1965).
  * Khrapak (2020) S. A. Khrapak, “Sound velocities of Lennard-Jones systems near the liquid-solid phase transition,” Molecules 25, 3498 (2020).
  * Stephan  _et al._ (2019) S. Stephan, M. Thol, J. Vrabec, and H. Hasse, “Thermophysical properties of the Lennard-Jones fluid: Database and data assessment,” J. Chem. Informat. Model. 59, 4248–4265 (2019).
  * Stephan  _et al._ (2020) S. Stephan, J. Staubach, and H. Hasse, “Review and comparison of equations of state for the Lennard-Jones fluid,” Fluid Phase Equilibria 523, 112772 (2020).
  * Thol  _et al._ (2016) M. Thol, G. Rutkai, A. Köster, R. Lustig, R. Span, and J. Vrabec, “Equation of state for the Lennard-Jones fluid,” J. Phys. Chem. Ref. Data 45, 023101 (2016).
  * Nasrabad  _et al._ (2006) A. E. Nasrabad, R. Laghaei, and B. C. Eu, “Molecular theory of thermal conductivity of the Lennard-Jones fluid,” J. Chem. Phys. 124, 084506 (2006).
  * Galliero and Boned (2009) G. Galliero and C. Boned, “Thermal conductivity of the Lennard-Jones chain fluid model,” Phys. Rev. E 80, 061202 (2009).
  * Heyes  _et al._ (2023) D. M. Heyes, D. Dini, S. Pieprzyk, and A. C. Brańka, “Departures from perfect isomorph behavior in Lennard-Jones fluids and solids,” J. Chem. Phys. 158, 134502 (2023).
  * Hirschfelder  _et al._ (1954) J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird,  _The Molecular Theory of Gases and Liquids -_ (Wiley, New York, 1954).
  * Hirschfelder  _et al._ (1948) J. O. Hirschfelder, R. B. Bird, and E. L. Spotz, “The transport properties for non-polar gases,” J. Chem. Phys. 16, 968–981 (1948).
  * Smith and Munn (1964) F. J. Smith and R. J. Munn, “Automatic calculation of the transport collision integrals with tables for the Morse potential,” J. Chem. Phys. 41, 3560–3568 (1964).
  * Khrapak (2014a) S. A. Khrapak, “Classical scattering in strongly attractive potentials,” Phys. Rev. E 89, 032145 (2014a).
  * Khrapak (2014b) S. A. Khrapak, “Accurate transport cross sections for the Lennard-Jones potential,” Eur. Phys. J. D 68, 276 (2014b).
  * Kim and Monroe (2014) S. U. Kim and C. W. Monroe, “High-accuracy calculations of sixteen collision integrals for Lennard-Jones (12-6) gases and their interpolation to parameterize neon, argon, and krypton,” J. Comput. Phys. 273, 358–373 (2014).
  * Kristiansen (2020) K. R. Kristiansen, “Transport properties of the simple Lennard-Jones/spline fluid I: Binary scattering and high-accuracy low-density transport coefficients,” Frontiers Phys. 8, 271 (2020).