Logarithmic finite-size effects on interfacial free energies: Phenomenological theory and Monte Carlo studies

The computation of interfacial free energies between coexisting phases (e.g., saturated vapor and liquid) by computer simulation methods is still a challenging problem due to the difficulty of an atomistic identification of an interface and interfacial fluctuations on all length scales. The approach to estimate the interfacial tension from the free-energy excess of a system with interfaces relative to corresponding single-phase systems does not suffer from the first problem but still suffers from the latter. Considering $d$-dimensional systems with interfacial area $L^{d-1}$ and linear dimension $L_z$ in the direction perpendicular to the interface, it is argued that the interfacial fluctuations cause logarithmic finite-size effects of order $ln\left(L\right)/L^{d-1}$ and order $ln\left(L_z\right)/L^{d-1}$, in addition to regular corrections (with leading-order $const/L^{d-1}$). A phenomenological theory predicts that the prefactors of the logarithmic terms are universal (but depend on the applied boundary conditions and the considered statistical ensemble). The physical origin of these corrections are the translational entropy of the interface as a whole, “domain breathing” (coupling of interfacial fluctuations to the bulk order parameter fluctuations of the coexisting domains), and capillary waves. Using a new variant of the ensemble switch method, interfacial tensions are found from Monte Carlo simulations of $d=2$ and $d=3$ Ising models and a Lennard-Jones fluid. The simulation results are fully consistent with the theoretical predictions.

Figure 1

Useful boundary conditions to study interfaces in Ising models, using a simulation box of linear dimension(s) $L$ parallel to the interface(s) and $L_z$ in the perpendicular direction. In the parallel direction(s), periodic boundary conditions (PBC) are used in all cases. Interfaces are schematically indicated by thick wavy lines, and the double arrows indicate the sign of the magnetization in the coexisting domains $(m_{+},m_{-})$. Note that $\langle m_{+}\rangle =m_{\mathrm{coex}},\langle m_{-}\rangle =-m_{\mathrm{coex}},m_{\mathrm{coex}}$ being the spontaneous magnetization of the Ising model in the thermodynamic limit, while in finite simulation boxes $m_{+},m_{-}$ may fluctuate around these average values. Case (a) assumes free surfaces at the first $(n=1)$ and last $(n=L_z)$ layer, with a fixed spin boundary condition for all spins in the adjacent layers, $S_i=+1$ for all spins in the layer $n=0,S_i=-1$ in the layer $n=L_z+1$. Case (b) indicates the antiperiodic boundary condition (APBC) and case (c) uses also a periodic boundary condition in the $z$ direction.

Figure 2

Density distribution $P_{L,L_z}\left(\rho \right)$ plotted vs $\rho$ for the two-dimensional Ising (lattice gas) model at temperature $k_{\mathrm{B}}T/J=2.0$ for the case (a) $L=20$ and varying $L_z,L_z=30,40,50,60,70,80,90$ and 100 (from bottom to top at $\rho =0.5)$, for the case (b) $L_z=120$, and varying $L,L=10,20,30,40,50,60,70,80$, and 90 (from top to bottom at $\rho =0.5)$, and for the case (c) $L=6$ and increasing $L_z,L_z=24,48,96,192,384$ (from bottom to top at $\rho =0.5)$, illustrating the change from the double-peak distribution for rather small $L_z$ $\left(L_z=24\right)$ to a three-peak distribution (e.g., $L_z=96,192)$ to a distribution with a single central peak $\left(L_z=384\right)$.

Figure 3

(a) Coarse-grained view of a single interface on a length scale $L_z$, to illustrate the counting of configurations to estimate its translational entropy. The scale $L_z$ is partitioned into cells of width $w_L$. Coarse-graining in the $x$ direction over a length $l_x$ of order $\xi$ (not shown), to eliminate overhangs, bubbles, etc., present in the original microscopic configurations [see parts (b), (c)], one obtains a smooth interface with intrinsic width $w_0$ comparable to $l_x$. This coarse-grained interface shows undulations on all length scales from $l_x$ to $L$ (capillary waves) and thereby exhibits a width $w_L>w_0$. [(b) and (c)] show two snapshots resulting from a simulation with $L=60,L_z=120,k_{\mathrm{B}}T/J=2.0$, APBC in the $z$ direction, and using the grand-canonical ensemble. The interface position “explores” the entire length $L_z$ of the system. Up-spins are shown as black dots, and down-spins are not shown.

Figure 4

Four snapshots of two-dimensional Ising systems at fixed magnetization $m=0$ to visualize domain breathing. Spontaneous fluctuations of the magnetization densities inside the two on-average equal-sized domains $(L_x=L_z/2=60)$ have an effect on the location of the interface, when the total magnetization is conserved. In the shown example the typical length $\sqrt{{\langle \Delta \rangle }^2}$ over which the interface position fluctuates is comparable to the lattice spacing, and the shift of the interface is barely visible in these snapshots. However, one can see rather large fluctuations in the bulk phases, which are the reason for domain breathing. All snapshots are for $k_{\mathrm{B}}T/J=2.0$. APBC in the $z$ direction (horizontal) and PBC in the $x$ direction (vertical) are used.

Figure 5

Schematic explanation of the “ensemble switch method” to find the interfacial free energy. A system is constructed as a linear combination of two Hamiltonians, $\mathcal{H}\left(\kappa \right)=\kappa {\mathcal{H}}_1+(1-\kappa ){\mathcal{H}}_0$, where ${\mathcal{H}}_1$ is the desired system of interest (containing two interfaces in the case shown here) and ${\mathcal{H}}_0$ consists of two separate systems of half the linear dimension $L_z/2$ each, and periodic boundary conditions. The mixing parameter $\kappa$ is in the interval $0\le \kappa \le 1$, and the free-energy difference between the states with $\mathcal{H}(\kappa =0)$ and $\mathcal{H}(\kappa =1)$ yields twice the interfacial free energy. If the state $\kappa =1$ is a state with an APBC, one can obtain the free energy associated with a single interface analogously.

Figure 6

Free-energy difference $\Delta F\left(\kappa \right)$ versus $\kappa$, as obtained from the ensemble switch method for $L\times L_z$ Ising systems. (a) Data for $d=2$ and $L=40$; (b) data for $d=3$ and $L=14$. Different choices of boundary conditions (PBC, APBC) and ensembles [grand-canonical (gc) and canonical (c)] are compared for several $L_z$. The order of the curves at $\kappa =1$ from top to bottom is as indicated, the bottom three curves being for APBC(gc), the middle three curves for PBC(c), and the top curves for APBC(c), with increasing $L_z$ from top to bottom. Note that only a small section of the whole variation of $\kappa$ and of $\Delta F\left(\kappa \right)$ is shown to display the finite-size effects on $\Delta F\left(\kappa \right)$ clearly.

Figure 7

Estimates for the interfacial tension ${\gamma }_{L,L}$ of the $d=2$ Ising model at $k_{\mathrm{B}}T/J=2.0$ plotted vs $1/L$. Two sets of data from the probability distribution method are included (one from Ref. [46], one from the present work), which agree within statistical errors with each other (the statistical error is smaller than the size of the symbols throughout). The third set of data is due to the ensemble switch method and has in this case slightly larger finite-size effects than the method based on Eq. (27).

Figure 8

(a) Interfacial tension ${\gamma }_{L,L_z}$ for the $d=2$ Ising model plotted at fixed $L$ versus $L^{-1}ln\left(L_z\right)$ for the cases APBC(c), PBC(c), and APBC(gc). The upper set of data refers to $k_{\mathrm{B}}T/J=1.2,L=10$, the lower set of data refers to $k_{\mathrm{B}}T/J=1.6,L=10,20$, and 30, as indicated. The straight lines shown have the theoretical slopes $x_{\perp }=1/2$, 3/4, and 1, respectively. (b) Same as (a), but showing only the PBC(c) case for the three temperatures $k_{\mathrm{B}}T/J=1.2,1.6$, and $2.0$, respectively. Always three choices of $L$ are shown, namely $L=10,20$, and 40 (from top to bottom). The straight lines illustrate the theoretical slope $x_{\perp }=3/4$ throughout.

Figure 9

(a) Interfacial tension ${\gamma }_{L,L_z}$ plotted vs $L$ for the $d=2$ Ising model at $k_{\mathrm{B}}T/J=1.2$ for two choices of $L_z,L_z=60$ and $L_z=120$; the horizontal straight line shows the known value of ${\gamma }_{\infty }$ [Eq. (6)] [111], while the curves are fits of Eq. (26) to the data (symbols) for the cases APBC(c), top set of curves; PBC(c), middle set of curves; APBC(gc), bottom set of curves. In each set $L_z$ increases from top to bottom. The theoretical values $x_{\perp },x_{\parallel }$ from Table 1 were used in the fit. (b) Reduced interfacial tension $\tilde{\gamma }$ [Eq. (34)] plotted vs $1/L$, for $k_{\mathrm{B}}T/J=1.2$, and three choices of boundary conditions and/or ensembles, as indicated [PBC(c), APBC(c), and APBC(gc)]. In each case two choices of $L_z$ are included, $L_z=60$ and $L_z=120$. Symbols represent the simulation results, and straight lines show the fits $\tilde{\gamma }={\gamma }_{\infty }+C/L$; the fitted values ${\gamma }_{\infty },C$ are quoted in the figure. The horizontal straight line shows the known exact result, ${\gamma }_{\infty }=1.284$ [from Eq. (6)] [111]. (c) Same as (b) but for $k_{\mathrm{B}}T/J=1.6$. Here the exact result is ${\gamma }_{\infty }=0.660$ [111]. (d) Same as (b) but for $k_{\mathrm{B}}T/J=2.0$. Here the exact result is ${\gamma }_{\infty }=0.228$ [111].

Figure 10

(a) Interfacial tension ${\gamma }_{L,L_z}$ for the $d=3$ Ising model at $k_{\mathrm{B}}T/J=3.0$ is plotted vs the variable $(1/L^2)ln\left(L_z\right)$ using PBC and the canonical ensemble. Several choices of $L$ are included, as indicated. The straight lines show the theoretical exponent $x_{\perp }=3/4$ (resulting from the entropy of interface translation and domain breathing). (b) Interfacial tension ${\gamma }_{L,L_z}$ for the $d=3$ Ising model at $k_{\mathrm{B}}T/J=3.0$ plotted vs $L$ for three choices of $L_z,L_z=20,40$ and 80. The horizontal straight line shows the previous result ${\gamma }_{\infty }=0.434$ due to Hasenbusch and Pinn [49], while the curves are fits of Eq. (26) to the data (symbols) for the cases APBC(c), top set of curves; PBC(c), middle set; APBC(gc), bottom set. In each set, $L_z$ increases from top to bottom. The theoretical values of $x_{\perp },x_{\parallel }$ from Table 1 were used, so each curve contains a single adjusted constant (the prefactor of the $1/L^2$ term) only.

Figure 11

Reduced interfacial tension $\tilde{\gamma }$ [Eq. (35)] plotted vs $1/L^2$ for the $d=3$ Ising model, using three chocies of $L_z$ and three choices of boundary conditions and/or ensembles, as indicated. Case (a) includes the parameter $C_1$, in the fit, while case (b) requires $C_1=0$. The resulting estimates for the parameters ${\gamma }_{\infty }$ and $C_2$ are quoted in the figure.

Figure 12

Interfacial tension of the LJ fluid at $k_{\mathrm{B}}T/\varepsilon =0.78$ plotted vs the scaling variable $L^{-2}ln\left(L_z\right)$, for three choices of the cross-sectional area $A=L^2$, as indicated. The slope of the straight lines is again the theoretical value, $x_{\perp }=3/4$. The circles show preliminary data for rather large $L_z$ with insufficient statistics (see text).

Figure 13

Interfacial tension of the LJ fluid at $k_{\mathrm{B}}T/\varepsilon =0.78$ plotted vs $L^{-2}$, using either local or nonlocal moves, for fixed $L_z$, as indicated. The logarithmic corrections have not been subtracted. The lines are fits of the form as in Eq. (35).

Figure 14

Reduced interfacial tension $\tilde{\gamma }$ [Eq. (35)] of the Lennard-Jones fluid at $k_{\mathrm{B}}T/\varepsilon =0.78$ plotted vs $L^{-2}$. In case (a) the parameter $C_1$ was forced to be zero, while in case (b) both parameters $C_1,C_2$ were fitted. The results of the fits are quoted in the figure. Note that both data dare included where only local moves of the particle were permitted, as well as data where randomly chosen particles were removed from their position and reinserted at a randomly chosen position anywhere in the box.

Figure 15

Interfacial tension of the symmetric LJ fluid at $k_{\mathrm{B}}T/\varepsilon =1.0$ plotted vs $L^{-2}$. Cubic boxes of volume $L^3$ have been used. The lines are fits of the form as in Eq. (35). The original data are from Ref. [119], in which a fit of the form $A+B/L$ (corresponding to fit with $C_2=0$ and $C_1=-0.13)$ suggests ${\gamma }_{\infty }\approx 0.722$. After subtracting the logarithmic contributions, the result is ${\gamma }_{\infty }\approx 0.717$. Note that the statistical errors of these data are hard to estimate precisely but could be of the order of a few percentages. Hence the scatter of the data points is well within the statistical errors.