## 2 Simulation Method
### Net Momenum in Steady-State Vacuum Simulation  
The net momentum problem resulting from evaporated atom transportation needs to be addressed as it is related to the accuracy and the correctness of the proposed method. The transportation of evaporated atoms introduces net momentum changes in the system \(\Delta p_{z}\), which consequently induces the force acting on the liquid phase in the \(z\) direction. This force is proportional to the momentum change \(\Delta p_{z}\) and can be approximately evaluated in the following way:  
\[F_{z}=ma_{z}=\frac{\Delta p_{z}}{\Delta t}=\frac{\Delta m\cdot v_{z}}{\Delta t }\approx\frac{m_{Ar}\Delta N_{sp}\overline{v}_{z}}{\Delta t}=m_{Ar}J_{sp}S_{xy} \overline{v}_{z} \tag{8}\]  
where \(\Delta p_{z}\) is the momentum change over time \(\Delta t\), \(\Delta m\) is the transporter mass over time \(\Delta t\), \(v_{z}\) is the velocity in \(z\) direction, \(\Delta N_{sp}\) is the number of spontaneously evaporated atoms at the right vapor boundary over time \(\Delta t\), \(\overline{v}_{z}\) is the average velocity of evaporated atoms in the \(z\) direction. Furthermore, the atoms evaporated to the left side are reflected by the mirror reflection boundary back to the liquid, thus also inducing the momentum change, which induces the force on the liquid phase equal to \(2m_{Ar}J_{sp}S_{xy}\overline{v}_{z}\) (assuming that all evaporated atoms were reflected by the boundary). However, a fraction of the atoms is reflected by the collisions with the incoming transported atoms; therefore, the total force acting on the liquid phase \(F_{z}\) is expect to be in the interval of \(\left[m_{Ar}J_{sp}S_{xy}\overline{v}_{z};\;3m_{Ar}J_{sp}S_{xy}\overline{v}_{z}\right]\). This force accelerates only the liquid phase, or, in our case, decelerates the evaporated atoms at the right-side interface with rate:  
\[a_{z}=\frac{F_{z}}{m_{liq}}\approx\frac{F_{z}}{Nm_{Ar}} \tag{9}\]  
when the center of mass velocity is subtracted from every atom of the system every simulation timestep to keep the liquid phase velocity at zero. Here, \(N\) is the total number of argon atoms in the system. As a result, the energetic characteristics of the collected ensemble of evaporated atoms in the \(z\) direction should be shifted to the lower energy side. However, the velocity deceleration of the evaporated atoms is negligible over the time periods it takes for the atoms to travel from the liquid surface to the vapor boundary. For example, the evaporated atom with an average velocity \(\overline{v}_{z}=171.53\) m/s is decelerated by negligible the amount of 0.31 m/s (or 0.18% of the initial value) by the time the atom travels the distance from the liquid surface to the right vapor boundary (\(l_{evap}=28.04\) A) with the representative spontaneous evaporation rate \(J_{sp}=5.67\times 10^{27}\times 1/\)m\({}^{2}\)s and the upper estimation of the deceleration rate \(a_{z}=1.91\times 10^{10}\) m/s\({}^{2}\). Consequently, the results show that no noticeable velocity distribution shifts were observed for evaporating atoms compared to the Maxwellian distribution at 90 K. The average velocity of evaporated atoms at 90 K was estimated assuming the Maxwellian distribution in the evaporated atoms [32]:  
\[\overline{v}_{z}=\sqrt{\frac{\pi k_{B}T}{2m_{Ar}}} \tag{10}\]  
Furthermore, the method error is expected to scale with the increasing simulation temperature since the spontaneous evaporation rate \(J_{sp}\) is an increasing function of temperature. Therefore, the deceleration rates and corresponding errors should be evaluated before analyzing the statistics of evaporated atoms at higher temperature simulations.