## 2 Simulation Method
### Liquid-Vapor Phase Equilibria Simulation  
The phase equilibria simulation of argon is performed by placing a condensed liquid film in a simulation box with a prolonged length in the \(z\) direction and the periodic boundary conditions applied in all three spatial directions. The liquid-vapor equilibria condition is achieved by equilibrating the system at the desired temperature while the initially unoccupied free volume is saturated by the evaporated atoms as the simulation evolves in time. The simulation snapshot is shown in Figure 2.  
The positions of the liquid and vapor boundaries were set in terms of the density profile obtained from the production run of the equilibria simulation as shown in Figure 1. The simulation data shows that the density continuously decreases from the liquid phase density to the vapor phase density in the left- and right-side transition regions of the phase equilibria simulation density profile. The transition region density can be well approximated by a hyperbolic tangent function for the left-side interface [29]:  
\[\rho_{le}(z)=\frac{1}{2}(\rho_{v}+\rho_{l})-\frac{1}{2}(\rho_{v}-\rho_{l})tanh \left(\frac{z-z_{0,low}}{\delta_{low}}\right) \tag{3}\]  
and for the right-side interface:  
\[\rho_{ri}(z)=\frac{1}{2}(\rho_{v}+\rho_{l})+\frac{1}{2}(\rho_{v}-\rho_{l}) tanh\left(\frac{z-z_{0,up}}{\delta_{up}}\right) \tag{4}\]  
Here, \(\rho_{l}=1275.8\) kg/m\({}^{3}\) is the simulation liquid phase density, \(\rho_{v}=20.9\) kg/m\({}^{3}\) is the simulation vapor phase density, \(z\) is the absolute coordinate in the z axis, \(z_{0}\) is the center position of the density transition region, and \(\delta\) is the measure of the transition region thickness. Subscripts \(le\) and \(ri\) denote the left- and right-side interfaces, respectively. The transition region thickness \(d\), a distance over which the density at the interphase drops from 90% of liquid density to 10% of liquid density, is related to the measure of interphase thickness by relation \(d=2.1972\delta\)[31]. The values of the fitted transition region parameters are listed in Table 1. With the obtained transition region parameter values, the normalized boundary positions were set according to the previous argon study at 90 K [29]: the liquid boundary position was set to \(z^{*}=-0.9\), while the vapor boundary was set to \(z^{*}=3\). Here, \(z^{*}\) is the normalized coordinate of the density transition region  
\[z^{*}=\frac{z-z_{0}}{d} \tag{5}\]  
Figure 2: Phase equilibria simulation snapshot at a random moment of time. The periodic boundary conditions are applied to \(x\), \(y\) and \(z\) axes.  
The normalized and absolute positions for phase equilibria simulation boundaries (listed in Figure 1) are given in Table 2.  
The atoms that cross the vapor boundaries in an outward direction in equilibria simulation represent the sum of the spontaneously evaporated atoms and the reflected from the interface atoms, as illustrated in Figure 1. Furthermore, the number of argon atoms that crossed the vapor boundaries \(N_{out}(t)\) is a linear function of time, for which the slope coefficient can be found using the linear regression analysis and used to evaluate the outward mass flux [32]:  
\[J_{out}=\frac{m_{Ar}b}{S} \tag{6}\]  
Here, \(m_{Ar}\) is the mass of an argon atom, \(b\) is the slope coefficient of function \(N_{out}(t)\), and \(S=2S_{xy}\) is the surface area of total interface region. The outgoing atom mass flux estimated from the phase equilibria simulation is \(J_{out}=105.9\) gcm\({}^{-2}\)s\({}^{-1}\)).