Molecular Dynamics Simulation of Tolman Length and Interfacial Tension of Symmetric Binary Lennard–Jones Liquid

: The Tolman length and interfacial tension of partially miscible symmetric binary Lennard– Jones (LJ) ﬂuids (A, B) was revealed by performing a large-scale molecular dynamics (MD) simulation with a sufﬁcient interfacial area and cutting distance. A unique phenomenon was observed in symmetric binary LJ ﬂuids, where two surfaces of tension existed on both sides of an equimolar dividing surface. The range of interaction ε AB between the different liquids and the temperature in which the two LJ ﬂuids partially mixed was clariﬁed, and the Tolman length exceeded 3 σ when ε AB was strong at higher temperatures. The results show that as the temperature or ε AB increases, the Tolman length increases and the interfacial tension decreases. This very long Tolman length indicates that one should be very careful when applying the concept of the liquid–liquid interface in the usual continuum approximation to nanoscale droplets and capillary phase separation in nanopores.


Introduction
The surface tension of liquid has been examined widely and is well known. In a gas-liquid interface, two surfaces exist: the surface of tension and equimolar dividing surface. Tolman expressed the decrease in surface tension of a nano-scale droplet by the existence of the two surfaces [1][2][3]. This theory has been improved by many thermodynamic studies [4][5][6][7][8]. The distance between the surface of tension and the equimolar dividing surface is called the Tolman length. Even now, the Tolman length itself and the relationship between the Tolman length and droplet size, evaporation velocity, and other parameters are the focus of many theoretical studies [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. In addition, in studies of flat surfaces of liquids, the surface tension of Lennard-Jones (LJ) liquids has been simulated by several researchers including the author. The relationship between the particle position and local surface tension was derived by Nijmeijer et al. [24]. The minimum system size, liquid thickness, and equilibration time necessary for describing the surface tension were determined by Holcomb et al. [25]. A Tolman length with a cutoff distance of 2.5 σ for various temperatures was simulated by Haye and Bruin [26]. Chen discussed the dependence of surface tension on the surface area. These studies indicated that a surface area of 10 σ × 10 σ is sufficient to provide the correct surface tension [27]. Trokhymchuk and Alejandre discussed the dependence of surface tension on the cutoff distance. They indicated that a cutoff distance of 4.5 σ is insufficient to mimic the full potential, but a length of approximately 5 σ is sufficient [28]. In a previous study, following previous research and by using a cutoff distance of 5 σ, the surface tension and Tolman's length for a LJ liquid were obtained for various temperatures by the author [29]. Recently, Baidakov and Khotienkova have attempted to describe the surface tension of methane-nitrogen mixtures using van der Waals gradient theory and compared this with experimental results [30]. However, there are no examples of research on the Tolman length and interfacial tension for liquid-liquid interfacial tension.
The behavior of capillary phase separation in pores should be compared with the behavior of phase separation in bulk, but surprisingly, there are few examples of studies on the phase separation behavior of LJ liquids in bulk. Only Díaz-Herrera et al. examined the bulk phase diagram and temperature dependence of interfacial tension for various interaction strengths with a cutoff of 3 σ [48,49]. However, they did not examine the surface of tension and equimolar dividing surface in the interfacial region of a binary LJ fluid mixture. In this problem, two surfaces of tension on both sides of an equimolar dividing surface were found by the authors [50]. However, in the simulation, only one condition was considered; i.e., the LJ parameter ε AB = 0.8 ε AA = 0.8 ε BB for symmetric binary LJ liquids "A" and "B", and the reduced temperature T* was 0.8. Moreover, the cutoff distance was 3.5 σ, which was too short to describe the full LJ potential. Therefore, the results of the authors' previous study could not provide useful information on the interfacial physical properties of symmetric binary LJ liquids. In 2020, researchers have focused on the importance of the cutting distance, and gas-liquid equilibrium and critical points and gas-liquid surface tensions have been studied for LJ mixtures of low and high-boiling components with long cutting distances. However, they do not deal with the interfacial tension and Tolman length of the liquid-liquid interface in the partially mixed state [51,52].
Despite recent improvements in computational power and the ability to use more complex simulations than before, there has been no study of the Lennard-Jones fluid interface in two-component systems. Even now, specific values for interfacial tension and the Tolman length are not given for a sufficient interfacial area or sufficient cutting distance.
In this study, molecular dynamics simulations of symmetric binary LJ liquids with a sufficient interfacial area and cutting distance have been performed to provide fundamental insights into their interfacial tension and Tolman length. The purpose of molecular dynamics simulation is not to reproduce reality but to provide the interfacial properties between the LJ liquid and the liquid treated by molecular dynamics simulation as an ideal material, which is required by the thermodynamic equation to explain the various complex phenomena in the nanopores simulated by the molecular dynamics simulation.

Setting
An NVE molecular dynamics simulation of symmetric binary Lennard-Jones LJ fluids was used in this study. We used a self-developed molecular dynamics simulation code. A sandwich-structured liquid film consisting of three layers of LJ (12,6) fluids in a rectangular cell was simulated as shown in Figure 1. A membrane of LJ liquid was placed in the center of the y-direction and sandwiched by membranes of the other LJ liquid from both ends. Truncated and shifted LJ argon was used for the particles. A system was set up comprising 3200 LJ argon "A" particles and 6400 LJ argon "B" particles. The intermolecular interaction parameter for homologous LJ argon "A" and "B" were set as ε AA/ k B = ε BB /k B = 120.0 K, σ AA = σ BB = σ AB = 0.34 nm, where k B is Boltmann's constant. The cutoff distance was 5.0 σ AA , which was sufficiently large to represent particles using the full LJ potential [27]. The masses of LJ argon "A" and "B" were 6.636 × 10 -26 kg [50]. The MD simulations were conducted in a box with dimensions of L x × L y × L z = 14.0 σ AA × 120.0 σ AA × 15.2 σ AA , which was sufficiently large to give a correct value for the surface tension [25,27].
The intermolecular interaction parameter ε AB between the heterogeneous LJ argon was set as ε AB = α × ε AA , where α was between 0.50 and 0.85. The intermolecular interaction parameter σ AB between the heterogeneous LJ argon was set as σ AB = σ AA = σ BB . The leapfrog Verlet method was used to integrate the equations of motion numerically. Each run comprised integration steps with a time increment ∆t of 5 fs. The total time required for each simulation was 1 µs, which was confirmed as a sufficiently long time in a previous study [29]. The reduced temperature T* = k B T/ε AA was set as between 0.7 and 1.0. The lowest temperature T* = 0.7 was set considering that the triple point of pure LJ liquid is T* = 0.685-0.687 [53]. The temperature of the system was controlled by velocity scaling in the usual manner once every 100 steps. Thermal control was applied up to 0.9 µs, after which there was almost no fluctuation in temperature, which was almost equal to the temperature set in the simulation. We did not observe any anomalous increase in the velocity of some particles associated with the adoption of the temperature change algorithm instead of Nosé-Hoover for the NVT ensemble. In addition, we did not adopt the data from the period of the temperature change algorithm (0-0.9 µs) but only analyzed the data from the subsequent period of the NVE ensemble (0.9-1.0 µs). The number of particles, the cutting distance, the interface area, and the total simulation time were very large and long compared to previous studies.  The intermolecular interaction parameter εAB between the heterogeneous LJ argon was set as εAB = α × εAA, where α was between 0.50 and 0.85. The intermolecular interaction parameter σAB between the heterogeneous LJ argon was set as σAB = σAA = σBB. The leapfrog Verlet method was used to integrate the equations of motion numerically. Each run comprised integration steps with a time increment Δt of 5 fs. The total time required for each simulation was 1 μs, which was confirmed as a sufficiently long time in a previous study [29]. The reduced temperature T* = kBT/εAA was set as between 0.7 and 1.0. The lowest temperature T * = 0.7 was set considering that the triple point of pure LJ liquid is T * = 0.685-0.687 [53]. The temperature of the system was controlled by velocity scaling in the usual manner once every 100 steps. Thermal control was applied up to 0.9 μs, after which there was almost no fluctuation in temperature, which was almost equal to the temperature set in the simulation. We did not observe any anomalous increase in the velocity of some particles associated with the adoption of the temperature change algorithm instead of Nosé-Hoover for the NVT ensemble. In addition, we did not adopt the data from the period of the temperature change algorithm (0-0.9 μs) but only analyzed the data from the subsequent period of the NVE ensemble (0.9-1.0 μs). The number of particles, the cutting distance, the interface area, and the total simulation time were very large and long compared to previous studies.

Analysis
The reduced local density * and the reduced local interfacial tension L * on sliced spaces of dimensions Lx × 0.1 σAA × Lz were determined by calculating the running average of 0.95-1.00 μs once every 10 steps [25]. In all steps in the calculation of * and L * , the center of gravity in the y direction of the entire system was calculated and the instantaneous distribution of * and L * was shifted so that the center of gravity was y = 0, meaning that the oscillations of the liquid film in the y direction did not affect * and L * . The value of * was calculated for the sliced spaces as follows: where <Ny> is the average number of argon particles "A" or "B" in the sliced space of

Analysis
The reduced local density ρ * y and the reduced local interfacial tension γ * L on sliced spaces of dimensions L x × 0.1 σ AA × L z were determined by calculating the running average of 0.95-1.00 µs once every 10 steps [25]. In all steps in the calculation of ρ * y and γ * L , the center of gravity in the y direction of the entire system was calculated and the instantaneous distribution of ρ * y and γ * L was shifted so that the center of gravity was y = 0, meaning that the oscillations of the liquid film in the y direction did not affect ρ * y and γ * L . The value of ρ * y was calculated for the sliced spaces as follows: where <N y > is the average number of argon particles "A" or "B" in the sliced space of volume V y . The equimolar dividing surface was determined as the position in the y direction where the ρ * y of argon "A" was equal to the ρ * y of argon "B". Furthermore, the bulk densities of each component ρ * A and ρ * B were determined from the ρ * y of argon "A" and "B" in the ranges of 0 < |y| < 2.5 σ AA , where y = 0 is the center of the three-phase LJ liquid film and was sufficiently distant to be unaffected by the anisotropic structure of the interface. The density ratio of argon "A" and "B" represented the saturation solubility.
The reduced local interfacial tension γ * L was calculated using the following statistical mechanical equation [24]: where r ij is the distance between particles "i" and "j", and y ij is the distance in the ydirection. u ij is the LJ interaction potential between particles "i" and "j". This equation shows the anisotropy of the stresses, and if the same stresses are applied in all three directions, the local interfacial tension is zero. This equation is an approximation and is known to generate artifacts in the local normal pressure in the y-direction and the profiles of the local surface tension, such as local negative tension at the gas-liquid interface, at a position closer to the gas than the equimolar dividing surface. However, by integrating these local artifacts in the thickness direction, the correct overall interfacial tension can be obtained [54]. In addition, this effect is very weak for the liquid-liquid interface, and the effect on the position of the tension surface seems to be limited. The two peaks of local tension were integrated in the y-direction to account for the thickness of the sliced space (0.1 σ AA ), and then the interfacial tension γ* of the liquid-liquid interface was obtained. In the symmetric binary LJ system, two surfaces of tension exist for an equimolar dividing surface, and the sum of the two surfaces of tension should be recognized macroscopically as the interfacial tension for each equimolar plane, meaning that the range of integration includes the two tension surfaces [29]. The Tolman length δ* = δ/σ AA is the distance between the equimolar dividing surface and the surface of tension, which is the peak position of the local interfacial tension in the y-direction. In this study, two equimolar dividing surfaces and four surfaces of tension occurred, so the Tolman length was the average of four distances between the two equimolar dividing surfaces and the four surfaces of tension.

Results
MD simulations were conducted for several combinations of α and T* to obtain the LJ liquid-liquid phase equilibrium. Typical distributions of local density and local surface tension are shown in Figure 2. The solid line shows the local density distribution of argon "A", the dotted line is argon "B", and the dashed line is the total density of the two argons "A" and "B". In the bulk region of the argon "A" phase, the background is filled in with gray. The position where the local densities of argon "A" and "B" coincide is the equimolar dividing surface (blue lines). As shown in Figure 2, the local surface tension was almost zero at the equimolar dividing surfaces and, in the bulk phase of LJ argon "A" and "B", was sufficiently far from the equimolar dividing surfaces. Two tension surfaces were present on both sides of each equimolar dividing surface (red lines). Incidentally, the large local interfacial tensions that occurred at the gas-liquid interface, as shown in Figure 2, were not considered in this study because they corresponded to the surface tensions of the gas-liquid interface. The four peaks integrated in the y-direction for the calculation of the interfacial tension γ* of two interfaces are filled in with shaded lines.
The obtained ρ * A and ρ * B are listed in Table 1 and shown in Figures 3 and 4. At α = 0.50, each of the three argon LJ liquid films was isolated and did not contact with each other at all ("complete separation" is indicated in Table 1). When α was greater than 0.85, LJ argon "A and "B" were completely mixed without creating interfaces, even at T* = 0.70 ("complete mixture" is indicated in Table 1). Between α = 0.55 and 0.80, the partial mixing of LJ argon "A" and "B" was observed. Between α = 0.55 and 0.80, the greater the value of α, the lower the temperature required for partial mixing. For all α, when T* was higher than 1.05, no interface was observed at all due to the perfect mixing of heterogeneous LJ argons. For a smaller value of α, the weak interactions between the LJ heterogeneous argons prevented them from mixing at lower temperatures. In the high temperature region, the increase in entropy caused some mixing. The greater the value of α, the greater the depression in temperature for complete mixing. Therefore, when dealing with symmetric binary LJ liquids and conducting capillary phase separation in pores, it is necessary to use settings within this limited range of α and T*.     Next, the obtained γ* and δ* are listed in Table 1 and shown in Figures 5 and 6. A higher T* or bigger ε AB resulted in a lower interfacial tension γ* and longer Tolman length δ*. In particular, δ* was greater than 3 in the longest case. In other words, the liquid-liquid equilibrium changed to complete mixing when the δ* was slightly larger than 3. In contrast, the δ* of the single component LJ liquid-gas interface was about 1 at the triple point and slightly more than 2 at T* = 1 [29]. The δ* of the binary LJ liquid-liquid interface was obviously larger than that of the single LJ liquid-gas interface.   Figure 5 shows that γ* changed strongly under the influence of T*, but at the same time, ρ * B also changed strongly under the influence of T*, as shown in Figure 4. In other words, since γ* may also be strongly affected by ρ * B , the effects of ρ * B and T* on γ* are shown in Figure 7. Comparing γ* between similar ρ * B , γ* tends to be slightly larger with a higher T*, but the effect of T* on γ* is small. In contrast, the effect of the difference in ρ * B on γ* is very large, and the larger ρ * B is, the smaller γ* becomes.

Discussion
The γL* was nearly plane-symmetric against the position of the peak of the surface of tension, as shown in Figure 2. In other words, even at a distance of δ to 2 δ from the equimolar dividing surface, local interfacial tension occurred; that is, the thickness of the transition region of the LJ liquid-liquid interface was 4 δ, which implies a thickness of roughly 5 σ to 13 σ, based on the values of δ* shown in Table 1.
It is well known that the Tolman length at the vapor-liquid interface of Lennard-Jones fluid is generally between 1 σ and 2 σ [29], depending on the temperature, while the Tolman length at the liquid-liquid interface is much longer than that at the vapor-liquid interface. It is also a peculiar characteristic of the symmetric binary Lennard-Jones system that tension surfaces occur on both sides of the equimolar dividing surface. Since local tension reflects the local stress anisotropy, the tension surfaces should shift to either side of the phase between materials with significantly different densities. In other words, this phenomenon is unique to the symmetric binary system as the bulk density of the two phases is the same, and great care must be taken when using this as a model system for molecular simulations.
This implies that the usual concept of the interfacial tension of the binary LJ liquidliquid at flat interfaces cannot be applied when nano-scale LJ droplets are surrounded by other LJ liquids. For example, as shown in Figure 8, when the diameter of the Ar "A" droplet is 6 σAA and δ is slightly smaller than 3 σAA, the internal surface of tension of the Ar "A" phase shrinks to about 1 σAA or less. Furthermore, inside the surface of tension with a diameter of σ, there is still a weak local interfacial tension, and the local interfacial tensions overlap in the face-to-face direction inside the surface of tension.

Discussion
The γ L * was nearly plane-symmetric against the position of the peak of the surface of tension, as shown in Figure 2. In other words, even at a distance of δ to 2 δ from the equimolar dividing surface, local interfacial tension occurred; that is, the thickness of the transition region of the LJ liquid-liquid interface was 4 δ, which implies a thickness of roughly 5 σ to 13 σ, based on the values of δ* shown in Table 1.
It is well known that the Tolman length at the vapor-liquid interface of Lennard-Jones fluid is generally between 1 σ and 2 σ [29], depending on the temperature, while the Tolman length at the liquid-liquid interface is much longer than that at the vapor-liquid interface. It is also a peculiar characteristic of the symmetric binary Lennard-Jones system that tension surfaces occur on both sides of the equimolar dividing surface. Since local tension reflects the local stress anisotropy, the tension surfaces should shift to either side of the phase between materials with significantly different densities. In other words, this phenomenon is unique to the symmetric binary system as the bulk density of the two phases is the same, and great care must be taken when using this as a model system for molecular simulations.
This implies that the usual concept of the interfacial tension of the binary LJ liquidliquid at flat interfaces cannot be applied when nano-scale LJ droplets are surrounded by other LJ liquids. For example, as shown in Figure 8, when the diameter of the Ar "A" droplet is 6 σ AA and δ is slightly smaller than 3 σ AA , the internal surface of tension of the Ar "A" phase shrinks to about 1 σ AA or less. Furthermore, inside the surface of tension with a diameter of σ, there is still a weak local interfacial tension, and the local interfacial tensions overlap in the face-to-face direction inside the surface of tension. The interfacial tension results also indicate that it is difficult to bring the usual concept of interfacial tension to bear on capillary phase separation inside nanopores. A schematic diagram of the application of the usual concept of interfacial tension to capillary phase separation, also in nanopores, is shown in Figure 9. Nanopores are pores with a size at the lower mesopore limit, at 2 to 4 nm [29]. This means that the reduced pore size is approximately 5 σ to 10 σ. Inside the nanopore, the liquid with a strong interaction with the pore wall is capillary condensed as phase B, and the liquid with a weak interaction with the pore wall exists outside of it as phase A. There is also a phase B covering the surface of the pore wall [50]. As shown in Figure 9, one tension surface will be embedded inside the pore wall, and the other tension surface will be located close to another pore wall. Since the liquid-liquid interface forms a meniscus in capillary phase separation, the position of the tension surfaces inside phase A becomes strange. This means that the concept of the Tolman length at a flat interface may not be applicable to the capillary phase separation of symmetric binary liquids.  The interfacial tension results also indicate that it is difficult to bring the usual concept of interfacial tension to bear on capillary phase separation inside nanopores. A schematic diagram of the application of the usual concept of interfacial tension to capillary phase separation, also in nanopores, is shown in Figure 9. Nanopores are pores with a size at the lower mesopore limit, at 2 to 4 nm [29]. This means that the reduced pore size is approximately 5 σ to 10 σ. Inside the nanopore, the liquid with a strong interaction with the pore wall is capillary condensed as phase B, and the liquid with a weak interaction with the pore wall exists outside of it as phase A. There is also a phase B covering the surface of the pore wall [50]. As shown in Figure 9, one tension surface will be embedded inside the pore wall, and the other tension surface will be located close to another pore wall. Since the liquid-liquid interface forms a meniscus in capillary phase separation, the position of the tension surfaces inside phase A becomes strange. This means that the concept of the Tolman length at a flat interface may not be applicable to the capillary phase separation of symmetric binary liquids. The interfacial tension results also indicate that it is difficult to bring the usual concept of interfacial tension to bear on capillary phase separation inside nanopores. A schematic diagram of the application of the usual concept of interfacial tension to capillary phase separation, also in nanopores, is shown in Figure 9. Nanopores are pores with a size at the lower mesopore limit, at 2 to 4 nm [29]. This means that the reduced pore size is approximately 5 σ to 10 σ. Inside the nanopore, the liquid with a strong interaction with the pore wall is capillary condensed as phase B, and the liquid with a weak interaction with the pore wall exists outside of it as phase A. There is also a phase B covering the surface of the pore wall [50]. As shown in Figure 9, one tension surface will be embedded inside the pore wall, and the other tension surface will be located close to another pore wall. Since the liquid-liquid interface forms a meniscus in capillary phase separation, the position of the tension surfaces inside phase A becomes strange. This means that the concept of the Tolman length at a flat interface may not be applicable to the capillary phase separation of symmetric binary liquids.  In the liquid phase consisting of a single component, the surface tension-temperature relationship is well known, and the Guggenheim-Katayama equation accurately shows the experimental results for neon, argon, oxygen, and nitrogen over a wide range of temperatures [55]: where V is the molar volume and the subscript C indicates a critical point. Since molecular dynamics simulations of liquid-liquid interfacial tensions, such as the present study, are unprecedented, we investigated whether the Guggenheim-Katayama could can be directly applied to symmetric binary Lennard-Jones liquid-liquid interfacial tensions as a practical approximation.
Here, we defined T c for a symmetric binary Lennard-Jones liquid-liquid system as the temperature at which the two liquid phases mix completely with each other. Equation (4) can be used to determine the critical temperature in the single-component system [56]: where ρ is density, and the subscripts L and V indicate liquid and vapor phases, respectively. C is the adjustable parameter used to determine the critical point. Thus, for the binary Lennard-Jones system, we assume the following similar equation: According to equation 5, using C determined from ρ 1 and ρ 2 by the least-squares method, ρ C values were obtained as 0.374-0.387, 0.344-0.359, 0.337-0.351, and 0.322-0.337 for values of α of 0.80, 0.75, 0.70, and 0.65, respectively.
In Figure 10, the γ* and (1 − T/T C ) 11/9 relationship is shown in order to examine the applicability of the Guggenheim-Katayama equation for the interfacial tension of a binary LJ liquid system. In cases in which T is close to T c and α = 0.65 or 0.70, the decay pattern becomes nonlinear and does not follow the Guggenheim-Katayama equation. Surprisingly, the proportionality of the Guggenheim-Katayama equation was almost valid for the MD simulation results for α = 0.75 or 0.80, even though it is a rather rough approximation. As is well known, the phase diagram and physical properties of a one-component LJ fluid are non-dimensionalized by the LJ parameters ε and σ, such as T* and γ* [54]. Even for LJ materials with different ε and σ, the phase diagram and physical properties will be the same due to non-dimensionalization. On the other hand, T C indicates the temperature at which perfect mixing occurs and depends on α, the interaction strength between different LJ Ar. α is not involved in the non-dimensionalization of T* and γ*, and so such non-dimensionalization is not suitable for T/T C ; thus, non-dimensionalization by ε and σ is not used here.

Conclusions
The interfacial tension and Tolman length of the liquid-liquid interface in symmetric binary LJ fluid systems, which have not been studied before, are revealed in this work with a sufficient cutting distance and sufficient interfacial area. Two surfaces of tension were observed on both sides of the equimolar dividing surface of a partially miscible LJ mixture. The higher the mutual solubility, the longer the Tolman length becomes, exceeding 3 σ. Therefore, it is necessary to recognize that applying the concept of the liquid-liquid interface by continuum approximation or the concept of the Tolman length to the curved liquid-liquid interface at the nanoscale may require unreasonable assumptions to be made.