## 4 Conclusions  
In the present paper, a constant virtual vacuum state simulation method was proposed, in which the virtual vacuum condition is achieved at the right-side interface by transporting the evaporated atoms to the left-side vapor boundary and recentering the systems center of mass at each timestep to adjust for the transportation effects. The proposed method was used to perform virtual vacuum simulations of argon at 90 K in order to study the effects of the chosen temperature control approach in virtual vacuum simulations on the system's temperature profiles, spontaneous evaporation rates and corresponding evaporation coefficient values, and the energetic characteristics of the evaporated atoms. The following conclusions can be made from the results presented in this study:  
* In both the deletion vacuum and the steady-state vacuum simulations, the temperature control approach, in which the thermostat is simply applied only to the liquid phase and the interfaces are simulated in the NVE ensemble, provides non-uniform parabolic temperature profiles across the liquid phase due to the evaporative heat sink at the evaporating surface. Consequently, the temperature at the liquid surface position is obtained lower by several K compared to the targeted simulation temperature of 90 K.
* In steady-state vacuum simulations, the division of the liquid phase into 2 bins, to which the thermostat is applied separately, already flattens the non-uniform temperature profiles and raises the liquid surface temperature; however, the surface temperature reaches saturation value below the targeted temperature with the further reduction of bin width below 9.2 A. Consequently, the division into 2 bins leads to  
Figure 11: Velocity distribution functions of evaporating argon atoms at the right-side vapor boundary in (**a**) \(x\), (**b**) \(y\) and (**c**) \(z\) axes with different thermostat region boundary positions \(z_{t}^{*}\).  
increased spontaneous evaporation rates, with the evaporation coefficient values increasing by 11%, and the spontaneous evaporation rate increases by another 2.5% and reaches saturation value when bin width is reduced to 9.2 A.
* The liquid surface temperature increases linearly with increasing normalized (or absolute) thermostat region boundary position values at the right-side interface; however, the surface temperature reaches the targeted simulation temperature of 90 K only when the thermostat boundary is placed beyond the liquid surface. Furthermore, the spontaneous evaporation rate is an increasing linear function of the liquid surface temperature in a temperature range of 87.55-90 K. In this narrow range, the spontaneous evaporation rate and corresponding evaporation coefficient estimation from MD simulations can vary up to 21%.
* The velocity distributions of evaporating argon atoms in \(x\), \(y\) and \(z\) directions follow the Maxwellian distribution for all simulated cases, and the thermostat has little to no noticeable effect on the energy characteristics and velocity distributions of the evaporating argon atoms, even when the thermostat boundary is placed above the liquid surface.
* The proposed temperature control strategies can be used to obtain the evaporation/condensation coefficient values at desired liquid surface temperatures in practical applications.  
G.S.: Methodology, Investigation, Formal analysis, Writing--Original Draft, Visualization; A.D.: Conceptualization, Supervision, Writing--Review & Editing; E.M.: Investigation, Formal analysis, Writing--Review & Editing; R.N.: Data Curation, Visualization, Writing--Review & Editing. All authors have read and agreed to the published version of the manuscript.  
This research received no external funding.  
Data Availability Statement: Data sharing not applicable.  
The authors declare no conflict of interest.