Also at ]Center for Applied Physics and Technology, HEDPS, Peking University,
Beijing, China.

# On the computation of thermo-relaxing multi-component flows with the Baer-
Nunziato model

Chao Zhang<EMAIL_ADDRESS>Lifeng Wang  [<EMAIL_ADDRESS>Institute
of Applied Physics and Computational Mathematics, Beijing, China.

###### Abstract

In inertial confinement fusion (ICF) implosions, mixing the ablator into the
fuel and the hot spot is one of the most adverse factors that lead to ignition
degradation. Recent experiments in the Marble campaign at the Omega laser
facility and the National Ignition Facility (NIF) demonstrate the significance
of the temperature separation in heterogeneous mixing flows [Haines B.M., et
al. Nature Communications, 2020]. In the present work we provide an approach
to deal with thermally disequilibrium multi-component flows with the ultimate
aim to investigate the temperature separation impact on mixing and fusion
burn. The present work is two-fold: (a) we derive a model governing the multi-
component flows in thermal disequilibrium with transport terms, (b) we use the
derived model to study the Rayleigh–Taylor (RT) instability in thermally
relaxing multi-component systems. The model is reduced from the full
disequilibrium multi-phase Baer-Nunziato model in the limit of small Knudsen
number $Kn<<1$. Velocity disequilibrium is closed with the diffusion laws and
only one mass-weighted velocity is retained formally. Thus, the complex wave
structure of the original Baer-Nunziato model is simplified to a large extent
and the obtained model is much more computationally affordable. Moreover, the
capability to deal with finite temperature relaxation is kept. Efficient
numerical methods for solving the proposed model are also presented. Equipped
with the proposed model and numerical methods, we further investigate the
impact of thermal relaxation on the RT instability development at the inertial
confinement fusion (ICF) deceleration stage. On the basis of numerical
simulations, we have found that for the RT instability at an interface between
the high-density low-temperature component and the low-density high-
temperature component, the thermal relaxation significantly suppresses the
development of the instability.

††preprint: APS/123-QED

##  I Introduction

Mixing of the ablator into the fuel and the hot spot is considered to be one
of the most adverse factors that lead to ignition degradation in inertial
confinement fusion (ICF) implosions. Mixing takes place at different scales
driven by different mechanisms. At the macroscopic scale where $Kn<<1$, the
hydrodynamic instabilities (such as the Rayleigh–Taylor instability (RTI),
Richtmyer–Meshkov instability (RMI) and Kelvin–Helmholtz instability (KHI))
play predominate roles in causing mixing. With the continuous development of
the hydrodynamic instabilities, the flows transition into turbulence where
mixing happens at very different scales. At the small scale where
$Kn=\mathcal{O}(1)$ the mass diffusion is proceeding all the time as result of
random molecular motions (i.e., the kinetic effect). These two mechanisms have
fundamental difference and scaling laws. For example, the mixing length caused
by RTI and further induced turbulence can be scaled as $L_{mix}=\alpha
Agt^{2}$, where $\alpha$ is a constant, $A$ is the Atwood number, $g$ is the
acceleration, and $t$ the time. In contrast, the mixing length caused by mass
diffusion is expressed as $L_{mix}=\beta\sqrt{Dt}$, where $D$ is the mass
diffusivity and $\beta$ is a constant. Recent works demonstrate that the the
mass diffusion maybe the leading mixing mechanism for the implosion
experiments under moderate temperature and convergence[1].

The ICF mixing is usually categorized into two types, i.e., the atomic mixing
and the chunk mixing. In the former the mixing of the components takes place
on the atomic scale, while in the latter the constituents are separated from
each other as either at a rippled interface or in discrete clumps [2]. The ICF
mixing is intrinsically heterogeneous as a combination of atomic mixing and
non-atomic mixing. In the course of hydrodynamic development towards
turbulence, the chunk mixing dominate the early stage and the atomic mixing
takes over at later times. Meanwhile the mass diffusion is always producing
atomic mixing. From the perspective of numerical modeling, to discriminate
different types of mixing with direct numerical simulation is a formidable
task due to their very different characteristic scales. Therefore, one has to
rely on mixing models that allow the coarse-grained description of the atomic
mixing on an affordable grid (Figure 1). The grid does not resolve the small
scale mixing topology and the component mass fraction span over several
computational cells in the macroscopic description. In the mixing cells some
closure relation is needed to make the governing equations solvable. One of
the most frequently used closure is the temperature equilibrium, which is
assumed in the ICF-relevant mixing models such as the k-L model [3, 4] and the
BHR model [5, 6]. These models have been successfully used to simulate some
experiments on the Omega laser facility and the National Ignition Facility
(NIF). However, they give numerical results that deviate from experimental
measurements in the case of temperature separation. The root of such failure
is attributed to the thermal equilibrium assumption, which has been proved in
subsequent experiments at the Marble platform at NIF [7, 8] (referred to as
“Marble” henceforth).

Figure 1: The typical mass fraction distribution on the Eulerian grid in the
case of heterogeneous mixing.

As noted by Haines et al. [9, 7], strong ion temperature separation arises
during the shock flash. Mechanisms that drive the component temperatures
towards equilibrium include electron thermal conduction and local collision in
atomic mixing. They happen in a timescale comparable to that of the ICF
implosion (at least 1.6ns for the Marble campaign on the Omega facility [9]).
Thus, the ion temperatures of each species are not fully equilibrated in the
course of ICF implosions. This means that in numerical simulation the
components should have their own temperature in a mixing cell.

We aim to develop a multicomponent model to describe the temperature
separation in mixing. Note that the phenomenlohgical two-field turbulence
models [10, 11, 12, 13] has the potential to deal with turbulence mixing with
finite temperature relaxation. However, here we adopt a more direct and strict
derivation from a fully disequilibrium model - the BN model[14]. Due to the
complexity and relaxation stiffness of the original BN model, a hierarchy of
reduced models has been established (Figure 2). Reduced models are derived in
the limit of instantaneous relaxation of corresponding variable (chemical
potential, pressure, velocity, and temperature). For example, in the case
where the phase velocities relax instantaneously, one can derive the $u$-eq
model via the asymptotic analysis. Starting from the $u$-eq model, one can
further derive the $up$-eq model (i.e., the Kapila’s five equation model [15])
on the basis of the instantaneous phase pressure relaxation. Following such
procedure, a full hierarchy of reduced models can be obtained. A similar
hierarchy has been described by Lund [16], however, the velocity
disequilibrium that is vital for modeling mass diffusion is neglected there.

The applicability of a particular reduced model is determined by the
corresponding assumption on relaxation rates. In most applications the thermal
relaxation time is smaller than the mechanical (pressure and velocity)
relaxation times, thus models in the right subsidiary are scarcely used. In
term of the current work, we need a model that retains the thermal
disequilibrium. Moreover, since mass diffusion is related to the species
velocity difference, the velocity disequilibrium should also be maintained.
Thus, the possible models to satisfy these requirements are the original BN
model and the $p$-eq model. Here we focus on reformulation the former by
invoking the diffusion laws to close the velocity difference and abandoning
the terms of order $\mathcal{O}(\text{Kn}^{2})$. Formally the derived model
has one velocity, i.e., the mass weighted velocity. In fact, hyperbolic sub-
system of the model coincides with the six equation model presented in [17],
which is more robust than the Kapila’s five equation model in numerical
implementation. This reduction significantly simplifies the wave structure of
the hyperbolic subsystem of the model, thus improving the computational
efficiency.

Figure 2: The hierarchy of the reduced models of the Baer-Nunziato model.

We propose second-order methods for solving the hyperbolic-parabolic-
relaxation system. The numerical methods are validated against some benchmark
problems. Being equipped with the model and numerical methods, we then
investigate the RTI in the thermally relaxing multi-component flows,
especially the dependence of the mixing length on the thermal relaxation rate
and the initial temperature separation.

The rest of the article is organized as follows. In Section II we derive the
reduced temperature-disequilibrium model with diffusions. In Section III the
numerical methods for solving the proposed model are briefly described. In
Section IV some numerical results for validating the numerical methods are
presented. Analysis of the RTI problem is performed in particular detail to
investigate the impact of the thermal relaxation.

##  II Model formulation

###  II.1 The BN-type seven-equation model

The starting point of the following model formulation is the complete BN-type
seven-equation model [14, 18, 19, 20]. It reads:

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}}}{\partial{t}}+\nabla\bm{\cdot}(\alpha_{k}\rho_{k}\bm{u}_{k})=0,$
(1a)
$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}\bm{u}_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}\bm{u}_{k}\bm{\otimes}\bm{u}_{k}-\alpha_{k}\overline{\overline{T}}_{k}\right)=$
(1b)
$\displaystyle-\overline{\overline{T}}_{I}\bm{\cdot}\nabla{\alpha_{k}}+\mathcal{M}_{k},$
$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}E_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}E_{k}\bm{u}_{k}-\alpha_{k}\overline{\overline{T}}_{k}\bm{\cdot}\bm{u}_{k}\right)=$
(1c)
$\displaystyle-\bm{u}_{I}\bm{\cdot}\left(\overline{\overline{T}}_{I}\bm{\cdot}\nabla{\alpha_{k}}\right)+\bm{u}_{I}\mathcal{M}_{k}-p_{I}\mathcal{F}_{k}+\mathcal{Q}_{k}+q_{k}+\mathcal{I}_{k},$
$\displaystyle\frac{\partial{\alpha_{k}}}{\partial{t}}+\bm{u}_{I}\bm{\cdot}\nabla\alpha_{k}=\mathcal{F}_{k},$
(1d)

where the notations used are standard:
$\alpha_{k},\;\rho_{k},\;\bm{u}_{k},\;p_{k},\;\overline{\overline{T}}_{k},\;E_{k}$
are the volume fraction, phase density, velocity, pressure, stress tensor, and
total energy of phase $k$.

The BN-type formulation can also be derived by averaging the one-phase
governing equations[21, 22, 23]. Within such framework the volume fraction
$\alpha_{k}$ is the spatially average characteristic function $X_{k}$, which
takes the value of 1(or 0) in the presence (or absence) of the $k$-th
component, i.e.,

$\displaystyle\begin{split}X_{k}(\bm{r},t)=\left\\{\begin{array}[]{ll}1,&\text{the\;point\;}\bm{r}\text{\;in\;the\;material\;$k$},\\\
0,&\text{otherwise.}\end{array}\right.\end{split}$

Then, we have

$\alpha_{k}=\frac{1}{V}\int_{V}X_{k}\text{d}V=\frac{V_{k}}{V}.$

The volume fraction $\alpha_{k}$ can be perceived as a macroscopic measure of
spatial mixing extent when the microscopic interfaces can not be accurately
resolved. Moreover, the BN model can be used for describing both miscible and
immiscible flows according to different interpretations of $\alpha_{k}$ and
other variables [24].

For the sake of clarity we constrict our discussions within the scope of two-
phase flows, $k=1,2$. The phase density $\rho_{k}$ is defined as the mass per
unit volume occupied by $k$-th phase. The mixture density $\rho$ is the sum of
the partial densities $\alpha_{k}\rho_{k}$, i.e.,
$\rho=\sum{\alpha_{k}\rho_{k}}$. The last equation eq. 1d is written for only
one component thanks to the saturation constraint for volume fractions
$\sum_{k=1}^{2}\alpha_{k}=1$. The total energy is
$E_{k}=e_{k}+\mathcal{K}_{k}$ where $e_{k}$, and
$\mathcal{K}_{k}=\frac{1}{2}\bm{u}_{k}\bm{\cdot}\bm{u}_{k}$ are the internal
energy and kinetic energy, respectively.

The variables with the subscript “I” represent the variables at interfaces,
for which there are several possible definitions [18, 20, 25]. Here we choose
the following

$\displaystyle\bm{u}_{I}={\bm{u}}={\sum y_{k}\bm{u}_{k}},\quad
p_{I}=\sum\alpha_{k}p_{k},$
$\displaystyle\quad{\overline{\overline{\tau}}}_{I}=\sum\alpha_{k}{\overline{\overline{\tau}}}_{k},\quad\overline{\overline{T}}_{I}=-p_{I}\overline{\overline{I}}+\overline{\overline{\tau}}_{I},$

where $y_{k}$ denotes the mass fraction $y_{k}=\alpha_{k}\rho_{k}/\rho$, and
${\bm{u}}$ is the mass-fraction weighted mean velocity. Here, we have used the
mass weighted velocity $\bm{u}$ to approximate the interface velocity
$\bm{u}_{I}$. Other possible definitions for $\bm{u}_{I}$ as convex
combination of component velocities will result in the same reduced model
since $(\bm{u}_{I}-\bm{u})=\mathcal{O}(\text{Kn})$, as to be shown in Section
II.3.

The interfacial stress $\overline{\overline{\tau}}_{I}$ is defined in such way
that the thermodynamical laws are respected. In the present work we do not
expand the issue on defining the interfacial variables since they are included
in the terms of order $\mathcal{\text{Kn}}^{2}$ that are to be omitted.

The inter-phase exchange terms include the velocity relaxation
$\mathcal{M}_{k}$, the pressure relaxation $\mathcal{F}_{k}$, and the
temperature relaxation $\mathcal{Q}_{k}$. They are as follows:

$\displaystyle\mathcal{M}_{k}=\vartheta\left(\bm{u}_{k^{*}}-\bm{u}_{k}\right),$
$\displaystyle\mathcal{F}_{k}=\varsigma\left({p}_{k}-{p}_{k^{*}}\right),$
$\displaystyle\mathcal{Q}_{k}=\eta\left(T_{k^{*}}-T_{k}\right),$ (2)

where $k^{*}$ denotes the conjugate component of the $k$-th component, i.e.,
$k=1,\;k^{*}=2$ or $k=2,\;k^{*}=1$. The relaxation velocities are all positive
$\vartheta>0,\;\eta>0,\;\varsigma>0$.

The phase stress tensor, $\overline{\overline{T}}_{k}$, can be written as

$\overline{\overline{T}}_{k}=-p_{k}\overline{\overline{I}}+\overline{\overline{\tau}}_{k}.$
(3)

For the viscous part we use the Newtonian approximation

$\overline{\overline{\tau}}_{k}=2\mu_{k}\overline{\overline{D}}_{k}+\left(\mu_{b,k}-\frac{2}{3}\mu_{k}\right)\nabla\bm{\cdot}\bm{u}_{k},$
(4)

where $\mu_{k}>0$ is the coefficient of shear viscosity and $\mu_{b,k}>0$ is
the coefficient of bulk viscosity. The tensor $\overline{\overline{D}}_{k}$ is
the deformation rate, which takes the following form

$\overline{\overline{D}}_{k}=\frac{1}{2}\left(\nabla\bm{u}_{k}+\left(\nabla\bm{u}_{k}\right)^{\text{T}}\right).$

The heat conduction is represented by
$q_{k}=\nabla\bm{\cdot}\left(\lambda_{k}\nabla T_{k}\right)$. The term
$\mathcal{I}_{k}$ is a volume energy source term, which can be used to model
the laser energy adsorption in direct drive ICF.

Performing straightforward calculus similar to that of Murrone [26], one can
derive the following equation for the internal energy

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}e_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{u}_{k}\right)=-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{u}_{k}$
$\displaystyle-
p_{I}\mathcal{F}_{k}+\left(\bm{u}-\bm{u}_{k}\right)\bm{\cdot}\mathcal{M}_{k}$
$\displaystyle+\left(\bm{u}_{k}-\bm{u}\right)\bm{\cdot}\left(\overline{\overline{T}}_{I}\bm{\cdot}\nabla\alpha_{k}\right)+\mathcal{G}_{k},$
(5)

where
$\mathcal{G}_{k}=q_{k}+\mathcal{I}_{k}+\mathcal{S}_{k}+\mathcal{Q}_{k},\;\;\mathcal{S}_{k}=\alpha_{k}\overline{\overline{\tau}}_{k}:\overline{\overline{D}}_{k}$.

###  II.2 Splitting diffusion-related terms

In continuum mechanics the mass diffusion strength is described through the
difference between the phase velocity and the velocity of the mass center. We
introduce the following definition of the diffusion velocity

$\bm{w}_{k}=\bm{u}_{k}-\bm{u}.$ (6)

For future use, we gather the diffusion related terms separately. To do this,
we replace $\bm{u}_{k}$ by $\bm{u}+\bm{w}_{k}$ to obtain the following
reformulation of eq. 1:

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}}}{\partial{t}}+\nabla\bm{\cdot}(\alpha_{k}\rho_{k}{\bm{u}})=-\nabla\bm{\cdot}\bm{J}_{k},$
(7a)
$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}\bm{u}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}\bm{u}\bm{u}-\alpha_{k}\overline{\overline{T}}_{ak}\right)=-\overline{\overline{T}}_{I}\bm{\cdot}\nabla\alpha_{k}$
$\displaystyle+\mathcal{M}_{k}+\frac{\partial{\alpha_{k}\rho_{k}\bm{w}_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(2\alpha_{k}\rho_{k}\bm{u}\bm{w}_{k}+\alpha_{k}\rho_{k}\bm{w}_{k}\bm{w}_{k}\right)$
$\displaystyle+\nabla\bm{\cdot}\left(\alpha_{k}\overline{\overline{T}}_{wk}\right),$
(7b)
$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}{{E}}_{ak}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}{{E}}_{ak}\bm{u}-\alpha_{k}\overline{\overline{T}}_{ak}\bm{\cdot}\bm{u}\right)=$
$\displaystyle-\bm{u}_{I}\bm{\cdot}\left(\overline{\overline{T}}_{I}\bm{\cdot}\nabla{\alpha_{k}}\right)+\bm{u}\mathcal{M}_{k}-p_{I}\mathcal{F}_{k}+\mathcal{Q}_{k}+q_{k}+\mathcal{I}_{k}$
$\displaystyle{\color[rgb]{0,0,0}-}\frac{\partial{\alpha_{k}\rho_{k}{E}_{wk}}}{\partial{t}}-\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}\left(\bm{u}E_{wk}+{E}_{ak}\bm{w}_{k}+E_{wk}\bm{w}_{k}\right)\right)$
$\displaystyle+\nabla\bm{\cdot}\left(\alpha_{k}\overline{\overline{T}}_{ak}\bm{\cdot}\bm{w}_{k}+\alpha_{k}\overline{\overline{T}}_{wk}\bm{\cdot}\bm{u}+\alpha_{k}\overline{\overline{T}}_{wk}\bm{\cdot}\bm{w}_{k}\right),$
(7c)
$\displaystyle\frac{\partial{\alpha_{k}}}{\partial{t}}+\bm{u}\bm{\cdot}\nabla\alpha_{k}=\mathcal{F}_{k},$
(7d)

where the diffusion flux

$\bm{J}_{k}=\alpha_{k}\rho_{k}\bm{w}_{k}.$

We have used the following decomposition

$E_{k}=E_{ak}+E_{wk},\quad\overline{\overline{T}}_{k}=\overline{\overline{T}}_{ak}+\overline{\overline{T}}_{wk},$

The quantities with subscript “$ak$” contains only the mass centered velocity
$\bm{u}$, i.e.,

$E_{ak}=e_{k}+\frac{1}{2}\bm{u}\bm{\cdot}\bm{u},$
$\overline{\overline{T}}_{ak}=-p_{k}\overline{\overline{I}}+2\mu_{k}\overline{\overline{D}}_{ak}+\left(\mu_{b,k}-\frac{2}{3}\mu_{k}\right)\nabla\bm{\cdot}\bm{u},$
$\overline{\overline{D}}_{k}=\overline{\overline{D}}_{ak}+\overline{\overline{D}}_{wk},\quad\overline{\overline{D}}_{ak}=\frac{1}{2}\left(\nabla\bm{u}+\left(\nabla\bm{u}\right)^{\text{T}}\right).$

The terms with $\bm{w}_{k}$ can be regarded as diffusion-induced. Similarly,
eq. 5 can be reformulated as

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}e_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{u}\right)=-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{u}-p_{I}\mathcal{F}_{k}$
$\displaystyle+\bm{w}_{k}\bm{\cdot}\mathcal{M}_{k}+\bm{w}_{k}\bm{\cdot}\left(\overline{\overline{T}}_{I}\bm{\cdot}\nabla\alpha_{k}\right)-\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{w}_{k}\right)$
$\displaystyle-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{w}_{k}+\mathcal{G}_{k}.$ (8)

###  II.3 Reduction of the BN model

Continuum assumption is usually accepted in the case of a small Knudsen
number, i.e.,

$Kn=\frac{\lambda}{\Delta}\leq 0.001<<1.$ (9)

For moderate Mach number and collision of ions with comparable masses, it can
be shown that[27]

$\frac{|\bm{w}_{k}|}{u_{shock}}\approx\frac{\lambda}{\Delta}<<1.$ (10)

With such scale estimation, in the following we will drop the terms of order
$\mathcal{O}(|\bm{w_{k}}^{2}|)$. The velocity relaxation $\mathcal{M}_{k}$ in
the model (1) is related to the ion friction
$\mu_{kk^{*}}\nu_{i}n_{k}\left(\bm{w}_{k}-\bm{w}_{k^{*}}\right)$, where
$\mu_{kk^{*}}$ is the reduced mass, $\nu_{i}$ is the ion collision frequency,
and $n_{k}$ is the number density. This term is finite under the concerned
scenario. Thus, the velocity relaxation rate is estimated to be
$\mathcal{O}\left(1\right)$. The general idea is to reduce the Baer-Nunziato
model (1) via order analysis. More concretely, we will drop
$\mathcal{O}(|\bm{w_{k}}^{2}|)$ terms and close $\mathcal{O}(|\bm{w_{k}}|)$
terms with the established diffusion laws.

The diffusion velocity is defined with the diffusion laws such as the Fick’s
law, resulting in the redundance of the model (7b). Thus we can retain only
one momentum equation, i.e., the sum of eq. 7b

$\frac{\partial{\rho\bm{u}}}{\partial{t}}+\nabla\bm{\cdot}\left(\rho\bm{u}\bm{u}-\sum\alpha_{k}\overline{\overline{T}}_{ak}\right)=\sum\nabla\bm{\cdot}\left(\alpha_{k}\overline{\overline{T}}_{wk}\right).$
(11)

With such evaluation, the equation for the internal energy section II.2 is
reduced to

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}e_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{u}\right)=-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{u}-p_{I}\mathcal{F}_{k}$
$\displaystyle-\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{w}_{k}\right)-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{w}_{k}+\mathcal{G}_{k},$
(12)

By summing eq. 7c and abandoning terms of $\mathcal{O}(|\bm{w}_{k}|^{2})$, one
can obtain

$\displaystyle\frac{\partial{\rho{{E}}_{a}}}{\partial{t}}+\nabla\bm{\cdot}\left(\rho{{E}}_{a}\bm{u}-\overline{\overline{T}}_{a}\bm{\cdot}\bm{u}\right)=\sum\left(\mathcal{Q}_{k}+q_{k}+\mathcal{I}_{k}\right)$
$\displaystyle-\nabla\bm{\cdot}\left(\sum\alpha_{k}\rho_{k}e_{k}\bm{w}_{k}\right)$
$\displaystyle+\nabla\bm{\cdot}\left(\sum\alpha_{k}\overline{\overline{T}}_{ak}\bm{\cdot}\bm{w}_{k}+\sum\alpha_{k}\overline{\overline{T}}_{wk}\bm{\cdot}\bm{u}\right),$
(13)

where

$E_{a}=\sum y_{k}E_{ak}.$

Note that the enthalpy diffusion flux comes from the second and third terms on
the right hand side of eq. 13.

After the above reduction we can have the following closed system consisting
of the phase mass equations (7a), the mixture momentum equation (11), the
phase internal energy equation (II.3) and the volume fraction equation (1d).
For clarity, we present the obtained model as follows

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}}}{\partial{t}}+\nabla\bm{\cdot}(\alpha_{k}\rho_{k}{\bm{u}})=-\nabla\bm{\cdot}\bm{J}_{k},$
(14a)
$\displaystyle\frac{\partial{\rho\bm{u}}}{\partial{t}}+\nabla\bm{\cdot}\left(\rho\bm{u}\bm{u}-\sum\alpha_{k}\overline{\overline{T}}_{ak}\right)=\sum\nabla\bm{\cdot}\left(\alpha_{k}\overline{\overline{T}}_{wk}\right),$
(14b)
$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}e_{k}}}{\partial{t}}+\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{u}\right)=-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{u}$
$\displaystyle-
p_{I}\mathcal{F}_{k}-\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}e_{k}\bm{w}_{k}\right)-\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{w}_{k}+\mathcal{G}_{k},$
(14c)
$\displaystyle\frac{\partial{\alpha_{k}}}{\partial{t}}+\bm{u}\bm{\cdot}\nabla\alpha_{k}=\mathcal{F}_{k}.$
(14d)

Note that the energy flux caused by mass diffusion is represented by
$\alpha_{k}\rho_{k}e_{k}\bm{w}_{k}$. Energy diffusion effects of heat
conduction and viscous dissipation is included in the term $\mathcal{G}_{k}$
(see the definition in eq. 5).

The absence of a mass diffusion term in the volume fraction equation (14d) is
due to the particular choice of interfacial velocity $\bm{u}_{I}=\bm{u}$. In
more general case we can reformulate eq. 1d as follows

$\frac{\partial{\alpha_{k}}}{\partial{t}}+\bm{u}\bm{\cdot}\nabla\alpha_{k}=(\bm{u}-\bm{u}_{I})\bm{\cdot}\nabla\alpha_{k}+\mathcal{F}_{k},$

Here, one can see that the first term on the right hand side is actually of
order $\mathcal{O}(\text{Kn}^{2})$ as long as the $\bm{u}_{I}$ is of the same
order as $\bm{u}_{k}$. This means that even in this general situation eq. 14d
holds in the context our approximation.

In the framework of the fractional step method, the hydrodynamic subsystem of
(14) is non-conservative due to the equation for internal energy and the
volume fraction. To alleviate the non-conservativeness, we use the mixture
energy equation (13) as an auxiliary equation in solving the hyperbolic sub-
system as the six equation model [28]. An instantaneous pressure relaxation
follow after solving the hyperbolic sub-system, thus resulting in a pressure-
equilibrium model. Such approach ensures robustness in numerical
implementation with marginal sacrifice of computation efficiency.

Note that eq. 13 is not a strict consequence of eq. 14 since the latter loses
a momentum equation after reduction. However, eq. 13 is only invoked in
solving the hyperbolic part where the consistency remains.

We now check the entropy dissipative property of the model. The full entropy
equation before reducetion can be derived in a way similar to [29, 26], which
read

$\displaystyle\alpha_{k}\rho_{k}T_{k}\frac{\mathrm{D}_{k}s_{k}}{\mathrm{D}t}=\left(p_{k}-p_{I}\right)\mathcal{F}_{k}+\mathcal{G}_{k}$
(15)
$\displaystyle+p_{k}\bm{w}_{k}\cdot\nabla\alpha_{k}+\bm{w}_{k}\bm{\cdot}\overline{\overline{T}}_{I}\bm{\cdot}\nabla\alpha_{k}+\left(\bm{u}-\bm{u}_{k}\right)\cdot\mathcal{M}_{k},$

The last three terms are evaluated to be of $\mathcal{O}(\text{Kn}^{2})$ and
therefore abandoned. Thus, the entropy equation is reduced to the following

$\displaystyle\alpha_{k}\rho_{k}T_{k}\frac{\mathrm{D}_{k}s_{k}}{\mathrm{D}t}=\left(p_{k}-p_{I}\right)\mathcal{F}_{k}+\mathcal{G}_{k}.$
(16)

The model (14) keeps the equations for the internal energy, and abandoning a
momentum equation has no impact on the entropy equation. Therefore, the above
reduced entropy equation can also be derived from the reduced model (14). With
eq. 16, one can readily prove the entropy inequality in the absence of
external heat flux and energy source

$\displaystyle\sum_{k=1}^{2}\alpha_{k}\rho_{k}\frac{\mathrm{D}_{k}s_{k}}{\mathrm{D}t}+\sum_{k=1}^{2}\nabla\bm{\cdot}\left(\frac{\bm{q}_{k}}{T_{k}}\right)-\sum_{k=1}^{2}\frac{\mathcal{I}_{k}}{T_{k}}\geq
0.$ (17)

The proof can be performed in a totally similar manner to that in our previous
publication[29], and we provide a brief proof in the appendix.

In the model (14) all the mass diffusion effects are gathered into the terms
containing $\bm{w}_{k}$. In comparision with the model in literature[30], our
model contains the mass diffusion contribution to the viscous stress.
Moreover, as the original Baer-Nunziato model, this model includes all the
disequilibrium effects in phase pressure and temperature, which is driven
towards equilibrium with the corresponding relaxation terms.

The derived model has certain advantages over the widely used model of Cook
[30] in the following aspects: (a) It retains the pressure/temperature
disequilibria between phases through the relaxation terms, thus allowing us to
consider the disequilibrium effects between components. Such disequilibrium is
significant in some phenomena of plasma physics. (b) It frees the hydrodynamic
step of the temperature equilibrium constraint, which is the root of spurious
oscillations for diffuse interface problems[31]. Some special numerical cures
exist in literature[32, 33]. (c) As to the constitutive law for viscous
stress, the viscous stress of each component is determined with their own
velocity (derived by virture of the diffusion law), which is compatible with
thermodynamic relations. In the model of Cook [30], only the mass weighted
velocity is used for the evaluation of viscous stress, which is decoupled from
the diffusion effect. These two approaches may lead to noticeable difference
in simulation. A numerical test is to be considered in Section IV to
demonstrate this issue.

###  II.4 Diffusion models

To make the model solvable, some closure relations are needed to relate the
diffusion terms to basic variables (including density, temperature, and
pressure). The diffusion coefficients are usually derived from kinetic
theories or from experimental measurements. For ordinary neural flows without
plasma, there are some well accepted equations to describe the diffusion
processes, for example, the Fourier’s law for heat flux, the Newton’s law for
viscous stress, and the Fick’s law for mass diffusion. However, for plasma
flows, the diffusion laws are far more complicated. The electron-ion
disequilibrium may also play an important role. For simplicity, in the present
paper we concentrate on the case of electron-ion equilibrium and show the
feasibility of our approach to consider multi-component plasma flows with
diffusions. This assumption should not limit the present model’s applicability
to the disequilibrium ion-electron temperature flows. Further work on the ion-
electron disequilibrium plasma flows will follow in our next publication.
Moreover, one can also regard our model as the limit of instantaneous ion-
electron relaxation time of some ion/electron temperature disequilibrium
model.

Based on the above discussions, we use the ion-electron equilibrium Spitzer-
Harm model[34] for the heat flux. The commonly used models to calculate the
plasma viscosity and mass diffusion include the Clerouin’s model[35] for
viscosity and the Paquette’s model[36] for mass diffusion, respectively. In
recent years, driven by the need to evaluate species mixing/separation in
thermonuclear inertial confinement fusion plasmas, a series of transport
models [37, 38, 39] for multicomponent plasmas are derived. In these models
the mass diffusion is driven by the gradients of field variables such as
species concentration, ion/electron pressure (baro-diffusion) and ion/electron
temperature (thermal diffusion). For Simakov[40], the closure law for the
diffusion velocity is written as follows:

$\bm{w}_{k}=-\sum
D_{kj}\bm{d}_{j}+D_{k}^{T}\nabla\left(\text{ln}T_{i}\right),$ (18)

and

$\displaystyle\bm{d}_{k}=\nabla
x_{k}+\left(x_{k}-y_{k}\right)\nabla\left(\text{ln}p\right)+\left(z_{k}-y_{k}\right)\frac{\nabla
p_{e}}{p}$
$\displaystyle+\left(\frac{Z_{k}n_{k}}{n_{e}}-\frac{Z_{k}^{2}n_{k}}{\sum_{j}Z_{j}^{2}n_{j}}\right)\frac{n_{e}}{n}\frac{\beta_{0}\nabla
T_{e}}{T},$ (19)

where $x_{k}$, $y_{k}$, $z_{k}$, $Z_{k}$ and $n_{k}$ are the number fraction,
the mass fraction, the charge fraction, the charge number and the number
density of component $k$, respectively.

The parameter $\beta_{0}$ is a function of effective charge number,

$\beta_{0}\left(Z_{eff}\right)=\frac{30Z_{eff}\left(11Z_{eff}+15\sqrt{2}\right)}{217Z_{eff}^{2}+604\sqrt{2}Z_{eff}+288}.$
(20)

In the case of two components, the above representation is identical to that
of Kagan & Tang [27]:

$\bm{w}_{k}=-\frac{D}{y_{k}}\left(\nabla
y_{k}+D_{pk}\nabla\text{log}p\right),$ (21)

where we have omitted the electro- and thermo-diffusions. The baro-diffusion
coefficient $D_{pk}$ is determined as follows

$D_{pk}=y_{k}y_{k*}\left(M_{k*}-M_{k}\right)\left(\frac{y_{1}}{M_{1}}+\frac{y_{2}}{M_{2}}\right).$
(22)

where $M_{k}$ is the ion mass of the $k$ component.

##  III Numerical method

The model (14) is solved by using the splitting procedure. According to
physical processes, the model is split into three sub-systems: the hyperbolic
sub-system, the diffusion sub-system and the the relaxation sub-system. These
sub-systems are solved within each time step in order. The solution of one
sub-system serves as the initial condition for the next. The numerical methods
for the hyperbolic and parabolic part are similar to those described in our
previous works [29], thus are omitted here. Particularities here consist in
the temperature relaxation and mass diffusion, which are to be described in
this section.

###  III.1 Thermal relaxation

Here we describe some details for the solution of the thermal relaxation sub-
system, which reads

$\frac{\partial{\alpha_{k}\rho_{k}e_{k}}}{\partial{t}}=\mathcal{Q}_{k}-p_{I}\frac{\partial{\alpha_{k}}}{\partial{t}}.$
(23)

To account for the interface motion due to arising pressure disequilibrium, we
adopt the following assumption

$\frac{\partial{\alpha_{k}}}{\partial{t}}=\Lambda\mathcal{Q}_{k}$ (24)

where the exchanged energy $\mathcal{Q}_{k}=\eta\left(T_{k}^{*}-T_{k}\right)$
and the interfacial pressure is approximated as $p_{I}=\sum{\alpha_{k}p_{k}}$.

The volume fraction varies in such a way that the phasic pressure equilibrium
is maintained, implicitly meaning that the pressure relaxation rate is much
larger than the thermal relaxation rate. By using the pressure equilibrium
condition

$\frac{\partial{p_{1}}}{\partial{t}}=\frac{\partial{p_{2}}}{\partial{t}},$

one can derive the evolution equation for $\alpha_{k}$ in the case of
polytropic gas EOS as follows:

$\Lambda=\frac{1}{p_{I}+\frac{\sum{p_{k}/\alpha_{k}}}{\sum{\left(\gamma_{k}-1\right)/\alpha_{k}}}}.$
(25)

Note that partial densities and momentum remain unchanged during the thermal
relaxation stage. Reformulation of eqs. 23, 24 and 25 gives

$\frac{\partial{T_{k}}}{\partial{t}}=\widehat{\eta}_{k}\left(T_{k}^{*}-T_{k}\right),$
(26)

where

$\widehat{\eta}_{k}=\frac{\frac{\sum{p_{k}/\alpha_{k}}}{\sum{\left(\gamma_{k}-1\right)/\alpha_{k}}}}{\left(\frac{\sum{p_{k}/\alpha_{k}}}{\sum{\left(\gamma_{k}-1\right)/\alpha_{k}}}+p_{I}\right)m_{k}C_{vk}}\eta.$

For the solution of the nonlinear eqs. 26 and 24, iterative methods should be
used, where for each iteration one solves the following linearized ODEs in
each time step $\Delta t$

$\frac{\partial{T_{k}^{(s+1)}}}{\partial{t}}=\widehat{\eta}_{k}^{(s)}\left(T_{k}^{*(s+1)}-T_{k}^{(s+1)}\right),$
(27)

where $(s)$ denotes the iteration index. The analytical solution for the above
ODEs are as follows

$\displaystyle
T_{1}^{(s+1)}=A+B\widehat{\eta}_{1}^{(s)},\;\;T_{2}^{(s+1)}=A-B\widehat{\eta}_{2}^{(s)}$
$\displaystyle
A=\frac{T_{10}\widehat{\eta}_{2}^{(s)}+T_{20}\widehat{\eta}_{1}^{(s)}}{\widehat{\eta}_{2}^{(s)}+\widehat{\eta}_{1}^{(s)}},\;\;B=\frac{\left(T_{10}-T_{20}\right)C}{\widehat{\eta}_{2}^{(s)}+\widehat{\eta}_{1}^{(s)}},$
$\displaystyle
C=e^{-\left(\widehat{\eta}_{1}^{(s)}+\widehat{\eta}_{2}^{(s)}\right)\Delta
t},$

where $T_{10}$ and $T_{20}$ are the initial component temperature at the
beginning of the temperature relaxation stage.

The iterations are performed until the convergence condition
$|T_{k}^{(s+1)}-T_{k}^{(s)}|<\epsilon$ is satisfied.

###  III.2 Mass diffusion

The sub-system for the mass diffusion reads:

$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}}}{\partial{t}}$
$\displaystyle=$
$\displaystyle-\nabla\bm{\cdot}\left(\alpha_{k}\rho_{k}\bm{w}_{k}\right):=\mathcal{C}_{k},$
$\displaystyle\frac{\partial{\alpha_{k}\rho_{k}e_{k}}}{\partial{t}}$
$\displaystyle=$ $\displaystyle-\alpha_{k}\rho_{k}\bm{w}_{k}\bm{\cdot}\nabla
e_{k}$ $\displaystyle-$
$\displaystyle\alpha_{k}p_{k}\nabla\bm{\cdot}\bm{w}_{k}-p_{I}\frac{\partial{\alpha_{k}}}{\partial{t}}:=\mathcal{E}_{k}-p_{I}\frac{\partial{\alpha_{k}}}{\partial{t}}.$

In similar manner to that for defining the evolution equation in Section
III.1, the equation describing the volume fraction variation under the
pressure equilibrium takes the following form:

$\frac{\partial{\alpha_{1}}}{\partial{t}}=\frac{\left(\frac{G_{1}\mathcal{E}_{1}}{\alpha_{1}}-\frac{G_{2}\mathcal{E}_{2}}{\alpha_{2}}\right)+\left(\frac{p_{1}\mathcal{C}_{1}}{\alpha_{1}\rho_{1}}-\frac{p_{2}\mathcal{C}_{2}}{\alpha_{2}\rho_{2}}\right)}{\left(\frac{G_{1}}{\alpha_{1}}+\frac{G_{2}}{\alpha_{2}}\right)p_{I}+\left(\frac{p_{1}}{\alpha_{1}}+\frac{p_{2}}{\alpha_{2}}\right)},$
(28)

where $G_{k}=\gamma_{k}-1$.

Note that the diffusion term in eq. 28 is to maintain the obtained pressure
equilibrium. This does not contradicts eq. 14d, but is a result of the
sequence in solving the split sub-systems. Similar numerical strategy is
adopted in solving phase transition problems[41].

The involved spatial derivatives are approximated with the central difference
scheme.

##  IV Numerical results

In this section we present some numerical results to validate the proposed
model and numerical methods. Moreover, with the aid of the proposed method we
investigate the impact of thermal relaxation on RT instability development.

###  IV.1 The mass diffusion problem

#### The pure diffusion problem

Let us consider a pure diffusion problem as in [42, 43]. The two components
are characterized by the polytropic EOS with the adiabatic coefficients
$\gamma_{1}=2.0$ and $\gamma_{2}=1.4$ and densities $\rho_{1}=20.0$ and
$\rho_{2}=1.0$, respectively. In the computational domain $[0,1]$ the fluids
are initially in temperature and pressure equilibrium which means
$\left(\gamma_{k}-1\right)\rho_{k}C_{vk}=\text{const}$. This relation gives a
constraint for prescribing the heat capacities $C_{vk}$.

The initial mixture density and partial density are given as follows:

$\displaystyle\rho=\frac{1}{2}\left(\rho_{1}+\rho_{2}\right)-\frac{1}{2}\left(\rho_{1}-\rho_{2}\right)\text{erf}\left(z\right),$
$\displaystyle\rho
y_{1}=\frac{1}{2}\rho_{1}-\frac{1}{2}\rho_{1}\text{erf}\left(z\right),$ (29)

where

$z=\frac{x-x_{0}}{\sqrt{4Dt+h_{0}^{2}}}.$

In the present test we use $t=0.0$, $x_{0}=0.5$ and $h_{0}=0.02$ for
prescribing the initial condition.

With the pressure as large as $p=1\times 10^{5}$, the Mach number is so small
that the compressibility effect can be neglected. In this case given the
mixture density profile one can derive the mass weighted mean velocity

$u=-\frac{D}{\rho}\frac{\partial{\rho}}{\partial{x}}.$

Moreover, the analytical solution for density to this pure diffusion problem
is given by section IV.1 [42, 44].

Computations are performed to the time moment $t=0.5$ on a series of refining
grids of 32, 64, 128, 256, 512 and 1024 cells. Reflective boundary conditions
are imposed on both sides. The errors for density are defined as its distance
to the analytical solution section IV.1. The convergence performance is
displayed in Figure 3. One can see that the second order is reached as
expected. The corresponding convergence performance for different variables
are demonstrated in Figure 4.

#### The advection-diffusion problem

We continue to consider the advection-diffusion problem where the multi-
component fluid is transported by a uniform velocity $u=4.0$ while diffusing.
By choosing a reference moving at the transport velocity, one can see that the
analytical solution is still described by section IV.1, only with a
transported interface center. The computational domain is enlarged to [0,4]
and the interface center is transported to $x=2.5$ at $t=0.5$. Similarly, we
obtain the convergence rate for this problem as displayed in Figure 3. The
convergence rate is somewhat smaller than that for the pure diffusion problem
since the numerical resolution of the advection part adds to some error.

The comparison between the pure-diffusion and advection-diffusion problems are
displayed in Figure 5. One can observe that on the coarse grid the numerical
results (for the density, velocity and mass/volume fraction) of the advection-
diffusion problem suffer from more deviations from the exact solutions. This
can be explained by the fact that extra numerical dissipation is needed in
solving the advection part. On the other hand, the introduced numerical
dissipation also smooths the pressure/temperature oscillations, as can be seen
in Figure 5(c-d).

(a)

(b)

Figure 3: The convergence rate for the mass diffusion problem. Upper: the pure
diffusion problem, lower: the advection-diffusion problem.

(a) density (b) velocity

(c) pressure (d) temperature

(e) Mass fraction (f) Volume fraction

Figure 4: The numerical results in the pure diffusion problem for density,
velocity, pressure, temperature, mass fraction and volume fraction with the
refinement of grid.

(a) density (b) velocity

(c) pressure (d) temperature

(e) Mass fraction (f) Volume fraction

Figure 5: The numerical results in the pure-diffusion / advection-diffusion
problems for density, velocity, pressure, temperature, mass fraction and
volume fraction.

#### The baro-diffusion and viscous effect

The initial density and velocity profiles are the same as those in the pure-
diffusion problem in Section IV.1. To show the ability of our method to deal
with the baro-diffusion and viscous effect, we consider a problem with
considerable compressibility effect under the background pressure $p=100$.
Under such pressure the compressibility results in a non-uniform pressure
profile, i.e., non-zero pressure gradient. Thus, the baro-diffusion begins to
effect the results. The kinetic viscosity is taken to be $\nu=50D$, which is
large enough to show the impact of the viscous dissipation. Initial and
boundary conditions are the same as the test in Section IV.1.

As for the viscous part, the conventional way is to determine the viscous
stress by using the mass-weighted velocity as in [30]. Different from this
approach, we use the component velocity to calculate their respective viscous
stress. The numerical results obtained by these two approaches are compared in
Figure 6. One can see noticeable difference between these numerical results,
especially in pressure. New extreme arises in the pressure profile obtained
with the $u_{av}$ approach. The difference in pressure results in the
corresponding difference in density.

We assume that the two components are Carbon (C) and Deuterium (D),
respectively. According to eq. 22, the baro-diffusion coefficient for Carbon
is negative ($D_{p1}<0$) and that for the Deuterium is positive ($D_{p2}>0$).
As can be seen from Figure 7, in the neighborhood of the diffuse front, the
pressure gradient is along $-x$. Thus, for Carbon the baro-diffusion flux
$-\rho DD_{p1}\nabla\text{log}p<0$, which is opposite to the mass fraction
gradient driven flux $-\rho D\nabla y_{1}>0$. Therefore, in fig. 7, one can
observe that the Carbon mass fraction is less diffused with the baro-diffusion
effect being included.

(a) Density

(b) Mass fraction

(c) Pressure

Figure 6: The numerical results for the mass diffusion problem with viscous
effect. “Vis. stress with $u_{k}$/$u_{av}$” $-$ numerical results obtained
with the viscous stress being calculated with the component velocity/the mass
weighted velocity.

(a) Density

(b) Mass fraction

(c) Pressure

Figure 7: The numerical results for the mass diffusion problem with/without
the baro-diffusion effect.

###  IV.2 The shock passage through a mixing zone

In this test we demonstrate the capability of the proposed model (14) to deal
with temperature separation phenomenon when the shock travels through the
mixing zone. The mixture consists of two polytropic components with adiabatic
coefficients $\gamma=2.0$ and $\gamma=5/3$. The heat capacity is calculated by
$C_{v}={N_{0}k_{b}(1+Z)}/({\gamma A}),$ where $N_{0}$ is the Avogadro
constant, $k_{b}$ is the Boltzmann constant, $A$ and $Z$ are the atomic weight
and number, respectively. Again, we assume that the two components are Carbon
($A=12{\text{g/mol}},\;Z=6$) and Deuterium ($A=2{\text{g/mol}},\;Z=1$).

The one-dimensional computational domain is of length $L=$40.96$\mu$m. The
initial mixture density is characterized by section IV.1 with $x_{0}=0.6L$,
$h_{0}=L/10$ and $D=0.02$cm2/$\mu$s. A leftward shock of Mach number 5 hits
the mixing interface. The temperature disequilibrium in the initial postshock
zone is neglected since the concentration of the second component is
negligibly small ($y_{2}=1\times 10^{-6}$). The pre-shock mixture is in
temperature and pressure equilibrium with uniform profile $p=5000$Mbar and
$T=100$MK. The component densities can be determined as $\rho_{1}=2.0617$g/cm3
and $\rho_{2}=1.2027$g/cm3 via the EOSs. The initial profiles for the
component temperatures, volume fraction and mass fraction are demonstrated in
Figure 8(a). We first perform computation with the physical temperature
relaxation rate determined with the Coloumb collosion frequency [45]. At the
time moment $t=6\times 10^{-3}$ns the shock travels through the mixing
interface, temperature disequilibrium arises in postshock zone (Figure 8(b)).
The temperature difference can be as large as 525MK. We then increase the
physical temperature relaxation rate by 100 times and perform the same
computation. The corresponding results are displayed in Figure 8(c). It can be
observed that the temperature disequilibrium is reduced. For comparison
purpose we also present the numerical results with the temperature equilibrium
model (or $\eta\to\infty$) in Figure 8(d). The equilibrium temperature lies
between the component temperatures in Figure 8(b-c). The deviation between the
mass fraction and the volume fraction is more obvious in Figure 8(d). This is
because the temperature relaxation leads to the variation of the volume
fraction while it has no impact on the mass fraction, as analysed in Section
III.1. Comparison between the temperature-equilibrium and the temperature-
disequilibrium model indicates that the commonly used temperature-equilibrium
model maybe inadequate for evaluating the temperature relaxation effect in
mixing topology evolution.

(a) Initial data (b) Physical $\eta$

(c) Physical $\eta$ enlarged by 100 times (d) Infinite $\eta$

Figure 8: The variable (component temperatures, the volume fraction and the
mass fraction) distributions after the shock travels through the mixing zone.
Non-dimensional temperatures displayed are $T_{k}/(1500{\text{MK}})$. (a) The
initial profile, (b) the numerical results at $t=6\times 10^{-3}$ns with
physical relaxation rate $\eta_{Phys}$, (c) the numerical results with
$100\eta_{Phys}$, (d) the numerical results with $\eta\to\infty$.

###  IV.3 The RT instability problem under thermal relaxation

In this section we consider a planar RT instability problem in a computational
domain
$\left(x,y\right)\in[0\mu\text{m},10.24\mu\text{m}]\times[0\mu\text{m},40.96\mu\text{m}]$.
The initial condition along $y=L_{y}/2=5.12$cm is demonstrated in Figure 9.
The acceleration is set to be $2.1\times 10^{5}$cm$/\mu$s2, which is equal to
that of a realistic implosion during the deceleration stage on the Omega
facility [46]. The temperature at the interface zone is also within the
temperature range in ICF implosion.

Note that we assume a smeared interface, which maybe a result of various
mixing mechanisms such as molecular diffusion or turbulence. The mixing
interface center is perturbed with a cosine profile
$x_{0}=0.4L_{x}-0.03L_{y}\text{cos}(2\pi y/L_{y})$. Inside the smeared
interface the component temperatures relax towards a equilibrium one. Such
thermal relaxation mechanism has a significant impact on the RT instability
development, as to be shown below by the following direct simulations. The
relaxation rate is determined with the formula in the NRL plasma
formulary[45].

To ensure grid-independence, we first compare the numerical results obtained
on a series of refining grids from 320$\times$80 cells to 2560$\times$640
cells. The corresponding results are displayed in Figure 10. One can see that
the numerical results tend to converge with physical diffusions being
included. The mixing length evolution with time displayed in Figure 12 also
confirms the convergence.

Then we investigate the sensitivity of the mixing length evolution to the
relaxation rate $\eta$ on the $1280\times 320$ grid. The relaxation rate
$\eta$ is taken to be $1\times 10^{5}$, $1\times 10^{6}$ and $\infty$. The
case $\eta=\infty$ correspond to the case where component temperatures relax
instantaneously, i.e., temperature equilibrium.

The density and mass fraction distributions at 0.2ns with different relaxation
rates are shown in Figure 11. The evolution of the mixing length (bubble-to-
spike distance) with time is demonstrated in Figure 12. It can be seen that
the thermal relaxation tend to suppress the growth of the mixing length. The
main mechanism here lies in the unsteady acceleration of the interface caused
by the thermal relaxation. The details of the physical mechanism will be dealt
with in a separate paper and here we focus on the model itself.

Based on the simulation results, we can see that under the ICF deceleration
condition the temperature equilibrium or disequilibrium assumptions lead to
underestimation or overestimation of the mixing length evolution. The
temperature relaxation timescale is comparable to the ICF-concerned timescale,
and the temperature separation is diminishing at a finite rate. This effect
should be considered for the accurate evaluation of the mixing in ICF.

(a)

(b)

Figure 9: The initial condition for the RT instability problem. Top: mixture
density and pressure, bottom: component temperature. Figure 10: The
distribution of the mixture density $\rho$ (Left) and mass fraction $y_{1}$
(right) on series of refining grids (From top to bottom : 320$\times$80,
640$\times$160, 1280$\times$320, 2560$\times$640 ). 10 uniform contours from
0.01 to 0.99 for $y_{1}$ is display on the right. Figure 11: The distribution
of the mixture density $\rho$ (Left) and mass fraction $y_{1}$ (right) with
different relaxation rates $\eta$ at the time moment $t=0.2$ns. From top to
bottom : $\eta=0,\;1\times 10^{6},\;\infty,$ and $\eta$ determined by physical
model (last row). 10 uniform contours from 0.01 to 0.99 for $y_{1}$ is display
on the right.

(a)

(b)

(c)

Figure 12: The evolution of the mixing length with time. Top and middle sub-
figures are the results corresponding to different relaxation rates $\eta$.
The bottom sub-figure shows the grid independence of the numerical results.

## Conclusion

In the present paper we have presented a temperature disequilibrium diffuse-
interface model for compressible multicomponent flows with interphase heat
transfer and diffusions (including viscous, heat conduction and mass
diffusion). The model is reduced from the BN model in the limit of small
Knudsen number ($Kn<<1$) and consists of six equations including phase density
equations, mixture momentum equation, phase internal energy equations and the
volume fraction equation. Velocity difference is closed by the mass diffusion
laws, thus the velocity of each component is available. The viscous stress is
determined by component velocities rather than mixture velocity in literature.
Moreover, the model has included the effect of finite thermal relaxation. We
have described second-order numerical methods for solving the advection-
diffusion part of the proposed model. As for the thermal relaxation, we
propose method to keep the pressure equilibrium. Being equipped with this
model and its solution methods, we have considered a RT instability problem
under ICF deceleration condition with finite thermal relaxation rate. We have
performed a parametric study on the dependence of the mixing length
development on the relaxation rate. Direct numerical simulations demonstrate
that for the RT instability at an interface between the high-density low-
temperature component and the low-density high-temperature component, the
thermal relaxation tend to suppress the development of the instability.
Further details of this mechanism will appear in our future publications.

## Acknowledgement

The present study is supported by National Natural Science Foundation of China
(grant numbers 12205022,11975053).

## Appendix

Here we provide the proof of the entropy inequality eq. 17. Inserting the
pressure relaxation eq. 2 into eq. 17 and performing some reformulations, one
obtains

$\displaystyle\alpha_{1}\rho_{1}\frac{\text{D}_{1}s_{1}}{\text{D}t}+\alpha_{2}\rho_{2}\frac{\text{D}_{2}s_{2}}{\text{D}t}+\nabla\bm{\cdot}\left(\frac{\bm{q}_{1}}{T_{1}}\right)+\nabla\bm{\cdot}\left(\frac{\bm{q}_{2}}{T_{2}}\right)-\frac{\mathcal{I}_{1}}{T_{1}}-\frac{\mathcal{I}_{2}}{T_{2}}$
$\displaystyle=\frac{\alpha_{2}\varsigma(p_{1}-p_{2})^{2}}{T_{1}}+\frac{\alpha_{1}\varsigma(p_{2}-p_{1})^{2}}{T_{2}}+\frac{\eta(T_{2}-T_{1})^{2}}{T_{1}T_{2}}$
$\displaystyle+\frac{\mathcal{S}_{1}}{T_{1}}+\frac{\mathcal{S}_{2}}{T_{2}}+\bm{q}_{1}\bm{\cdot}\nabla\left(\frac{1}{T_{1}}\right)+\bm{q}_{2}\bm{\cdot}\nabla\left(\frac{1}{T_{2}}\right).$
(30)

The first three terms on the right hand side are obviously non-negative. The
term
$\mathcal{S}_{k}=\alpha_{k}\overline{\overline{\tau}}_{k}:\overline{\overline{D}}_{k}$
represents the kinetic energy dissipation due to viscous friction and is non-
negative as long as the viscosity coefficient $\mu_{k}\geq 0$. The non-
negativity of the terms
$\bm{q}_{k}\bm{\cdot}\nabla\left(\frac{1}{T_{k}}\right)$ are ensured by the
Fourier’s law $\bm{q}_{k}=-\alpha_{k}\lambda_{k}\nabla T_{k}$. Thus entropy
does not decrease in the absence of external heat flux and energy source,
which does not contradicts the second law of thermodynamics.

## References

  * Zylstra  _et al._ [2018] A. Zylstra, N. M. Hoffman, H. W. Herrmann, M. Schmitt, Y. Kim, K. Meaney, A. Leatherland, S. Gales, C. Forrest, V. Y. Glebov, _et al._ , Diffusion-dominated mixing in moderate convergence implosions, Physical Review E 97, 061201 (2018).
  * Wilson  _et al._ [1011] D. C. Wilson, P. S. Ebey, T. C. Sangster, W. T. Shmayda, V. Y. Glebov, and R. A. Lerche, Atomic mix in directly driven inertial confinement implosions, Physics of Plasmas 18, 5 (1011).
  * Dimonte and Tipton [2006] G. Dimonte and R. Tipton, K-l turbulence model for the self-similar growth of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities, Physics of Fluids 18, 085101 (2006).
  * Schilling [2021] O. Schilling, Self-similar reynolds-averaged mechanical–scalar turbulence models for Rayleigh-Taylor, Richtmyer-Meshkov, and Kelvin-Helmholtz instability-induced mixing in the small atwood number limit, Physics of Fluids 33, 085129 (2021).
  * Besnard  _et al._ [1992] D. Besnard, F. Harlow, R. M. Rauenzahn, and C. Zemach, Turbulence transport equations for variable-density turbulence and their relationship to two-field models, LAUR-12303 (1992).
  * Grinstein  _et al._ [2021] F. F. Grinstein, J. A. Saenz, and M. Germano, Coarse grained simulations of shock-driven turbulent material mixing, Physics of Fluids 33, 035131 (2021).
  * Olson  _et al._ [2020] R. E. Olson, T. J. Murphy, B. M. Haines, M. R. Douglas, B. J. Albright, M. A. Gunderson, Y. Kim, T. Cardenas, C. E. Hamilton, and R. B. Randolph, Development of the marble experimental platform at the national ignition facility, Physics of Plasmas 27, 102703 (2020).
  * Murphy  _et al._ [2021] T. J. Murphy, B. Albright, M. Douglas, T. Cardenas, J. Cooley, T. Day, N. Denissen, R. Gore, M. Gunderson, J. Haack, _et al._ , Results from single-shock marble experiments studying thermonuclear burn in the presence of heterogeneous mix on the national ignition facility, High Energy Density Physics 38, 100929 (2021).
  * Haines  _et al._ [2020] B. M. Haines, R. Shah, J. Smidt, and B. Abright, Observation of persistent species temperature separation in inertial confinement fusion mixtures, Nature communications 11, 1 (2020).
  * Llor [2005] A. Llor,  _Statistical Hydrodynamic Models for Developed Mixing Instability Flows: Analytical” 0D” Evaluation Criteria, and Comparison of Single-and Two-Phase Flow Approaches_ , Vol. 681 (Springer Science & Business Media, 2005).
  * Llor and Bailly [2003] A. Llor and P. Bailly, A new turbulent two-field concept for modeling Rayleigh-Taylor, Richtmyer-Meshkov, and Kelvin-Helmholtz mixing layers, Laser and Particle Beams 21, 311 (2003).
  * Youngs [1994] D. L. Youngs, Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities, Laser and Particle Beams 12, 725 (1994).
  * Saurel  _et al._ [2012] R. Saurel, G. Huber, G. Jourdan, E. Lapébie, and L. Munier, Modelling spherical explosions with turbulent mixing and post-detonation, Physics of Fluids 24, 115101 (2012).
  * BAE [1986] A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials, International Journal of Multiphase Flow 12, 861 (1986).
  * Kapila  _et al._ [2001] A. Kapila, R. Menikoff, J. Bdzil, S. Son, and D. S. Stewart, Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations, Physics of Fluids 13, 3002 (2001).
  * Lund [2012] H. Lund, A hierarchy of relaxation models for two-phase flow, SIAM Journal on Applied Mathematics 72, 1713 (2012).
  * Saurel  _et al._ [2009a] R. Saurel, F. Petitpas, and R. Berry, Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures, Journal of Computational Physics 228, 1678 (2009a).
  * Saurel and Abgrall [1999] R. Saurel and R. Abgrall, A multiphase godunov method for compressible multifluid and multiphase flows, Journal of Computational Physics 150, 425 (1999).
  * Petitpas and Le Martelot [2014] F. Petitpas and S. Le Martelot, A discrete method to treat heat conduction in compressible two-phase flows, Comput. Therm. Sci. 6, 251 (2014).
  * Perigaud and Saurel [2005] G. Perigaud and R. Saurel, A compressible flow model with capillary effects, Journal of Computational Physics 209, 139 (2005).
  * Drew [1983] D. A. Drew, Mathematical modeling of two-phase flow, Annual Review of Fluid Mechanics 15, 261 (1983).
  * Chinnayya  _et al._ [2004] A. Chinnayya, E. Daniel, and R. Saurel, Modelling detonation waves in heterogeneous energetic materials, Journal of Computational Physics 196, 490 (2004).
  * Abgrall and Saurel [2003] R. Abgrall and R. Saurel, Discrete equations for physical and numerical compressible multiphase mixtures, Journal of Computational Physics 186, 361 (2003).
  * Hantke  _et al._ [2021] M. Hantke, S. Müller, and L. Grabowsky, News on baer–nunziato-type model at pressure equilibrium, Continuum Mechanics and Thermodynamics 33, 767 (2021).
  * Saurel and Pantano [2018] R. Saurel and C. Pantano, Diffuse-interface capturing methods for compressible two-phase flows, Annual Review of Fluid Mechanics 50, 105 (2018).
  * Murrone and Guillard [2005] A. Murrone and H. Guillard, A five equation reduced model for compressible two phase flow problems, Journal of Computational Physics 202, 664 (2005).
  * Kagan and Xianzhu [2014] G. Kagan and T. Xianzhu, Thermo-diffusion in inertially confined plasmas, Pyscial Letters A 378, 1531 (2014).
  * Saurel  _et al._ [2009b] R. Saurel, F. Petitpas, and R. Berry, Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures, Journal of Computational Physics 228, 1678 (2009b).
  * Zhang  _et al._ [2022] C. Zhang, I. Menshov, W. Lifeng, and S. Zhijun, Diffuse interface relaxation model for two-phase compressible flows with diffusion processes, Journal of Computational Physics , 111356 (2022).
  * Cook [2008] A. Cook, Enthalpy diffusion in multicomponent flows, Physics of Fluids 21 (2008).
  * Allaire  _et al._ [2002] G. Allaire, S. Clerc, and S. Kokh, A five-equation model for the simulation of interfaces between compressible fluids, Journal of Computational Physics 181, 577 (2002).
  * Johnsen and Ham [2012] E. Johnsen and F. Ham, Preventing numerical errors generated by interface-capturing schemes in compressible multi-material flows, Journal of Computational Physics 231, 5705 (2012).
  * Williams [2019] R. J. Williams, Fully-conservative contact-capturing schemes for multi-material advection, Journal of Computational Physics 398, 108809 (2019).
  * Spitzer Jr and Härm [1953] L. Spitzer Jr and R. Härm, Transport phenomena in a completely ionized gas, Physical Review 89, 977 (1953).
  * Clérouin  _et al._ [1998] J. Clérouin, M. Cherfi, and G. Zérah, The viscosity of dense plasmas mixtures, EPL (Europhysics Letters) 42, 37 (1998).
  * Paquette  _et al._ [1986] C. Paquette, C. Pelletier, G. Fontaine, and G. Michaud, Diffusion coefficients for stellar plasmas, The Astrophysical Journal Supplement Series 61, 177 (1986).
  * Kagan and Baalrud [2018] G. Kagan and S. Baalrud, Transport formulas for multi-component plasmas within the effective potential theory framework, arXiv , 1611.09872v2 (2018).
  * Simakov and Molvig [2016a] A. Simakov and K. Molvig, Hydrodynamic description of an unmagnetized plasma with multiple ion species. I. general formulation, Phys. Plasmas 23, 032115 (2016a).
  * Balashov and Savenkov [2018] V. Balashov and E. Savenkov, Classical transport equations for burning gas-metal plasmas, Journal of Applied Mechanics & Technical Physics 59, 434 (2018).
  * Simakov and Molvig [2016b] A. N. Simakov and K. Molvig, Hydrodynamic description of an unmagnetized plasma with multiple ion species. i. general formulation, Physics of Plasmas 23, 032115 (2016b).
  * ZEI [2010] Modeling phase transition for compressible two-phase flows applied to metastable liquids, Journal of Computational Physics 229, 2964 (2010).
  * Thornber  _et al._ [2018] B. Thornber, M. Groom, and D. Youngs, A five-equation model for the simulation of miscible and viscous compressible fluids, Journal of Computational Physics 372 (2018).
  * Kokkinakis  _et al._ [2015] I. Kokkinakis, D. Drikakis, D. Youngs, and R. Williams, Two-equation and multi-fluid turbulence models for Rayleigh–Taylor mixing, International Journal of Heat and Fluid Flow 56, 233 (2015).
  * Livescu [2013] D. Livescu, A multiphase model with internal degrees of freedom: application to shock-bubble interaction, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 371, 283 (2013).
  * Richardson [2019] A. Richardson,  _2019 NRL plasma formulary_ (2019).
  * Vold  _et al._ [2021] E. Vold, L. Yin, and B. Albright, Plasma transport simulations of Rayleigh–Taylor instability in near-ICF deceleration regimes, Physics of Plasmas 28, 092709 (2021).


# Note: Gas-liquid crossover in the Lennard-Jones system

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

###### Abstract

It is demonstrated that the crossover between gas- and liquid-like regions on
the phase diagram of the Lennard-Jones system occurs at a fixed value of the
density divided by its value at the freezing point, $\rho/\rho_{\rm fr}\simeq
0.35$. This definition is consistent with other definitions proposed recently.
As a result a very simple practical expression for the gas-to-liquid crossover
line emerges.

In recent years, the possibility to define a demarcation line between liquid-
like and gas-like behaviors of supercritical fluids has been a topic of major
interest Gorelli  _et al._ (2006); Simeoni  _et al._ (2010); McMillan and
Stanley (2010); Brazhkin  _et al._ (2011, 2012a, 2012b, 2013); Gorelli  _et
al._ (2013); Yang  _et al._ (2015); Bryk  _et al._ (2017); Brazhkin  _et al._
(2018); Bryk  _et al._ (2018); Bell  _et al._ (2020); Proctor  _et al._
(2019); Ploetz and Smith (2019); Banuti  _et al._ (2020); Ha  _et al._ (2019);
Maxim  _et al._ (2019); Sun  _et al._ (2020); Khrapak and Khrapak (2021a);
Bell  _et al._ (2021); Khrapak and Khrapak (2021b); Cockrell  _et al._ (2021a,
b). This problem is not just a matter of curiosity. The structural and
dynamical properties are very different in these regimes and quite different
approaches are required for their description. For example, transport
properties in low density (gas-like) fluid regime are dominated by interatomic
collisions and the Chapman-Enskog theory is appropriate Chapman and Cowling
(1990); Lifshitz and Pitaevskii (1995). Atoms move along straight trajectories
between collisions, the collision properties being determined by the pairwise
interaction potential. In dense fluids (liquid-like) regime the transport
properties are dominated by collective interactions between the atoms. As a
result, individual atoms spend a substantial amount of time in local minima of
the multidimensional potential energy surface (a phenomenon often referred to
as “caging”). These local minima migrate with time to allow liquids to flow,
but on considerably slower time scales. It is of great fundamental and
practical interest to understand where the crossover between the gas-like and
liquid-like regimes takes place. This problem is also related to a long-
standing debate about the nature of the supercritical fluid and a more general
question “What is liquid?” Barker and Henderson (1976); Sengers (1979); L.
Woodcock (2017).

Multiple definitions for the gas-to-liquid crossover have been proposed in
recent years Brazhkin  _et al._ (2012a, b); Bell  _et al._ (2020); Khrapak and
Khrapak (2021a). Not all of these definitions were consistent or universal,
which generated a significant amount of debate. From the isomorph theory
perspective Dyre (2014), it seems not unreasonable to assume that a
demarcation line between gas-like and fluid-like regimes should itself be an
approximate isomorph. As such, it should be characterized by a quasi-universal
value of the excess entropy (actual entropy minus entropy of an ideal gas at
the same density and temperature). As we shall see, the excess entropy is
indeed approximately constant along the crosssover line, although it passes
through a near-critical region, where the Pearson correlation coefficient
(which should be greater than $\simeq 0.9$ for Roskilde-simple systems) drops
considerably Dyre (2014); Bell  _et al._ (2019). Properly reduced structural
and dynamical properties should also be quasi-invariant. This point of view
has been elaborated in recent works Bell  _et al._ (2020, 2021). One of the
possible definitions is based on the location of the minimum of the
macroscopically scaled shear viscosity coefficient. This minimum corresponds
to the crossover between the gas-like and liquid-like mechanisms of momentum
transfer, which have quite different nature and hence different asymptotes. It
has been observed that for several different model systems the location of the
minimua in the macroscopically reduced shear viscosity coefficients occurs at
approximately the same value of the excess entropy per particle ($s_{\rm
ex}\simeq-2k_{\rm B}/3$) Bell  _et al._ (2020), and the minimuma values
themselves are also quasiuniversal Khrapak and Khrapak (2022). This is also
where kinetic and potential contributions to the shear viscosity coefficients
are equal to a good accuracy Bell  _et al._ (2021). Thus, the three considered
definitions of a gas-like to liquid-like demarcation line (minima of reduced
viscosity, constancy of excess entropy, and equality of kinetic and potential
contributions to viscosity) nearly coincide with each other Bell  _et al._
(2021).

Independently, it has been recently observed that in various simple fluids
(such as Lennard-Jones, Coulomb, Yukawa and hard-sphere fluids) there exist
two clear asymptotes for the product $D\eta(\Delta/k_{\rm B}T)$, where $D$ is
the self-diffusion coefficient, $\eta$ is the shear viscosity coefficient,
$\Delta=\rho^{-1/3}$ is the mean interatomic separation, $\rho$ is the atomic
density, and $T$ is the temperature Khrapak and Khrapak (2021b). In dense
fluids near the freezing point this product approaches a slightly system-
dependent constant value. This is the regime where the Stokes-Einstein
relation without the hydrodynamic diameter is satisfied (in fact, interatomic
separation plays the role of hydrodynamic diameter) Costigliola  _et al._
(2019); Khrapak (2019); Khrapak and Khrapak (2021b). Far away from the
freezing point, in a rarefied fluid, this product decreases with increasing
density. The intersection of the two asymptotes has been suggested as a
convenient practical condition for the crossover between the gas-like and
liquid-like regions on the phase diagram Khrapak and Khrapak (2021a, b). For
the systems considered, intersection is characterised by very close values of
the reduced excess entropy, $s_{\rm ex}/k_{\rm B}=-0.9\pm 0.1$ Khrapak and
Khrapak (2021b). Similar value was obtained earlier in Ref. Bell  _et al._
(2021), where it was also pointed out that this is nearly the critical point
entropy for simple fluids exhibiting a critical point.

Let us now focus on the Lennard-Jones (LJ) system, a very important model in
condensed matter studies. This model combines relative simplicity with
adequate semi-quantitative approximation of interatomic interactions in real
substances (such as condensed noble gases). The LJ potential is

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

where $\epsilon$ and $\sigma$ are the energy and length scales (LJ units),
respectively. The density and temperature expressed in LJ units are
$\rho_{*}=\rho\sigma^{3}$ and $T_{*}=k_{\rm B}T/\epsilon$.

It has been recently demonstrated that properly reduced transport coefficients
of LJ fluids along isotherms exhibit quasi-universal scaling on the density
divided by its value at the freezing point Khrapak and Khrapak (2021a). In
particular, transport coefficients such as self-diffusion, shear viscosity,
and thermal conductivity, expressed in Rosenfeld’s units Rosenfeld (1999)
remain approximately constant for fixed ratios $\rho/\rho_{\rm fr}$, where
$\rho_{\rm fr}$ is the atomic number density at the freezing point (at a given
temperature). Combined with the excess entropy scaling of the transport
coefficients Rosenfeld (1977, 1999); Dyre (2018); Bell  _et al._ (2019) this
implies that the lines of constant ratio of $\rho/\rho_{\rm fr}$ are likely
characterized by fixed values of excess entropy and thus are isomorphs. The
minima of the reduced shear viscosity and thermal conductivity coefficients
along isotherms correspond to $\rho/\rho_{\rm fr}\simeq 0.2-0.3$ Khrapak and
Khrapak (2021a). The asymptotes of the product $D\eta(\Delta/k_{\rm B}T)$ are
intersecting at $\rho/\rho_{\rm fr}\simeq 0.35$ and this can be considered as
an appropriate definition of the crossover between the gas-like and the
liquid-like regimes Khrapak and Khrapak (2021a, b). The Stokes-Einstein
relation does not apply immediately after the crossover; even higher densities
are required, such that $\rho/\rho_{\rm fr}\gtrsim 0.6$ for the LJ fluid
Khrapak and Khrapak (2021b).

Figure 1: (Color online) Phase diagram of the LJ system on the density-
temperature plane. The squares and stars correspond to the fluid-solid
coexistence boundaries as tabulated in Refs. Sousa  _et al._ (2012) and Hansen
(1970), respectively; the corresponding curves are simple fits from Ref.
Khrapak and Morfill (2011). The liquid-vapour boundary is plotted using the
formulas provided in Ref. Heyes  _et al._ (2019). The circles correspond to
the location of minima of kinematic viscosity Bell  _et al._ (2021). The
crosses denote the points where the contribution to viscosity due to atomic
translation is exactly equal to that due to interatomic interactions Bell  _et
al._ (2021). The dashed curve is plotted using Eq. (2) obtained from the
constant density ratio $\rho/\rho_{\rm fr}=0.35$. A solid line emanating from
the critical point corresponds to constant reduced excess entropy of $s_{\rm
ex}=-0.9k_{\rm B}$, as evaluated from Thol’s et al. equation of state Thol
_et al._ (2016). Inset shows a low temperature portion of the phase diagram
with the focus on near-critical region. All considered definitions of the
crossover between the gas-like and liquid-like regimes of atomic dynamics are
in excellent agreement.

This hypothesis is verified in Fig. 1. It shows the phase diagram of the LJ
system in the ($\rho_{*}$, $T_{*}$) plane. The fluid-solid coexistence data
points are taken from Refs. Sousa  _et al._ (2012); Hansen (1970). The curves
are the fits of the form $T_{*}^{\rm L,S}={\mathcal{A}}^{\rm
L,S}\rho_{*}^{4}-{\mathcal{B}}^{\rm L,S}\rho_{*}^{2}$ (superscripts L and S
correspond to liquid and solid, respectively ). This shape of the fluid-solid
coexistence of LJ systems with constant (or very weakly $\rho_{*}$-dependent)
constants ${\mathcal{A}}$ and ${\mathcal{B}}$ is a very robust result
reproduced in a number of various theories and approximations Pedersen  _et
al._ (2016); Rosenfeld (1976); Khrapak and Morfill (2011); Heyes  _et al._
(2015); Khrapak and Ning (2016); Costigliola  _et al._ (2016); Khrapak (2020);
Heyes  _et al._ (2021). Here we take constant values of ${\mathcal{A}}=2.29$
and ${\mathcal{B}}=0.71$ at freezing proposed in Ref. Khrapak and Morfill
(2011). The liquid-vapour 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.694$ Sousa  _et al._ (2012) and
$T_{*}^{\rm cr}\simeq 1.326$ Heyes  _et al._ (2019), respectively. Additional
symbols appearing in the supercritical region are: The circles correspond to
the location of minima of kinematic viscosity Bell  _et al._ (2021); the
crosses denote the points where the contribution to viscosity due to kinetic
and potential contributions are equal Bell  _et al._ (2021); the solid curve
emanating from vary nearly the critical point and terminating at $T_{*}=9$
corresponds to constant excess entropy of $s_{\rm ex}=-0.9k_{\rm B}$ as
calculated from the equation of state proposed by Thol et al. Thol  _et al._
(2016) in the domain of its applicability. The dashed curve corresponds to the
condition $\rho/\rho_{\rm fr}=0.35$. In view of the relation between the
freezing temperature and density a particularly simple expression for the
crossover line emerges:

$T_{*}=2.29(\rho_{*}/0.35)^{4}-0.71(\rho_{*}/0.35)^{2}.$ (2)

Excellent agreement with other definitions of the gas-to-liquid crossover are
observed from the critical point temperature up to the highest temperatures
investigated ($T_{*}=100$). While different definitions agree, the definition
based on the freezing density is generally more practical, because freezing
density is usually relatively well known (as compared to the exact location of
the minima in reduced transport coefficients, or lines of constant excess
entropy).

We observe in Fig. 1 that the crossover line originates in the vicinity of the
critical point. While this turns out to be quite convenient, this by no means
represents a necessary requirement. The crossover occurs also in systems with
purely repulsive interactions, where liquid-gas phase transition and the
critical point are absent at all Brazhkin  _et al._ (2013); Yang  _et al._
(2015). The crossover line would appear nearly parallel to the freezing line
in the log-log plot. This relationship was pointed out and explained in Ref.
Yang  _et al._ (2015).

To conclude, the crossover between the gas-like and liquid-like regions on the
phase diagram of the Lennard-Jones system corresponds to fixed values of the
density divided by its value at the freezing point and of the reduced excess
entropy. The value $\rho/\rho_{\rm fr}=0.35$ and $s_{\rm ex}\simeq-0.9k_{\rm
B}$, suggested previously from the analysis of the applicability of the
Stokes-Einstein relation, provides very good agreement with other definitions
(based on the minima of kinematic viscosity and equal contribution of the
kinetic and potential terms to the viscosity coefficient). A very simple
explicit expression for the demarcation line is derived. The obtained results
can be useful in identifying gas-to-liquid crossover in liquified noble gases.
Similar arguments are likely applicable to other simple systems beyond LJ.
While particular values of $\rho/\rho_{\rm fr}$ may be different, the values
of $s_{\rm ex}$ are expected to be approximately the same.

I would like to thank V. Brazhkin for useful comments on the manuscript.

## References

  * Gorelli  _et al._ (2006) F. Gorelli, M. Santoro, T. Scopigno, M. Krisch, and G. Ruocco, “Liquidlike behavior of supercritical fluids,” Phys. Rev. Lett. 97, 245702 (2006).
  * Simeoni  _et al._ (2010) G. G. Simeoni, T. Bryk, F. A. Gorelli, M. Krisch, G. Ruocco, M. Santoro, and T. Scopigno, “The Widom line as the crossover between liquid-like and gas-like behaviour in supercritical fluids,” Nature Phys. 6, 503–507 (2010).
  * McMillan and Stanley (2010) P. F. McMillan and H. E. Stanley, “Going supercritical,” Nature Phys. 6, 479–480 (2010).
  * Brazhkin  _et al._ (2011) V. V. Brazhkin, Yu. D. Fomin, A. G. Lyapin, V. N. Ryzhov, and E. N. Tsiok, “Widom line for the liquid–gas transition in Lennard-Jones system,” J. Phys. Chem. B 115, 14112–14115 (2011).
  * Brazhkin  _et al._ (2012a) V. V. Brazhkin, Yu. D. Fomin, A. G. Lyapin, V. N. Ryzhov, and K. Trachenko, “Two liquid states of matter: A dynamic line on a phase diagram,” Phys. Rev. E 85, 031203 (2012a).
  * Brazhkin  _et al._ (2012b) V. V. Brazhkin, A.G. Lyapin, V. N. Ryzhov, K. Trachenko, Y. D. Fomin, and E. N. Tsiok, “Where is the supercritical fluid on the phase diagram?” Phys.-Usp. 182, 1137–1156 (2012b).
  * Brazhkin  _et al._ (2013) V. V. Brazhkin, Yu. D. Fomin, A. G. Lyapin, V. N. Ryzhov, E. N. Tsiok, and K. Trachenko, “Liquid-gas” transition in the supercritical region: fundamental changes in the particle dynamics,” Phys. Rev. Lett. 111, 145901 (2013).
  * Gorelli  _et al._ (2013) F. A. Gorelli, T. Bryk, M. Krisch, G. Ruocco, M. Santoro, and T. Scopigno, “Dynamics and thermodynamics beyond the critical point,” Sci. Rep. 3, 1203 (2013).
  * Yang  _et al._ (2015) C. Yang, V. V. Brazhkin, M. T. Dove, and K. Trachenko, “Frenkel line and solubility maximum in supercritical fluids,” Phys. Rev. E 91, 012112 (2015).
  * Bryk  _et al._ (2017) T. Bryk, F. A. Gorelli, I. Mryglod, G. Ruocco, M. Santoro, and T. Scopigno, “Behavior of supercritical fluids across the Frenkel line,” J. Phys. Chem. Lett. 8, 4995–5001 (2017).
  * Brazhkin  _et al._ (2018) V. V. Brazhkin, C. Prescher, Yu. D. Fomin, E. N. Tsiok, A. G. Lyapin, V. N. Ryzhov, V. B. Prakapenka, J. Stefanski, K. Trachenko, and A. Sapelkin, “Comment on Behavior of supercritical fluids across the Frenkel line,” J. Phys. Chem. B 122, 6124–6128 (2018).
  * Bryk  _et al._ (2018) T. Bryk, F. A. Gorelli, I. Mryglod, G. Ruocco, M. Santoro, and T. Scopigno, “Reply to comment on behavior of supercritical fluids across the Frenkel line,” J. Phys. Chem. B 122, 6120–6123 (2018).
  * Bell  _et al._ (2020) I. H. Bell, G. Galliero, S. Delage-Santacreu, and L. Costigliola, “An entropy scaling demarcation of gas- and liquid-like fluid behaviors,” J. Chem. Phys. 152, 191102 (2020).
  * Proctor  _et al._ (2019) J. E. Proctor, C. G. Pruteanu, I. Morrison, I. F. Crowe, and J. S. Loveday, “Transition from gas-like to liquid-like behavior in supercritical n2,” J. Phys. Chem. Lett. 10, 6584–6589 (2019).
  * Ploetz and Smith (2019) E. A. Ploetz and P. E. Smith, “Gas or liquid? The supercritical behavior of pure fluids,” J. Phys. Chem. B 123, 6554–6563 (2019).
  * Banuti  _et al._ (2020) D.T. Banuti, M. Raju, and M. Ihme, “Between supercritical liquids and gases – reconciling dynamic and thermodynamic state transitions,” J. Supercrit. Fluids 165, 104895 (2020).
  * Ha  _et al._ (2019) M.Y. Ha, T. J. Yoon, T. Tlusty, Y. Jho, and W. B. Lee, “Universality, scaling, and collapse in supercritical fluids,” J. Phys. Chem. Lett. 11, 451–455 (2019).
  * Maxim  _et al._ (2019) F. Maxim, C. Contescu, P. Boillat, B. Niceno, K. Karalis, A. Testino, and C. Ludwig, “Visualization of supercritical water pseudo-boiling at Widom line crossover,” Nature Commun. 10, 4114 (2019).
  * Sun  _et al._ (2020) P. Sun, J. B. Hastings, D. Ishikawa, A. Q. R. Baron, and G. Monaco, “Two-component dynamics and the liquidlike to gaslike crossover in supercritical water,” Phys. Rev. Lett. 125, 256001 (2020).
  * Khrapak and Khrapak (2021a) S. A. Khrapak and A. G. Khrapak, “Transport properties of Lennard-Jones fluids: Freezing density scaling along isotherms,” Phys. Rev. E 103, 042122 (2021a).
  * Bell  _et al._ (2021) I. H. Bell, S. Delage-Santacreu, H. Hoang, and G. Galliero, “Dynamic crossover in fluids: From hard spheres to molecules,” J. Phys. Chem. Lett. 12, 6411–6417 (2021).
  * Khrapak and Khrapak (2021b) S. A. Khrapak and A. G. Khrapak, “Excess entropy and Stokes-Einstein relation in simple fluids,” Phys. Rev. E 104, 044110 (2021b).
  * Cockrell  _et al._ (2021a) C. Cockrell, V. V. Brazhkin, and K. Trachenko, “Universal interrelation between dynamics and thermodynamics and a dynamically driven “c” transition in fluids,” Phys. Rev. E 104, 034108 (2021a).
  * Cockrell  _et al._ (2021b) C. Cockrell, V.V. Brazhkin, and K. Trachenko, “Transition in the supercritical state of matter: Review of experimental evidence,” Phys. Rep. 941, 1–27 (2021b).
  * Chapman and Cowling (1990) S. Chapman and T. G. Cowling,  _The Mathematical Theory of Non-uniform Gases - An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases_ (Cambridge University Press, Cambridge, 1990).
  * Lifshitz and Pitaevskii (1995) E.M. Lifshitz and L. P. Pitaevskii,  _Physical Kinetics - Volume 10_ (Elsevier Science, Stanford, 1995).
  * Barker and Henderson (1976) J. A. Barker and D. Henderson, “What is ”liquid”? Understanding the states of matter,” Rev. Mod. Phys. 48, 587–671 (1976).
  * Sengers (1979) J.M.H. Levelt Sengers, “Liquidons and gasons: controversies about the continuity of states,” Phys. A 98, 363–402 (1979).
  * L. Woodcock (2017) L. Woodcock, “Percolation transitions and fluid state boundaries,” CMST 23, 281 (2017).
  * Dyre (2014) J. C. Dyre, “Hidden scale invariance in condensed matter,” J. Phys. Chem. B 118, 10007–10024 (2014).
  * 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).
  * Khrapak and Khrapak (2022) S. A. Khrapak and A. G. Khrapak, “Minima of shear viscosity and thermal conductivity coefficients of classical fluids,” Phys. Fluids 34, 027102 (2022).
  * 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 (2019) S. Khrapak, “Stokes–Einstein relation in simple fluids revisited,” Mol. Phys. 118, e1643045 (2019).
  * 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).
  * Dyre (2018) J. C. Dyre, “Perspective: Excess-entropy scaling,” J. Chem. Phys. 149, 210901 (2018).
  * 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).
  * Hansen (1970) J.-P. Hansen, “Phase transition of the Lennard-Jones system. II. High-temperature limit,” Phys. Rev. A 2, 221–230 (1970).
  * 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).
  * 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).
  * 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).
  * 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).
  * Rosenfeld (1976) Y. Rosenfeld, “Universality of melting and freezing indicators and additivity of melting curves,” Mol. Phys. 32, 963–977 (1976).
  * Heyes  _et al._ (2015) D. M. Heyes, D. Dini, and A. C. Brańka, “Scaling of Lennard-Jones liquid elastic moduli, viscoelasticity and other properties along fluid-solid coexistence,” Phys. Status Solidi (b) 252, 1514–1525 (2015).
  * Khrapak and Ning (2016) S. A. Khrapak and N. Ning, “Freezing and melting equations for the n-6 Lennard-Jones systems,” AIP Adv. 6, 055215 (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 (2020) S. A. Khrapak, “Lindemann melting criterion in two dimensions,” Phys. Rev. Research 2, 012040 (2020).
  * Heyes  _et al._ (2021) D. M. Heyes, S. Pieprzyk, and A. C. Branka, “Application of cell models to the melting and sublimation lines of the Lennard-Jones and related potential systems,” Phys. Rev. E 104, 044119 (2021).