# On the velocity distribution function of spontaneously evaporating atoms

Sergiu Busuioc  School of Engineering, The University of Edinburgh, Edinburgh,
EH9 3FB, United Kingdom  Livio Gibelli  School of Engineering, The University
of Edinburgh, Edinburgh, EH9 3FB, United Kingdom  Duncan A. Lockerby  School
of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom  James
E. Sprittles  Mathematics Institute, University of Warwick, Coventry CV4 7AL,
United Kingdom

###### Abstract

Numerical solutions of the Enskog-Vlasov (EV) equation are used to determine
the velocity distribution function of atoms spontaneously evaporating into
near-vacuum conditions. It is found that an accurate approximation is provided
by a half-Maxwellian including a drift velocity combined with different
characteristic temperatures for the velocity components normal and parallel to
the liquid-vapor interface. The drift velocity and the temperature anisotropy
reduce as the liquid bulk temperature decreases but persist for relatively low
temperatures corresponding to a vapor behaviour which is only slightly non-
ideal. Deviations from the undrifted isotropic half-Maxwellian are shown to be
consequences of collisions in the liquid-vapor interface which preferentially
backscatter atoms with lower normal-velocity component.

Enskog-Vlasov equation; Evaporation into vacuum; Kinetic boundary conditions;
Drifted anisotropic half-Maxwellian.

##  I Introduction

Evaporation and condensation are ubiquitous processes which play crucial roles
in different areas, from physics and chemistry to biology. The fundamental
understanding of the heat and mass transfer across the liquid-vapor interface
is not only of theoretical interest but also of considerable relevance to a
wide range of applications. Examples include the distillation process in high
vacuum conditions Safarian and Engh (2013); Li and Xu (2014), the heating and
evaporation of fuel droplets Sazhin (2017), the thermal management of modern
electronic devices John  _et al._ (2019), and the flash boiling in gasoline
direct-injection sprays Chang  _et al._ (2020).

The modelling of fluid flows with phase change is challenging since multiple
characteristic time and length scales are involved. The liquid and vapor bulk
phases are well described at the macroscopic scale by hydrodynamic equations,
while the separating region has a more complex structure, with sharp gradients
in flow variables which manifest as jumps at the macroscopic scale. In non-
equilibrium conditions, the region between the two bulk phases comprises two
sub-regions, namely the liquid-vapor interface, whose characteristic length is
of the order of the molecular diameter Rowlinson and Widom (2013), and the
tipically larger Knudsen layer that extends a few mean free paths into the
vapor phase Sone (2000).

The standard kinetic theory studies on evaporation and condensation processes
focus on the vapor dynamics in the Knudsen layer next to the liquid phase. The
liquid-vapor interface is simplified to a structureless surface bounding the
vapor and the molecular exchanges with the liquid phase are dealt with using a
phenomenological boundary condition. More specifically, the molecular flux
coming out of the liquid-vapor surface is assumed to have two contributions:
atoms spontaneously leaving the liquid bulk (evaporated atoms) and atoms which
are backscattered into the vapor after impinging on the liquid-vapor interface
(reflected atoms). The distribution function of evaporating atoms is usually
described by a half-Maxwellian while the Maxwell’s scattering kernel is the
most widely used choice to describe the molecular reflection from the liquid
phase. Two phenomenological coefficients are also introduced to weight the
relative importance of the evaporation and condensation fluxes.

In spite of the increasing number of experimental, theoretical and numerical
studies, the full understanding of evaporation processes is still lacking and
numerous issues remain unresolved. As an example, while the evaporation and
condensation coefficients are generally accepted to be lower than unity,
experimental measurements for water have provided values that span three
orders of magnitude Persad and Ward (2016). Furthermore, a temperature
discontinuity has been experimentally found across the water liquid-vapor
interface that is much larger in magnitude and in the opposite direction to
that predicted by classical kinetic theory or nonequilibrium thermodynamics
Fang and Ward (1999).

A lot of effort has been thus expended in assessing the physical
appropriateness of kinetic boundary conditions at the liquid-vapor interface
mostly through molecular dynamics (MD) simulations and mean-field kinetic
approximations of simple liquids. These studies can be roughly grouped into
two categories depending on whether the liquid evaporation occurs into vapor
Crout (1936); Bird and Liang (2019); Tsuruta  _et al._ (1999); Meland  _et
al._ (2004); Xie  _et al._ (2011, 2012); Kon  _et al._ (2014); Kobayashi  _et
al._ (2017) or into near vacuum conditions Zhakhovskii and Anisimov (1997);
Ishiyama  _et al._ (2004); Frezzotti  _et al._ (2005, 2018); Heinen and Vrabec
(2019). In this latter process the backscattered vapor component is virtually
absent and, therefore, the distribution function of evaporated atoms is
evaluated without any ambiguity.

Accurate MD simulations of liquid argon evaporating into vapor have shown that
atoms leaving the liquid-vapor interface are distributed with excellent
approximation according to a half-Maxwellian while, separately, evaporated and
reflected atoms are not Tsuruta  _et al._ (1999); Meland  _et al._ (2004).
More specifically, the velocity distribution function of evaporated atoms
shows a fatter tail compared to a Maxwellian while the one of reflected atoms
has a reduced tail and larger density in the low velocity region. However,
deviations are described differently, namely through the dependency of the
condensation/evaporation coefficient on the normal-velocity component of atoms
impinging on the liquid surface Tsuruta  _et al._ (1999) or by a drifted half-
Maxwellian Meland  _et al._ (2004).

Using MD simulations, deviations from a half-Maxwellian have also been
revealed in the case of n-dodecane, i.e. more complex molecules, but the
distribution function of evaporated atoms is fitted by an anisotropic half-
Maxwellian (also referred to as bi-Maxwellian) with the temperature normal to
the liquid-vapor interface larger than the parallel one Xie  _et al._ (2011,
2012). The hypothesis that an anisotropic half-Maxwellian provides a better
approximation of the distribution function close to the evaporating surface
was conjectured for the first time in Ref. Crout (1936) and the theoretical
analysis has been more recently extended to include the case of an evaporation
coefficient less than one Bird and Liang (2019).

According to review above, most of the MD studies of liquids evaporating into
vapor agree that the distribution function of evaporated atoms deviates from
the half-Maxwellian, albeit there is no consensus on the fitting function. The
simulation results of liquids evaporating into near vacuum conditions are even
less conclusive in that deviations have not always been found.

In some MD studies, it has been shown that, in the limit of a low-density
vapor, evaporated argon atoms are distributed according to a half-Maxwellian
Zhakhovskii and Anisimov (1997); Zhakhovsky  _et al._ (2019). This conclusion
has been explained by the fact that atoms leave the condensed phase due to
collisions from the edge of the liquid-vapor interface where the binding
energy is negligible, whence the half-Maxwellian shape.

In other MD studies, the velocity distribution function of spontaneously
evaporating atoms has been determined to be nearly a half-Maxwellian at low
temperatures, but a tail fatter than the one of a half-Maxwellian has been
found normal to the liquid-vapor interface at high temperatures Ishiyama  _et
al._ (2004); Heinen and Vrabec (2019).

The validity of standard kinetic boundary conditions at the liquid-vapor
interface has also been assessed based on the Enskog-Vlasov (EV) equation.
This kinetic equation has the capability of describing both the liquid and
vapor phases, including the interface region and has many attractive features.
Compare to MD simulations, its numerical solution requires less computational
effort and, unlike macroscopic approaches, such as diffuse interface models
and sharp interface methods, it can capture the nonequilibrium effects in the
liquid-vapor interface and in the Knudsen layer.

The EV equation has been used to study the evaporation into near vacuum of a
thin liquid film Frezzotti  _et al._ (2005), the steady
evaporation/condensation flow between two planar liquid films kept at
different temperatures Kon  _et al._ (2014), and, more recently, the
evaporation of multi-component substances into vapor and vacuum Kobayashi  _et
al._ (2017); Frezzotti  _et al._ (2018).

Remarkably, in spite of the simplified treatment of pair correlations in the
dense phase, the results obtained by the EV equation closely resemble the ones
provide by MD simulations. More specifically, it is generally concluded that
evaporated atoms are distributed according to a half-Maxwellian with liquid
bulk temperature even though a slight anisotropy between the velocity
components normal and parallel to the liquid-vapor interface has been found at
high evaporation temperatures.

The foregoing literature review shows that uncertainties remain on the
statistical features of atoms spontaneously emitted from the liquid phase.
Most of the studies observed that evaporated atoms are not distributed
according to a half-Maxwellian but none has systematically evaluated the
deviations as a function of the evaporation temperature. Furthermore, the
physical mechanism which leads the velocity distribution function of
evaporated atoms to differ from a half-Maxwellian has not been convincingly
explained.

The aim of the present work is thus twofold. First, to accurately determine
the distribution function of atoms evaporating into near-vacuum conditions by
numerically solving the EV equation. Compared to MD simulations, the numerical
solution of this kinetic equation requires less computational effort and,
therefore, permits one to get results with the required level of accuracy for
the analysis undertaken in this study Frezzotti and Barbante (2017). Second,
to provide a convincing explanation of the deviations from the half-
Maxwellian.

The rest of the paper is organised as follows. Section II outlines the
mathematical formulation of the evaporation of a monatomic single-component
liquid into near-vacuum condition based on the EV equation and outlines the
particle method of solution. Section III contains a brief description of the
computational setup and presents the main results of the paper. In particular,
it is shown that the reduced distribution function of spontaneously evaporated
atoms are well approximated by a drifted anisotropic half-Maxwellian.
Furthermore, it is shown, by numerical evidence and a simple mathematical
model, that deviations from an isotropic half-Maxwellian are due to the atoms’
collisions in the liquid-vapor interface. Section IV summarises and comments
on the main findings of the paper.

##  II Mean-field kinetic theory approach

###  II.1 The Enskog-Vlasov equation

Let us consider a fluid composed of spherical and identical atoms of mass $m$
and diameter $a$ interacting through the Sutherland potential given by a
superposition of a hard sphere potential and an attractive soft potential
tail:

$\phi(\rho)=\begin{cases}+\infty\,,\quad&\rho<a,\\\
\displaystyle{-\phi_{a}\left(\frac{\rho}{a}\right)^{-\gamma}}\,,\quad&\rho\geq
a,\end{cases}$ (1)

where $\rho=||\bm{r}_{1}-\bm{r}||$ is the distance between the atoms at
position $\bm{r}_{1}$ and $\bm{r}$ while the two positive constants $\phi_{a}$
and $\gamma$ define the depth of the potential well and the range of the soft
interaction, respectively. Karkheck et al Karkheck and Stell (1981) derived
the exact evolution equation for the one-particle distribution function of
this system of atoms, $f(\bm{r},\bm{v},t)$, but this equation is of little use
since it involves the two-particle distribution function. In order to obtain a
closed equation, two simplifying assumptions are introduced in the particle
dynamics, namely long-range particle correlations are neglected while short-
range particle correlations are approximated by following the Enskog theory
originally derived for dense gases. By adopting these assumptions, the
following closed equation for the evolution of $f(\bm{r},\bm{v},t)$ is
obtained:

$\frac{\partial f}{\partial
t}+\bm{v}\cdot\nabla_{\bm{r}}f+\frac{\mathcal{F}[n]}{m}\cdot\nabla_{\bm{v}}f=\mathcal{C}_{E}[f],$
(2a) where square brackets denote functional dependence. In Eq. (2a),
$\mathcal{F}(\bm{r},t)$ is the self-consistent force field generated by the
soft attractive tail which reads:
$\mathcal{F}[n]=\int_{||\bm{r}_{1}-\bm{r}||>a}\frac{d\phi_{a}(\rho)}{d\rho}\frac{\bm{r}_{1}-\bm{r}}{||\bm{r}_{1}-\bm{r}||}n(\bm{r_{1}})d\bm{r_{1}},$
(2b) and the hard-sphere collision integral $\mathcal{C}_{E}(f,f)$ is given
by:
$\mathcal{C}_{E}[f]=a^{2}\int(\bm{v}_{r}\cdot\bm{\hat{k}})^{+}d\bm{v}_{1}d^{2}\bm{\hat{k}}\left\\{\chi\left[n\left(\bm{r}+\frac{a}{2}\bm{\hat{k}},t\right)\right]f(\bm{r}+a\bm{\hat{k}},\bm{v}_{1}^{*},t)f(\bm{r},\bm{v}^{*},t)-\right.\\\
\left.\chi\left[n\left(\bm{r}-\frac{a}{2}\bm{\hat{k}},t\right)\right]f(\bm{r}-a\bm{\hat{k}},\bm{v}_{1},t)f(\bm{r},\bm{v},t)\right\\},$
(2c) where $(\cdot)^{+}$ indicates that the surface integral is restricted to
the half-sphere for which $\bm{v}_{r}\cdot\bm{\hat{k}}>0$ and $\chi[n]$ is the
contact value of the pair correlation function in a hard-sphere fluid in
equilibrium with number density $n$.

In the Standard Enskog Theory (SET), $\chi$ is approximated by using the value
of the pair correlation function in a fluid in uniform equilibrium with the
value of the density at the contact point of the two colliding atoms. An
approximate, but accurate expression for $\chi_{\mbox{\tiny SET}}$ can be
obtained from the equation of state of the hard-sphere fluid proposed by
Carnahan and Starling Carnahan and Starling (1969), as:

$\chi_{\mbox{\tiny
SET}}(n)=\frac{1}{nb}\left(\frac{p^{hs}}{nk_{B}T}-1\right)=\frac{1}{2}\frac{2-\eta}{(1-\eta)^{3}},\qquad
b=\frac{2\pi a^{3}}{3},\qquad\eta=\frac{\pi a^{3}n}{6}.$ (3)

where $p^{hs}$ is the pressure of a system of hard spheres and $k_{B}$ is the
Boltzmann constant. In the present work, following the Fischer-Methfessel
approach Fischer and Methfessel (1980), the actual value of the density at the
contact point is replaced with the value of the density field averaged over a
spherical volume of radius $a$, namely:

$\chi[n](\bm{r},\bm{r}-a\bm{\hat{k}})=\chi_{\mbox{\tiny
SET}}\left(\overline{n}\left(\bm{r}-a\frac{\bm{\hat{k}}}{2}\right)\right),$
(4a) where $\overline{n}(\bm{r},t)=\frac{3}{4\pi
a^{3}}\int_{\mathbb{R}^{3}}n(\bm{r}_{*},t)w(\bm{r},\bm{r}_{*})\,d\bm{r}_{*},\hskip
28.45274ptw(\bm{r},\bm{r}_{*})=\left\\{\begin{array}[]{cc}1,&\qquad\|\bm{r}_{*}-\bm{r}\|<a,\\\
0,&\qquad\|\bm{r}_{*}-\bm{r}\|>a.\end{array}\right.$ (4b)

The kinetic equation (2) is usually referred to as the Enskog-Vlasov (EV)
equation De Sobrino (1967); Grmela (1971); Karkheck and Stell (1981);
Frezzotti  _et al._ (2005); Benilov and Benilov (2019). This kinetic equation
has been applied to investigate a wide range of two-phase flows, including
condensation/evaporation processes Frezzotti  _et al._ (2005); Kon  _et al._
(2014); Frezzotti  _et al._ (2019), liquid menisci in nano-channels Barbante
_et al._ (2015), and its mathematical properties have been extensively studied
in connection with the liquid-vapor transition Takata  _et al._ (2018);
Benilov and Benilov (2018).

###  II.2 Particle method of solution

In this work, the EV equation is solved numerically by an extension of the
original Direct Simulation Monte-Carlo (DSMC) scheme to dense fluids Frezzotti
(1997). A thorough description of the numerical scheme and the analysis of its
computational complexity is given in Ref. Frezzotti  _et al._ (2019).

For EV simulations, the main framework of DSMC scheme used to solve the
Boltzmann equation is preserved, with modifications occurring in the collision
algorithm due to the nonlocal structure of the Enskog collision operator. The
distribution function is represented by $N$ computational particles:

$f(\bm{r},\bm{v},t)=\sum_{i=1}^{N}\delta{\left(\bm{r}-\bm{r}_{i}(t)\right)}\delta(\bm{v}-\bm{v}_{i}(t)),$
(5)

where $\bm{r}_{i}$ and $\bm{v}_{i}$ are the position and the velocity of the
$i$th particle at time $t$, respectively.

The distribution function is updated by a fractional-step method based on the
time-splitting of the evolution operator in two sub-steps, namely free
streaming and collision. In the first stage, the collisions between particles
are neglected and the distribution function is advanced from $t$ to $t+\Delta
t$ by solving the equation:

$\frac{\partial f}{\partial
t}+\bm{\bm{v}}\cdot\nabla_{\bm{r}}f+\frac{\mathcal{F}[n]}{m}\cdot\nabla_{\bm{\bm{v}}}f=0,$
(6)

which translates into updating the positions and velocities of the
computational particles according to:

$\displaystyle\bm{r}_{i}(t+\Delta t)$
$\displaystyle=\bm{r}_{i}(t)+\bm{v}_{i}\Delta
t+\frac{\mathcal{F}[n(t)]}{m}\frac{(\Delta t)^{2}}{2},$ (7a)
$\displaystyle\bm{v}_{i}(t+\Delta t)$
$\displaystyle=\bm{v}_{i}(t)+\frac{\mathcal{F}[n(t)]}{m}\Delta t.$ (7b)

In the second stage, the short range hard-sphere interactions are considered
and the updating rule is given by:

$f(\bm{r},\bm{v},t+\Delta t)=\tilde{f}(\bm{r},\bm{v},t+\Delta
t)+\mathcal{C}_{E}[\tilde{f}]\Delta t.$ (8)

During this stage, particles’ positions $\bm{x}_{i}$ are unchanged while their
velocities $\bm{v}_{i}$ are modified according to stochastic rules which
essentially correspond to the Monte Carlo evaluation of the collision integral
given by Eq. (2c).

The macroscopic quantities are obtained by time averaging the particles’
microscopic states. Note that for steady flows simulations, as the ones
considered below, the averaging time can be long enough to obtain accurate
results without the need of using a large number of computational particles.

##  III Evaporation of monatomic liquid into near vacuum

Figure 1: Schematic of the computational setup. The simulation domain is a
finite symmetric interval with perfectly absorbing boundary conditions at
$z=\pm 60a$. The initial condition consists of a liquid slab in equilibrium
with its vapor placed in the middle of the simulation domain.

###  III.1 Computational setup

The schematic of the computational setup is shown in Fig. 1. The simulation
domain is a finite symmetric interval $z\in[-L,L]$ with $L=60a$ and
homogeneous along the $x$ and $y$ axis. Perfectly absorbing boundary
conditions are assumed at $z=\pm L$, namely computational particles reaching
the boundaries are removed from the simulation domain. Initially, a liquid
slab in equilibrium with its vapor is considered in the middle of the
simulation domain. This initial state has been obtained by placing a liquid
slab in $z\in[-L_{s},L_{s}]$, with $L_{s}=50a$, gas elsewhere, and using
periodic boundary conditions at $z=\pm L$. Afterwards, the system has been let
to evolve until the equilibrium steady state is reached. The distance between
the liquid-vapor interfaces and the absorbing surfaces is about $10a$, namely
large enough not to influence the interface but, at the same time,
sufficiently small to minimise the backscattered flow.

The system is clearly symmetric with respect to the $z$ axis and, therefore,
in principle, the computational setup might be simplified by considering only
half of the domain. However, this would lead to the need of imposing a
specular boundary condition at the center of the liquid slab which, among
other things, would require a tricky treatment of the mean force field.
Furthermore, results in the two halves of the domain can be superimposed and,
therefore, the larger computational effort of simulating all the system is
used to reduce the statistical noise.

The simulation domain is divided into 2400 cells with size $\Delta z=a/20$ and
the time step is $\Delta t=2\times 10^{-4}\,a/(RT_{0})^{1/2}$, where
$R=k_{B}/m$ is the specific gas constant and $T_{0}$ is the reference
temperature. The number of computational particles is set to $1.2\times
10^{6}$ and made equal to the number of real atoms by a proper choice of the
cross section normal to the non-homogeneous direction $z$. The interaction
parameters were chosen to be $\phi_{a}/(k_{B}T_{0})=1$ and $\gamma=6$ so as to
match the same far field behavior as the 12-6 Lennard-Jones potential
Hirschfelder  _et al._ (1964). The study of evaporation of a liquid slab is
carried out at the following temperature values
$T_{\ell}/T_{c}=\\{0.53,0.596,0.663,0.729\\}$, where $T_{c}$ is the critical
temperature $T_{c}/T_{0}=0.754632$ as a consequence of the interaction
parameters chosen Frezzotti  _et al._ (2005). In this temperature range, the
vapor phase contains a number of particles sufficient to limit the statistical
noise of the results (which are obtained by averaging particles’ properties)
but, at the same time, is dilute enough to behave as an ideal gas. The
Andersen thermostat Andersen (1980) is applied in the central part of the
liquid slab, $10a$ wide, to prevent the progressive cooling of the system and
keep its temperature to the constant value of $T_{\ell}$. This thermostat is
chosen for its simplicity and computational efficiency.

Note that, in principle, during the evaporation into vacuum, the liquid slab
is consumed and the liquid-vapor interfaces slowly recede from the absorbing
surfaces. The flow is thus an unsteady process and the possibility to evaluate
the macroscopic quantities by time-averaging is jeopardized. In order to
circumvent this difficulty, in the present work the evaporation process is
studied in a frame of reference fixed relative to the liquid-vapor interface
using the following procedure. The simulation advances until the number of
particles reaching the absorbing surfaces equals or exceeds $2N_{\ell}$, where
$N_{\ell}$ is the average number of particles per cell in the liquid bulk at
temperature $T_{\ell}$. When this happens, the simulation stops. All particles
are then moved by $\Delta z$ towards the closest absorbing surface and the
empty gap, which forms in the centre of the simulation domain, is filled with
$2N_{\ell}$ particles sampled from a Maxwellian with temperature $T_{\ell}$.
Afterwards, the simulation is restarted. Note that the interface movement is
negligibly small during two successive applications of this procedure. In this
steady evaporation framework, simulations are run for the time duration
$2000\,a/(RT_{0})^{1/2}$ (equivalent to $10^{7}$ iterations).

###  III.2 Simulation results

####  III.2.1 Macroscopic quantities

(a) $T_{\ell}/T_{c}=0.530$

(b) $T_{\ell}/T_{c}=0.729$

Figure 2: Dimensionless density, $na^{3}$, mean velocity in the $z$-direction,
$V_{z}/(RT_{0})^{1/2}$, normal, $T_{\perp}/T_{c}$, parallel,
$T_{\parallel}/T_{c}$, and total, $T/T_{c}$, temperatures for (a) the lowest
and (b) the highest liquid bulk temperatures considered in the simulation
campaign.

A broad picture of a liquid slab evaporating into near-vacuum is provided by
observing the behaviour of the macroscopic quantities across the the liquid
and vapor phases shown in Fig. 2 for the lowest ($T_{\ell}/T_{c}=0.53$) and
the highest ($T_{\ell}/T_{c}=0.729$) temperatures considered. Due to the
symmetry of the system, only half of the simulation domain is shown.

The minimum density in the liquid phase is found in the central part of the
domain where the liquid is thermostated, while the maximum is reached close to
the liquid-vapor interface, with an almost linear profile in between them.
This behavior is due to the evaporation cooling of the liquid slab, which
causes a temperature decrease in the liquid region which is not thermostated.
After the drop in the liquid-vapor interface, the temperature in the vapor
phase levels off but still exhibits a small gradient. This clearly indicates
that some collisions between atoms occur in the gap region between the edge of
the interface at the vapor side and the absorbing surface.

The normal and parallel temperatures, which are defined based on the velocity
components normal and parallel to the liquid-vapor interface, superimpose up
to a point located at about half of the interface width. Their separation
indicates that the local thermodynamic quasi-equilibrium condition breaks
down, marking the beginning of a transition layer which extends a few
molecular diameters into the low density region. In the following, we refer to
the point at which the separation of temperatures occurs as the ‘separation
point’, denoted $z_{s}$, and the corresponding value of the temperature as the
‘separation temperature’, denoted $T_{s}$. After the separation point the
parallel temperature profile follows a similar slope as the one in the liquid
bulk, while the normal temperature exhibits a sharp drop coupled with an
increase in the mean velocity component normal to the liquid-vapor interface.

Note that the normal temperature and the mean velocity in the $z$-component
show gradients in a thin layer close to the boundary of the computational
domain, especially for the highest temperature case shown in Fig. 2(b). This
behaviour is not a numerical artifact but it is due to the the discontinuity
of the density field at the absorbing surface, which results in a mean force
field directed towards the liquid acting on the atoms in this boundary region.

####  III.2.2 Velocity distribution function of evaporated atoms

| Parameters | Undrifted anisotropic | Drifted isotropic | Drifted anisotropic
---|---|---|---|---
$T_{\ell}/T_{c}=0.530$ | $\xi/(RT_{0})^{1/2}$ | — | $0.026366\pm 0.00388$ | $0.03427\pm 0.01364$
| $\theta_{\perp}/T_{c}$ | $0.5278\pm 0.00254$ | $0.51238\pm 0.00176$ | $0.5072\pm 0.00852$
| $\theta_{\parallel}/T_{c}$ | $0.51238\pm 0.00176$ | $0.51238\pm 0.00176$ | $0.51238\pm 0.00176$
| $\chi^{2}$ | $0.03016$ | $0.02866$ | $0.02856$
$T_{\ell}/T_{c}=0.729$ | $\xi/(RT_{0})^{1/2}$ | — | $0.10903\pm 0.00246$ | $0.1628\pm 0.0056$
| $\theta_{\perp}/T_{c}$ | $0.66987\pm 0.00348$ | $0.59937\pm 0.00087$ | $0.55924\pm 0.00397$
| $\theta_{\parallel}/T_{c}$ | $0.59937\pm 0.00087$ | $0.59937\pm 0.00087$ | $0.59937\pm 0.00087$
| $\chi^{2}$ | $0.03911$ | $0.01163$ | $0.00688$
Table 1: Estimates of the fitting parameters, $\xi,\,\theta_{\perp}$,
$\theta_{\parallel}$ alongside the asymptotic standard error, and the residual
sum of squares, $\chi^{2}$, of the three considered fitting velocity
distribution functions for the lowest and the highest liquid bulk temperatures
considered in the simulation campaign.

The statistical features of spontaneously evaporated atoms have not been
systematically assessed until now. According to some studies their velocity
distribution function can be approximated by an anisotropic half-Maxwellian
Cheng  _et al._ (2011); Ishiyama  _et al._ (2004); Frezzotti  _et al._ (2018)
while others pointed out the presence of a velocity drift Meland  _et al._
(2004). Here, the velocity distribution function of evaporated atoms is
evaluated based on the particles collected at the absorbing surfaces. As
specified in Sec. III.1, these surfaces are placed at the edge of the liquid-
vapor interface where kinetic boundary conditions can be formulated. Three
different Maxwellian-like functional forms are tested for fitting:

Undrifted Anisotropic: $\displaystyle
f(v_{\parallel},v_{\perp})=\frac{C_{1}}{(2\pi
R)^{3/2}\theta_{\perp}\theta_{\parallel}^{1/2}}\exp{\left[-\frac{v_{\perp}^{2}}{2R\theta_{\perp}}-\frac{v_{\parallel}^{2}}{2R\theta_{\parallel}}\right]},\hskip
5.69046ptv_{\perp}>0,$ (9a) Drifted Isotropic: $\displaystyle
f(v_{\parallel},v_{\perp})=\frac{C_{2}}{(2\pi
R\theta_{\parallel})^{3/2}}\exp{\left[-\frac{\left(v_{\perp}-\xi\right)^{2}}{2R\theta_{\parallel}}-\frac{v_{\parallel}^{2}}{2R\theta_{\parallel}}\right]},\hskip
5.69046ptv_{\perp}>0,$ (9b) Drifted Anisotropic: $\displaystyle
f(v_{\parallel},v_{\perp})=\frac{C_{3}}{(2\pi
R)^{3/2}\theta_{\perp}\theta_{\parallel}^{1/2}}\exp{\left[-\frac{\left(v_{\perp}-\xi\right)^{2}}{2R\theta_{\perp}}-\frac{v_{\parallel}^{2}}{2R\theta_{\parallel}}\right]},\hskip
5.69046ptv_{\perp}>0,$ (9c)

where $C_{i}$ are constants which make the velocity distribution functions
normalised to unity, $v_{\perp}$, $v_{\parallel}$ are the velocity components
normal and parallel to the liquid-vapor interface, and $\xi$,
$\theta_{\perp}$, $\theta_{\parallel}$ are the fitting free parameters. Note
that $\xi$ cannot be identified with the mean velocity in the normal
$z$-direction since the velocity distribution function is defined only for
$v_{\perp}>0$. Likewise, $\theta_{\perp}$ and $\theta_{\parallel}$ cannot be
identified with the normal and parallel temperatures.

The values of the fitting parameters are listed in Table 1 for the lowest and
highest values of the liquid bulk temperatures considered in the simulation
campaign. The sum of squares of residuals, $\chi^{2}$, which is an indicator
of the goodness of the fit, is also reported. According to $\chi^{2}$, the
drifted anisotropic half-Maxwellian provides the best approximation,
especially for the highest temperature $T_{\ell}/T_{c}=0.729$.

(a) Normal molecular flux

(b) Residuals

(c) Q-Q plot

Figure 3: (a) Reduced normal molecular fluxes normalized to unity of
evaporated molecules, (b) residual distribution, and (c) corresponding Q-Q
plots of the residuals for the lowest liquid bulk temperature considered in
the simulation campaign ($T_{\ell}/T_{c}=0.53$). Dashed-dotted line: Fitted
undrifted anisotropic half-Maxwellian; Dashed: Fitted isotropic drifted half-
Maxwellian; Solid line: Fitted drifted anisotropic half-Maxwellian. Inset: The
tail of the distributions in logarithmic scale, with x axis tick values
matched to the ticks of the plot.

(a) Normal molecular flux

(b) Residuals

(c) Q-Q plot

Figure 4: (a) Reduced normal molecular fluxes normalized to unity of
evaporated molecules, (b) residual distribution, and (c) corresponding Q-Q
plots of the residuals for the highest liquid bulk temperature considered in
the simulation campaign ($T_{\ell}/T_{c}=0.729$). Dashed-dotted line: Fitted
undrifted anisotropic half-Maxwellian; Dashed: Fitted isotropic drifted half-
Maxwellian; Solid line: Fitted drifted anisotropic half-Maxwellian. Inset: The
tail of the distributions in logarithmic scale, with x axis tick values
matched to the ticks of the plot.

The direct comparison between the fitting curves and the data is shown in
Figs. 3 and 4 alongside their residuals and residuals Q-Q scatter plots which
provide a more quantitative assessment of the goodness of fit. More
specifically, the Q-Q scatter plot compares the quantiles of the residuals
distribution functions against one another. The alignment along the bisector
of first quadrant indicates that residuals are normally distributed as it
should be for the ideal fitting.

In the lowest temperature case shown in Fig. 3, all the proposed distributions
superimpose almost perfectly and their Q-Q scatter plots indicate a near
normal distribution of residuals. By contrast, in the highest temperature case
shown in Figs. 4, deviations can be observed between the different
distributions, with the drifted anisotropic half-Maxwellian clearly providing
the best fit. The direct inspection of residuals distributions, Fig. 3b and
Fig. 4b, shows that the undrifted anisotropic half-Maxwellian and the drifted
isotropic half-Maxwellian do not provide a good fit in the peak and tail
regions, while residuals of the drifted anisotropic half-Maxwellian are
randomly distributed except for the tail region where a weak pattern is
visible.

None of the fitting functions perfectly matches the data but the analysis
above indicates that the drifted anisotropic half-Maxwellian provides the best
approximation, albeit at the cost of an extra fitting parameter. In the
considered range of temperatures, using a distribution function with only one
fitting parameter leads to errors in the mean velocity and temperature within
one percent but errors rapidly increases if one considers higher-order moments
which are more sensitive to the accuracy of the fitting in the tail region.
Accordingly, for the remainder of the paper, we assume that the velocity
distribution function of evaporated atoms is an anisotropic drifted half-
Maxwellian.

(a) $T_{\ell}/T_{c}=0.530$

(b) $T_{\ell}/T_{c}=0.596$

(c) $T_{\ell}/T_{c}=0.663$

(d) $T_{\ell}/T_{c}=0.729$

Figure 5: Reduced velocity distribution function and molecular flux of
evaporated molecules, parallel and normal to the liquid-vapor interface,
respectively, at different liquid bulk temperatures. Coloured histograms are
the numerical results of the EV equation normalized to unity; solid and dashed
lines are their best fits based on a drifted anisotropic half-Maxwellian with
parameters $\xi$, $\theta_{\perp}$, and $\theta_{\parallel}$.

The reduced distribution function of the parallel velocity components $v_{x}$
and $v_{y}$, denoted $f_{x,y}$, and the reduced flux of the normal-velocity
component $v_{z}$, denoted $v_{z}f_{z}$, are shown in Fig. 5 while the values
of $\xi$, $\theta_{\perp}$, and $\theta_{\parallel}$ are summarised in Table 2
and plotted in Fig. 6. It can be observed from Fig. 6(b) that as the liquid
bulk temperature increases, the difference between the normal and parallel
fitting parameters becomes larger. The drift velocity also grows monotonically
with the bulk temperature, as it can be seen in Fig. 6(a).

It is worth stressing that, the compressibility factor of the vapor is only
slightly less than one in the entire range of evaporation temperatures
considered in the simulation campaign (see Table 2). Accordingly, the velocity
drift and temperature anisotropy pointed out in this study are expected to
show up even in vapors whose behaviour is only slightly non-ideal. Note that,
in the previous molecular dynamics studies Cheng  _et al._ (2011); Ishiyama
_et al._ (2004); Meland  _et al._ (2004) the deviations from the undrifted
isotropic half-Maxwellian were generically attributed to the collisions of
atoms in the liquid-vapor interface. We will examine this argument in more
depth in the next section.

$T_{\ell}/T_{c}$ |  $T_{s}/T_{c}$ |  $\xi/(RT_{0})^{1/2}$ |  $\theta_{\perp}/T_{c}$ |  $\theta_{\parallel}/T_{c}$ |  $n_{v}a^{3}$ |  $Z$
---|---|---|---|---|---|---
0.530 |  $0.518$ |  0.034 |  $0.507$ |  $0.512$ |  0.0008 |  0.990
0.596 |  $0.561$ |  0.081 |  $0.526$ |  $0.547$ |  0.0018 |  0.981
0.663 |  $0.597$ |  0.133 |  $0.541$ |  $0.579$ |  0.0028 |  0.971
0.729 |  $0.628$ |  0.163 |  $0.559$ |  $0.599$ |  0.0040 |  0.961
Table 2: Separation temperature, $T_{s}/T_{c}$, velocity drift,
$\xi/(RT_{0})^{1/2}$, normal and parallel temperature, $\theta_{\perp}/T_{c}$
and $\theta_{\parallel}/T_{c}$, of the velocity distribution function of
spontaneously evaporated atoms, number density in the vapor region,
$n_{v}a^{3}$, and vapor compressibility index, $Z=p_{v}/(n_{v}k_{B}T_{\ell})$,
being $p_{v}$ the pressure in the vapor phase and $k_{B}$ the Boltzmann
constant, as a function of the liquid bulk temperatures, $T_{\ell}/T_{c}$.

(a) Velocity drift.

(b) Parallel and normal temperatures.

Figure 6: Parameters of the velocity distribution function of spontaneously
evaporating atoms as a function of the liquid bulk temperature.

###  III.3 Velocity drift and temperature anisotropy

(a) $T_{\ell}/T_{c}=0.530$

(b) $T_{\ell}/T_{c}=0.729$

Figure 7: Histogram of the last collision cell of evaporated atoms for the
lowest and highest liquid bulk temperatures considered in the simulation
campaign. The peak of the distribution is close to the separation point
$z_{s}$ which is marked by the vertical black line.

####  III.3.1 Numerical study

In order to track the origin of the statistical properties of evaporated
atoms, their trajectories are traced backwards, to the points where they
collide for the last time. The spontaneously evaporating atoms can thus be
identified as atoms that, coming from different regions of the domain, reach
the absorbing surfaces in free-molecular motion under the action of the
conservative mean force field. These atoms are then divided into two groups,
i.e. those whose last collision takes place before and after the separation
point. Finally, the distribution functions of these two groups of atoms and
the corresponding fitting parameters are evaluated. Figure 7 shows the
histograms of the last collision cell of evaporated atoms for
$T_{\ell}/T_{c}=0.53$ and $T_{\ell}/T_{c}=0.729$. As reported in Table 3, most
of the atoms collected at the absorbing surfaces had their last collision
after the separation point, i.e. they represent more than $80\%$ of the total
evaporated atoms. Herein, the total number of evaporated atoms is denoted
$N_{e}$. Remarkably, the velocity drift and temperature anisotropy is much
larger for atoms whose last collision was before $z_{s}$, where the local
velocity distribution functions are undrifted and isotropic Maxwellians. It is
thus reasonable to focus on this group of atoms to shed light on the mechanism
that leads to the velocity drift and temperature anisotropy.

| Before $z_{s}$ | After $z_{s}$
---|---|---
$T_{\ell}/T_{c}$ | $\xi/(RT_{0})^{1/2}$ | $\theta_{\perp}/T_{c}$ | $\theta_{\parallel}/T_{c}$ | $N_{e,z<z_{s}}/N_{e}$ | $\xi/(RT_{0})^{1/2}$ | $\theta_{\perp}/T_{c}$ | $\theta_{\parallel}/T_{c}$ | $N_{e,z>z_{s}}/N_{e}$
0.530 | 0.154 | 0.455 | 0.495 | 17.8% | 0.036 | 0.512 | 0.518 | 82.2%
0.729 | 0.344 | 0.484 | 0.567 | 11.2% | 0.161 | 0.559 | 0.602 | 88.8%
Table 3: Parameters of the velocity distribution function of spontaneously
evaporating atoms before and after the separation temperature for the lowest
and largest liquid bulk temperatures considered.

(a) $T_{\ell}/T_{c}=0.53$

(b) $T_{\ell}/T_{c}=0.729$

Figure 8: Isocontours of the distribution functions of evaporated atoms
(dashed lines) and ‘potentially’ evaporating atoms (solid lines). Atoms
originate from locations before the separation point.

In Fig. 8 we plot the isocontours of the velocity distribution functions of
the evaporated atoms originating before the separation point (group A) and of
the ‘potentially’ evaporating atoms (group B). This latter group comprises all
the atoms coming from the same locations of the ones of group A under the
assumption that they can reach the absorbing surface without suffering any
backscattering due to collisions in the interface region. Both velocity
distribution functions are normalized to unity. The velocity distribution
function of the group B is an undrifted isotropic half-Maxwellian with the
weighted average temperature of the region where they originate. This is not
unexpected, since, as also proved in Section III.3.2, an undrifted Maxwellian
velocity distribution function of atoms moving in a conservative force field
remains Maxwellian with the same temperature. The comparison clearly shows
that the evaporated atoms have predominantly higher normal velocities and
lower parallel velocities. This result can be understood as the result of the
interplay between normal and parallel velocities due to collisions in the
liquid-vapor interface. Indeed, the higher is the normal-velocity component,
the less likely an atom suffers a collision in the liquid-vapor interface
because the lower is the time spent in that region. However, at the same time,
the larger is the speed, the greater becomes the probability of collisions
(see also Eq. (12b)). This interplay also explains why atoms coming from the
region before the separation point mostly contribute to the velocity drift and
temperature anisotropy. Indeed, these atoms must travel a larger distance
before reaching the absorbing surface and are thus more likely to be back-
scattered by collisions. Likewise, the liquid-vapor interface gets wider as
the liquid bulk temperature increases and, therefore, the effect of
preferential evaporation of atoms with large normal and lower parallel
velocity components is enhanced.

####  III.3.2 An instructive model

$T_{\ell}/T_{c}$ | $\bar{d}/a$ | $\bar{n}a^{3}$ | $\bar{T}/T_{c}$ | $\Delta\mathcal{U}/(k_{B}T_{s})$
---|---|---|---|---
0.530 | 9.8 | 0.00157 | 0.44 | 0.01686
0.729 | 12.4 | 0.00554 | 0.497 | 0.04936
Table 4: Numerical values of the parameters which enter in the model for the
lowest and highest liquid bulk temperatures considered in the simulation
campaign

(a) EV solution, $T_{\ell}/T_{c}=0.53$

(b) Model prediction, $T_{\ell}/T_{c}=0.53$

(c) EV solution, $T_{\ell}/T_{c}=0.729$

(d) Model prediction, $T_{\ell}/T_{c}=0.729$

Figure 9: Isocontours of distribution functions of evaporated particles for
two evaporation temperatures.

The results presented in the previous section suggest that the velocity drift
and the temperature anisotropy of spontaneously evaporating atoms is a
consequence of the collisions in the liquid-vapor interface which filter out
atoms with a lower normal-velocity component.

In order to support this conclusion, here we approximately evaluate the
distribution function of atoms which, being initially at the separation point,
cross the liquid-vapor interface, and are collected at the absorbing surface
placed at a distance $\bar{d}$. For simplicity, it is assumed that atoms move
in a region of constant density, $\bar{n}$, and temperature, $\bar{T}$. At the
separation point, atoms are supposed to be distributed according to an
isotropic Maxwellian with zero drift and temperature $T_{s}$:

$f_{s}\left(v_{\parallel},v_{\perp}\right)=\frac{n_{s}}{\left(2\pi
RT_{s}\right)^{3/2}}\exp\left(-\frac{v_{\parallel}^{2}+v_{\perp}^{2}}{2RT_{s}}\right).$
(10)

Only a fraction of these atoms can reach the absorbing surface, namely the
ones which (i) have a speed sufficiently large to overcome the potential
barrier of the mean force field which pushes atoms towards the liquid phase
and (ii) are not backscattered due to the collisions with other atoms. Note
that:

  * (i)

The conservation of mechanical energy implies:

$\frac{1}{2}mv_{s}^{2}+{\mathcal{U}}_{s}=\frac{1}{2}mv_{e}^{2}+{\mathcal{U}}_{e}\,\Longrightarrow\,v_{\perp,s}^{2}=v_{\perp,e}^{2}+v_{\perp,min}^{2},\hskip
14.22636ptv_{\perp,min}=\left(2\frac{\Delta\mathcal{U}}{m}\right)^{1/2},$ (11)

where $(v_{s},{\mathcal{U}}_{s})$ and $(v_{e},{\mathcal{U}}_{e})$ are the
atom’s velocity and the potential energy of the mean force field at the
separation point and at the absorbing surface, respectively, and
$\Delta\mathcal{U}=\mathcal{U}_{e}-\mathcal{U}_{s}$. In Eq. (11), it has been
used that $v_{\parallel,s}=v_{\parallel,e}$ since the mean force field acts
along the $z$-direction. It is plain that only atoms whose $z$-component of
the velocity is larger than $v_{\perp,min}$ can reach the absorbing surface,
all the others turn back, under the action of the force field.

  * (ii)

The number of collisions experienced by an atom in the time interval $\Delta
t$ can be assumed to follow a Poisson distribution:

$Pr(\mbox{$N$ Collisions in the time interval $\Delta
t$})=\frac{(\bar{\nu}\Delta t)^{N}e^{-\bar{\nu}\Delta t}}{N!},$ (12a) where
the collision rate $\bar{\nu}$ can be estimated as:
$\bar{\nu}(v)=\chi(\bar{n})\frac{\pi
a^{2}\bar{n}}{\sqrt{\pi}\beta}\left[e^{-\beta^{2}v^{2}}+\left(2\beta
v+\frac{1}{\beta v}\right)\frac{\pi}{2}\operatorname{erf}(\beta v)\right],$
(12b) where $\beta^{2}=1/(2R\bar{T})$ Kennard (1938). Accordingly, the
probability that an atom reaches the absorbing surface without suffering any
collision is: $Pr(\mbox{$0$ Collisions in the time interval $\Delta
t$})=e^{-\bar{\nu}(v)\Delta t}=e^{-\bar{\nu}(v)\frac{\bar{d}}{v_{\perp}}},$
(12c)

where it has been used that the time needed to cross the liquid-vapor
interface is $\Delta t=\bar{d}/v_{\perp}$.

By adopting cylindrical coordinates and using (i)-(ii), the net mass flux of
atoms across the liquid-vapor interface can be written as:

$\int_{0}^{+\infty}dv_{\parallel}\,v_{\parallel}\int_{v_{\perp,min}}^{+\infty}dv_{\perp}\,f_{s}\left(v_{\parallel},v_{\perp}\right)v_{\perp}=\int_{0}^{+\infty}d\tilde{v}_{\parallel}\,\tilde{v}_{\parallel}\int_{0}^{+\infty}d\tilde{v}_{\perp}f_{e}\left(\tilde{v}_{\parallel},\tilde{v}_{\perp}\right)\tilde{v}_{\perp},$
(13)

where $f_{e}$ is given by Eq. (10) and $f_{e}$ is the unknown distribution
function of evaporated atoms, i.e. atoms collected at the absorbing surface.
Note that in Eq. (13), the integral over the azimuthal angle cancels out due
to the symmetry of the system. The first integral on the right hand side of
Eq. (13) can be simplified by making the change of variables
$\left(v_{\perp},v_{\parallel}\right)\rightarrow\left((\tilde{v}_{\perp}^{2}+\tilde{v}_{\perp,min}^{2})^{1/2},\tilde{v}_{\parallel}\right)$
suggested by the conservation of mechanical energy, Eq. (11). The distribution
function of atoms collected at the absorbing surface can thus be readily
obtained:

$f_{e}\left(\tilde{v}_{\parallel},\tilde{v}_{\perp}\right)=\frac{n_{s}}{\left(2\pi
RT_{s}\right)^{3/2}}\exp\left(-\frac{\tilde{v}_{\perp,min}^{2}}{2RT_{s}}\right)\\\
\exp\left(-\frac{\tilde{v}_{\parallel}^{2}+\tilde{v}_{\perp}^{2}}{2RT_{s}}\right)\exp\left[-\frac{\bar{\nu}\left((\tilde{v}^{2}+v^{2}_{\perp,min})^{1/2}\right)\bar{d}}{\left(\tilde{v}_{\perp}^{2}+v_{\perp,min}^{2}\right)^{1/2}}\right].$
(14)

Note that if one disregards collisions, $\bar{\nu}=0$, Eq. (14) simplifies to
a half-Maxwellian with the parameter density given by $n_{s}$ reduced by the
Boltzmann factor $\exp{[-\Delta\mathcal{U}/(k_{B}T_{s})]}$.

The dimensionless values of the width of the liquid-vapor interface,
$\bar{d}/a$, the mean density, $\bar{n}a^{3}$, and temperature,
$\bar{T}/T_{c}$, of the liquid-vapor interface are estimated from the
simulation results presented in the previous sections. The potential jump
$\Delta\mathcal{U}/(k_{B}T_{s})$ is the one that occurs at the edge of a slab
having length $\bar{d}/a$ and density $\bar{n}a^{3}$ placed next to vacuum.
The numerical values of all these parameters are listed in Table 4. Figure 9
shows the comparison between the normalised velocity distribution function of
evaporated atoms obtained by numerically solving the Enskog-Vlasov equation,
and the predictions given by Eq. (14). The very good qualitative agreement
strongly suggests that the collisions in the liquid-vapor interface may be
responsible for deviations from the isotropic half-Maxwellian.

##  IV Conclusions

The Enskog-Vlasov (EV) equation has been used to study the one-dimensional
steady evaporation of a monatomic liquid into near vacuum conditions. The main
aim has been to elucidate the statistical features of atoms spontaneously
emitted by the liquid bulk. This is a key step in formulating kinetic boundary
conditions at the liquid-vapor interface.

The mean-field kinetic theory approach is used in this study since it is by
far less computationally demanding than molecular dynamics simulations and
permits one to get results with the required high level of accuracy.

The velocity distribution function of spontaneously evaporating atoms are
commonly assumed to be an half-Maxwellian at the temperature of the liquid-
vapor interface. By contrast, the main results of this work show that:

  * •

Evaporated atoms are distributed according to a drifted anisotropic half-
Maxwellian. Deviations from the isotropic half-Maxwellian become more
pronounced as the liquid bulk temperature increases.

  * •

The velocity drift and the temperature anisotropy are the results of
collisions in the liquid-vapor interface region which preferentially
backscatter atoms with a lower normal-velocity component.

Note that previous studies based on molecular dynamics simulations have
already pointed out that the distribution function of evaporating atoms
deviate from an isotropic half-Maxwellian but there hasn’t been unanimous
agreement on the functional form, i.e. an anisotropic half-Maxwellian Tsuruta
_et al._ (1999); Ishiyama  _et al._ (2004) or a drifted half-Maxwellian Meland
_et al._ (2004). Furthermore, deviations were noticed only for high
evaporation temperatures and, accordingly, attributed to the non-ideal vapor
behaviour.

The novelty of the present work is thus threefold. First, it establishes that
the drifted anisotropic half-Maxwellian provides the best fitting out of the
considered distributions. Second, it shows that deviations from the half-
Maxwellian occur even when the vapor is only slightly non-ideal, namely in
conditions where, in principle, the usual kinetic-theory treatment of the
evaporation process can still be used Frezzotti  _et al._ (2019). Third, it
brings evidence, through accurate numerical results and a simple model, that
collisions in the interface may explain the velocity drift and the temperature
anisotropy.

It is worth stressing that the results above are expected to be valid even in
presence of net condensation since it is commonly assumed that the evaporation
flux only depends on the state of the fluid in the liquid phase. However, a
more detailed study would be needed to assess to what extent the vapor
dynamics affects the structure of the liquid-vapor interface and, in turn, the
evaporation mass flux.

This work paves the way to a reformulation of boundary conditions usually
adopted at the liquid-vapor interface in kinetic theory studies of
evaporation/condensation processes. In this respect, an interesting research
perspective consists in determining the dependence of the drift-velocity and
temperature anisotropy on the properties of the liquid bulk by means of the
balance equations of mass, momentum, and energy derived from the EV equation.

###### Acknowledgements.

This work has been financially supported in the UK by EPSRC grants
(EP/N016602/1, EP/R007438/1, EP/S029966/1 & EP/P031684/1) and the Leverhulme
Trust.

## References

  * Safarian and Engh (2013) Jafar Safarian and Thorvald A Engh, “Vacuum evaporation of pure metals,” Metall. Mater. Trans. A 44, 747–753 (2013).
  * Li and Xu (2014) Yang Li and Song-Lin Xu, “DSMC simulation of vapor flow in molecular distillation,” Vacuum 110, 40–46 (2014).
  * Sazhin (2017) Sergei S Sazhin, “Modelling of fuel droplet heating and evaporation: recent results and unsolved problems,” Fuel 196, 69–101 (2017).
  * John  _et al._ (2019) Benzi John, Ryan Enright, James E Sprittles, Livio Gibelli, David R Emerson, and Duncan A Lockerby, “Numerical investigation of nanoporous evaporation using direct simulation monte carlo,” Phys. Rev. Fluids 4, 113401 (2019).
  * Chang  _et al._ (2020) Mengzhao Chang, Ziyoung Lee, Sungwook Park, and Suhan Park, “Characteristics of flash boiling and its effects on spray behavior in gasoline direct injection injectors: A review,” Fuel 271, 117600 (2020).
  * Rowlinson and Widom (2013) John Shipley Rowlinson and Benjamin Widom,  _Molecular theory of capillarity_ (Courier Corporation, 2013).
  * Sone (2000) Yoshio Sone, “Kinetic theoretical studies of the half-space problem of evaporation and condensation,” Transport Theor. Stat. 29, 227–260 (2000).
  * Persad and Ward (2016) Aaron H Persad and Charles A Ward, “Expressions for the evaporation and condensation coefficients in the hertz-knudsen relation,” Chem. Rev. 116, 7727–7767 (2016).
  * Fang and Ward (1999) G Fang and C A Ward, “Temperature measured close to the interface of an evaporating liquid,” Phys. Rev. E 59, 417 (1999).
  * Crout (1936) Prescott D Crout, “An application of kinetic theory to the problems of evaporation and sublimation of monatomic gases,” J. Math. Phys. 15, 1–54 (1936).
  * Bird and Liang (2019) Eric Bird and Zhi Liang, “Transport phenomena in the knudsen layer near an evaporating surface,” Phys. Rev. E 100, 043108 (2019).
  * Tsuruta  _et al._ (1999) Takaharu Tsuruta, Hiroyuki Tanaka, and Takashi Masuoka, “Condensation/evaporation coefficient and velocity distributions at liquid–vapor interface,” Int. J. Heat Mass 42, 4107–4116 (1999).
  * Meland  _et al._ (2004) Roar Meland, Aldo Frezzotti, Tor Ytrehus, and Bjørn Hafskjold, “Nonequilibrium molecular-dynamics simulation of net evaporation and net condensation, and evaluation of the gas-kinetic boundary condition at the interphase,” Phys. Fluids 16, 223–243 (2004).
  * Xie  _et al._ (2011) Jian-Fei Xie, Sergei S Sazhin, and Bing-Yang Cao, “Molecular dynamics study of the processes in the vicinity of the n-dodecane vapour/liquid interface,” Phys. Fluids 23, 112104 (2011).
  * Xie  _et al._ (2012) Jian-Fei Xie, Sergei S Sazhin, and Bing-Yang Cao, “Molecular dynamics study of condensation/evaporation and velocity distribution of n-dodecane at liquid-vapour phase equilibria,” J. Therm. Sci. Tech-Jpn. 7, 288–300 (2012).
  * Kon  _et al._ (2014) Misaki Kon, Kazumichi Kobayashi, and Masao Watanabe, “Method of determining kinetic boundary conditions in net evaporation/condensation,” Phys. Fluids 26, 072003 (2014).
  * Kobayashi  _et al._ (2017) Kazumichi Kobayashi, Kiyofumi Sasaki, Misaki Kon, Hiroyuki Fujii, and Masao Watanabe, “Kinetic boundary conditions for vapor–gas binary mixture,” Microfluid. Nanofluid. 21, 53 (2017).
  * Zhakhovskii and Anisimov (1997) VV Zhakhovskii and SI Anisimov, “Molecular-dynamics simulation of evaporation of a liquid,” JETP 84, 734–745 (1997).
  * Ishiyama  _et al._ (2004) Tatsuya Ishiyama, Takeru Yano, and Shigeo Fujikawa, “Molecular dynamics study of kinetic boundary condition at an interface between argon vapor and its condensed phase,” Phys. Fluids 16, 2899–2906 (2004).
  * Frezzotti  _et al._ (2005) Aldo Frezzotti, Livio Gibelli, and Silvia Lorenzani, “Mean field kinetic theory description of evaporation of a fluid into vacuum,” Phys. Fluids 17, 012102 (2005).
  * Frezzotti  _et al._ (2018) A Frezzotti, L Gibelli, Duncan A Lockerby, and James E Sprittles, “Mean-field kinetic theory approach to evaporation of a binary liquid into vacuum,” Phys. Rev. Fluids 3, 054001 (2018).
  * Heinen and Vrabec (2019) Matthias Heinen and Jadran Vrabec, “Evaporation sampled by stationary molecular dynamics simulation,” J. Chem. Phys. 151, 044704 (2019).
  * Zhakhovsky  _et al._ (2019) Vasily V Zhakhovsky, Alexei P Kryukov, Vladimir Yu Levashov, Irina N Shishkova, and Sergey I Anisimov, “Mass and heat transfer between evaporation and condensation surfaces: Atomistic simulation and solution of boltzmann kinetic equation,” PNAS 116, 18209–18217 (2019).
  * Frezzotti and Barbante (2017) Aldo Frezzotti and Paolo Barbante, “Kinetic theory aspects of non-equilibrium liquid-vapor flows,” Mech. Eng. Rev. 4, 16–00540 (2017).
  * Karkheck and Stell (1981) John Karkheck and George Stell, “Kinetic mean-field theories,” J. Chem. Phys. 75, 1475–1487 (1981).
  * Carnahan and Starling (1969) Norman F Carnahan and Kenneth E Starling, “Equation of state for nonattracting rigid spheres,” J. Chem. Phys. 51, 635–636 (1969).
  * Fischer and Methfessel (1980) Johann Fischer and Michael Methfessel, “Born-Green-Yvon approach to the local densities of a fluid at interfaces,” Phys. Rev. A 22, 2836 (1980).
  * De Sobrino (1967) Luis De Sobrino, “On the kinetic theory of a van der waals gas,” Can. J. Phys. 45, 363–385 (1967).
  * Grmela (1971) Miroslav Grmela, “Kinetic equation approach to phase transitions,” J. Stat. Phys. 3, 347–364 (1971).
  * Benilov and Benilov (2019) E S Benilov and M S Benilov, “The enskog–vlasov equation: a kinetic model describing gas, liquid, and solid,” J. Stat. Mech-Theory E 2019, 103205 (2019).
  * Frezzotti  _et al._ (2019) A Frezzotti, P Barbante, and L Gibelli, “Direct simulation monte carlo applications to liquid-vapor flows,” Phys. Fluids 31, 062103 (2019).
  * Barbante  _et al._ (2015) Paolo Barbante, Aldo Frezzotti, and Livio Gibelli, “A kinetic theory description of liquid menisci at the microscale,” Kinet. Relat. Mod. 8, 235–254 (2015).
  * Takata  _et al._ (2018) Shigeru Takata, Takuya Matsumoto, Anna Hirahara, and Masanari Hattori, “Kinetic theory for a simple modeling of a phase transition: Dynamics out of local equilibrium,” Phys. Rev. E 98, 052123 (2018).
  * Benilov and Benilov (2018) E S Benilov and M S Benilov, “Energy conservation and h theorem for the enskog-vlasov equation,” Phys. Rev. E 97, 062115 (2018).
  * Frezzotti (1997) Aldo Frezzotti, “A particle scheme for the numerical solution of the Enskog equation,” Phys. Fluids 9, 1329–1335 (1997).
  * Hirschfelder  _et al._ (1964) Joseph Hirschfelder, R Byron Bird, and Charles F Curtiss,  _Molecular theory of gases and liquids_ (Wiley, 1964).
  * Andersen (1980) Hans C Andersen, “Molecular dynamics simulations at constant pressure and/or temperature,” J. Chem. Phys. 72, 2384–2393 (1980).
  * Cheng  _et al._ (2011) Shengfeng Cheng, Jeremy B Lechman, Steven J Plimpton, and Gary S Grest, “Evaporation of Lennard-Jones fluids,” J. Chem. Phys. 134, 224704 (2011).
  * Kennard (1938) E. H. Kennard,  _Kinetic Theory of Gases_ (McGraw-Hill Book Company, New York and London, 1938).