Chasing the Critical Wetting Transition. An Effective Interface Potential Method

Abstract

Wettablity is one of the important characteristics defining a given surface. Here we show that the effective interface potential method of determining the wetting temperature, originally proposed by MacDowell and Müller for the surfaces exhibiting the first order wetting transition, can also be used to estimate the wetting temperature of the second order (continuous) wetting transition. Some selected other methods of determination of the wetting temperature are also discussed.

1. Introduction

Understanding and controlling wetting properties of materials are some of the most important factors in many industrial applications including oil recovery [1], mineral flotation [2] and design of superamphiphobic surfaces [3].

Wetting transition is a surface-induced transition in which the contact angle of a liquid deposited on a surface drops to zero from a non-zero value upon increasing temperature. This transition has been the subject of many theoretical and experimental studies [4,5,6,7,8,9,10,11,12]. The wetting transition can be either first order or continuous (second-order). In the case of the first-order wetting transition there is additional prewetting (thin-thick film) transition, which occurs off-coexistence and the prewetting line joins the binodal exactly at the wetting point. The critical (second-order) wetting transition is not accompanied by prewetting.

The nature of the wetting transition and its universality class depends on (among others) the dimensionality of the system and the range of interparticle interactions. While in most cases the nature of the wetting transition has been well established and documented, it turns out that the case of 3-dimensional critical wetting transition for short-ranged forces (i.e., decaying exponentially, or faster) posed significant problems.

One of the important early discoveries was that the critical wetting transition for short-range forces is nonuniversal. Using renormalization-group calculations based on an effective interfacial Hamiltonian Brezin et al. [8] predicted that the critical exponent for this transition depends on a dimensionless parameter

(1)$\omega =\frac{k_{B}T}{4\pi \Sigma {\xi }^2}$

where T is the temperature, $\xi$ is the correlation length of the bulk wetting phase and $\Sigma$ is the interfacial stiffness. Capillary-wave-like fluctuations give rise to a diverging transverse correlation length ${\xi }_{\parallel }$ and the relevant exponent is believed to behave

(2)${\nu }_{\parallel }=\left\{\begin{array}{cc}{(1-\omega )}^{-1}\hfill & \mathrm{for}0<\omega <1/2\hfill \\ {(\sqrt{2}-\sqrt{\omega })}^{-2}\hfill & \mathrm{for}1/2<\omega <2\hfill \\ \infty \hfill & \mathrm{for}\omega >2.\hfill \end{array}\right.$

Unfortunately, subsequent Monte Carlo calculations carried out for the Ising model [13,14,15] revealed that while the general features regarding the wetting transition for the Ising model do agree with theoretical predictions [9], the critical wetting transition is only very weakly nonuniversal. This disagreement has been the subject of ever-lasting efforts in order to bridge the gap between the theory and simulations. Using nonlocal effective interfacial Hamiltonian Parry et al. [16] argued that the spectrum of the interfacial fluctuations has a lower cutoff due to appearance of a new length scale ${\xi }_{NL}=\sqrt{l\xi }\propto \sqrt{ln{\xi }_{\parallel }}$. This gives rise to an effective wetting parameter, ${\omega }_{eff}$, of the form

(3)${\omega }_{eff}=\omega -\sqrt{2{\omega }^3}\frac{ln\left(\kappa l\right)}{\kappa l}$

where $\kappa ={\xi }^{-1}$, and l denotes film thickness. This leads to lowering of the value of the effective wetting parameter and yields lower effective critical exponent.

On the simulational front Albano and Binder [17] suggested that anisotropic finite size scaling (AFSS) theory should be suitable for studying wetting transitions. Using this approach Bryk and Binder [18] were able to recalculate the location of the wetting transition and confirm non-mean field character of critical wetting in 3D.

The main obstacle that hampered progress in our understanding of the nature of the critical wetting transition could be traced back to difficulties in accurate locating the critical wetting point from simulation. In the present work we discuss three methods of locating critical wetting from simulations. We show that the effective interface method can be used for locating the critical wetting transition. Two other alternative methods are also discussed.

2. Materials and Methods

Our model consists of cubic lattice of dimensions $L\times L\times D$ with two free boundary layers $L\times L$ located at $z=1$ and $z=D$, and periodic boundary conditions in the remaining directions. The pseudospin variable at a lattice site i takes values $s_i=\pm 1$. The Hamiltonian for the system is

(4)$\mathcal{H}=-J\sum _{bulk}{s}_i{s}_j-J_s\sum _{surf}{s}_i{s}_k-H\sum _{bulk}{s}_i-H_1\sum _{k\in z=1}{s}_k-H_D\sum _{k\in z=D}{s}_k.$

In the two free surface layers the exchange constant is $J_s$, otherwise the exchange constant is J throughout. The bulk field is H, and the surface fields acting on the first and last layer are $H_1$ and $H_D$, respectively. We consider three different types of systems, namely the symmetric systems with $H_1=H_D$, non-symmetric systems with $H_1\ne H_D$, $H_D=0$, and the anti-symmetric systems with $H_1=-H_D$. Throughout this study we restrict ourselves to the case $J_s=J$.

The systems were simulated using fast multispin coding algorithm [19]. In order to achieve better statistics we used the preferential sampling technique, so that, on average, 9 out of 10 samplings occurred in the region of interest (in the vicinity of the walls).

When simulating systems close to the critical wetting point the so-called critical slowing down hampers the statistics of the accumulated data. In order to overcome this drawback we applied hyper-parallel tempering technique [20] and simulated many systems at the same time. The swaps of spins between the systems m and n were accepted with the probability

(5)$P_{nm}=\min [1,exp(-\Delta \beta \Delta E)-\Delta \left(\beta H_1\right)\Delta m_1-\Delta \left(\beta H_D\right)\Delta m_D-\Delta \left(\beta H\right)\Delta m]$

where

(6)$\Delta \beta =\frac{1}{k_B{T}_n}-\frac{1}{k_B{T}_m},\Delta E=E_n-E_m$

(7)$\Delta \left(\beta H_1\right)=\frac{1}{k_B{T}_n}{H}_1^{\left(n\right)}-\frac{1}{k_B{T}_m}{H}_1^{\left(m\right)},\Delta m_1=m_1^{\left(n\right)}-m_1^{\left(m\right)}$

(8)$\Delta \left(\beta H_D\right)=\frac{1}{k_B{T}_n}{H}_D^{\left(n\right)}-\frac{1}{k_B{T}_m}{H}_D^{\left(m\right)},\Delta m_D=m_D^{\left(n\right)}-m_D^{\left(m\right)}$

(9)$\Delta \left(\beta H\right)=\frac{1}{k_B{T}_n}{H}^{\left(n\right)}-\frac{1}{k_B{T}_m}{H}^{\left(m\right)},\Delta m=m^{\left(n\right)}-m^{\left(m\right)}.$

In the above $E_n$ and $T_n$ are the energy and the temperature of the system n. $k_B$ is the Boltzmann constant. Among other quantities of interest accumulated during a simulation were the magnetization in the surface layer

(10)$m_1={\left(2L^2\right)}^{-1}\sum _{k\in surf1}<s_k>,$

total magnetization

(11)$m={\left(L^{2}D\right)}^{-1}\sum _i<s_i>,$

and “mixed” surface layer susceptibility

(12)${\chi }_1=\frac{\partial m_1}{\partial H}=L^{2}D(<m_{1}m>-<m_1><m>)/k_{B}T.$

3. The Effective Interface Potential Method

Tracking down the critical wetting transition point from computer simulation proved to be a challenge. MacDowell and Müller (MM), proposed a method [21,22] (hereinafter referred to as the effective interface potential (EIP) method) that relies on determining the distribution probability of magnetization from which the effective interface potential can be determined. We implemented this method with slight modifications. The key quantity is the probability $P\left(m_p\right)$ of finding a system with magnetization $m_p$ defined as

(13)$m_p={\left(D_r{L}^2\right)}^{-1}\sum _{k=1}^{D_r}<s_k>,$

where $D_r$ denotes the range of interest, $D/2\le D_r\le D$. We considered the non-symmetric system with $H_D=0$. Following [21,22] we split the calculations of the order parameter distribution in windows and applied the successive sampling technique but with windows on $m_p$ rather than on m. Since we used hyper-parallel tempering two or three windows were deemed sufficient. From the order parameter distribution we obtained the effective interaction potential (up to a constant) via

(14)$V_{eff}\left(m_p\right)/k_{B}T=-ln\left(P\left(m_p\right)\right).$

In our work, we used $P\left(m_p\right)$, while MacDowell and Müller [21,22] calculated the adsorption distribution $P(\Gamma )$. However the two distributions are closely related as one can recover the adsorption from $m_p$ by subtracting bulk magnetization.

It is important to note that the effective interface potential employed by us $V_{eff}\left(m_p\right)$ differs from that used traditionally in the literature [4,5,8]. Those authors considered $V_{eff}\left(l\right)$ where l is a distance of the interface from the wall, a local quantity. $V_{eff}\left(m_p\right)$ makes use of $m_p$—a global measure which translates into the mean distance from the wall only within mean field description. In other words our effective interface potential approach includes all the fluctuation effects, like interface overhangs and droplet excitations, which are not included in the local description. Since our interface potential carries a bulk term it cannot be directly used as an interface potential in conventional interface Hamiltonians. In such a case one can use a method proposed in [23].

MacDowell and Müller (and others [24,25,26,27,28,29,30]) have applied the outlined above approach to the long-range substrate potentials exhibiting the first-order wetting transition, where one expects the mean-field effective interface potential to be sufficient. It is natural to ask whether this approach works at all beyond mean-field. In order to answer this question we first carried out MC calculations for the 2D Ising model with short-range boundary fields. It is well established that for short-range forces in 2D wetting is completely dominated by the fluctuation effects [4]. An advantage of this procedure is that the simulational results can be compared with the exact solution for the wetting ordering field by Abraham [7]

(15)$exp(2J/k_{B}T)[cosh(2J/k_{B}T)-cosh(2H_{1c}/k_{B}T)]=sinh(2J/k_{B}T).$

These calculations provide a clear-cut, stringent test of the proposed approach.

Figure 1, Figure 2 and Figure 3 show the EIPs calculated for three system sizes, $L=210$, 630 and 1260 at $J/k_{B}T=1.0$. For smaller absolute values of the surface fields we observe a deep local minimum indicating that the system is in the partial wetting state. For these state points the successive sampling technique is particularly useful as it would be very difficult to obtain smooth effective interface potentials only from one simulation. As $H_1$ becomes more negative this minimum gradually becomes shallower and finally disappears completely indicating the critical wetting point for a system with the size L.

Figure 4 shows extrapolation of the critical surface field to the thermodynamic limit $L\to \infty$. We note, that the simulational value of $H_{1c}(\infty )/J=-0.935$ agrees very well with the exact result, −0.927. Another important observation is that the critical surface field $H_{1c}\left(L\right)$ depends quite visibly on L, indicating that finite size corrections are important. This is in line with theoretical predictions of 2D critical wetting [4], i.e., the interfacial fluctuations in 2D renormalize the temperature (or alternatively the critical surface field) of the critical wetting transition.

It is instructive to see how accurate is the EIP method at higher temperatures (recall, that the inverse bulk critical temperature $J/k_B{T}_c=ln(1+\sqrt{2})/2\approx 0.440687$). Figure 5 shows extrapolation of the critical surface field to the thermodynamic limit at $J/k_{B}T=0.625$ (the relevant effective potentials are not shown, for the sake of brevity). We observe that the simulational result $H_{1c}/J=-0.735$ now differs from the exact result ($-0.71717$) by almost 2.5%. This deterioration in accuracy can be traced back to the fact that in 2D (unlike in 3D) bulk fluctuations persist down to quite low temperatures. Since these fluctuations lead to broadening of $P\left(m_p\right)$ it is inevitable that the method becomes less robust as the temperature gets closer to the bulk critical temperature. Hence, one important restriction of the EIP method is that it is applicable to systems with well separated bulk and interfacial fluctuations.

Let us turn to the 3D Ising systems. Figure 6, Figure 7 and Figure 8 show the effective interface potentials calculated for $J/k_{B}T=0.35$, $D=60$, and for $L=63$, 126, and 252. We notice that even for the smallest system, $L=63$ (cf. Figure 6), the system with the surface field $H_1/J=-0.89$ is still in the partial wetting regime as there is a well visible minimum (as a side comment, this value was quoted as the critical surface field in Refs. [13,14,15]). We estimate that the wetting transition is located at $H_1/J=-0.908\pm 0.002$. Interestingly, for the larger systems (cf. Figure 7 and Figure 8) the transition is located also at $H_1/J=-0.908\pm 0.002$, i.e., there are no discernible finite-size effects.

We find similar behaviour of the EIPs for $J/k_{B}T=0.25$. Independent of the system size we estimate the critical surface field $H_{1c}/J=-0.608\pm 0.004$ (the EIPs for the largest considered system size, $L=256$ are presented in. Figure 9). This is again in line with our knowledge about the critical wetting transition, i.e., in 3D the interfacial fluctuations do not renormalize the temperature of the wetting transition [8].

4. Other Methods of Determination of the Critical Surface Field

4.1. Determination of the Critical Surface Field by Thermodynamic Integration

A new method for determination of the contact angle has been proposed in Refs. [31,32]. This technique can be interpreted as a variant of the thermodynamic integration method (TIN), whereby a series of calculations is carried out for anti-symmetric fields, $H_1<0$, $H_D=-H_1$, starting from $H_1=0$. The contact angle can be determined from the relation

(16)$cos\left(\theta \right)=-{\int }_0^{H_D}\frac{m_D+m_1}{{\gamma }_{lv}}d{H}_D^{\prime },$

where ${\gamma }_{lv}$ is the surface tension [33,34]. The advantage of this approach is that the system that has to be simulated does not undergo a critical wetting transition. Consequently there are no diverging length scales and no critical fluctuations. The downside is that for critical wetting the contact angle goes to 0 tangentially,

(17)$1-cos\left(\theta \right)\sim t^{2-{\alpha }_s},$

where ${\alpha }_s$ is the surface critical exponent for specific heat, $t=\frac{H_{1c}-H_D}{H_{1c}}$. This makes the finite-size analysis a delicate issue.

The calculational procedure consists of a series of calculations for the anti-symmetric systems with varying surface field with $\Delta H_1=0.001$. We carried out calculations for $J/k_{B}T=0.25$ and 0.35. The integrand of the right-hand-side of Equation (16) depends quite strongly on the system size close to the critical surface field, despite the fact that there is no wetting transition in the system. Therefore we carried out calculations for $D=160$ and $L=189$ for surface fields $H_1/J=0$ up to $H_1/J=-0.511$, while for the stronger surface fields $H_1/J<-0.511$ the calculations were carried out for $L=189$, 252 and 504. Similar to previous calculations hyper-parallel tempering was used with up to 120 systems simulated at the same time. The averages were accumulated over ${10}^7$ spin flips per site. Figure 10 and Figure 11 show the sum of the surface magnetizations divided by the surface tension, i.e., the integrand of Equation (16). We see that close to the critical surface field the integrand develops finite-size dependence. The insets of Figure 12 and Figure 13 show the cosine of the contact angle vs. the surface field. We see that the cosine goes to 1 tangentially. However, it is difficult to pinpoint exactly the tangential point. Moreover, the integral depends on the accuracy of the numerical values of the surface tension. Hence we adopt the following strategy. First we calculated numerically the derivative of the cosine of the contact angle with respect to the surface field (cf. main plots of Figure 12 and Figure 13). The derivative should be zero at the critical surface field, $H_{1c}\left(L\right)$. In order arrive at an estimate of $H_{1c}\left(L\right)$ we plotted two parallel lines $cos{\theta }^{\prime }=c_1$ and $cos{\theta }^{\prime }=c_2$ and constructed tangential lines via finite difference. Extrapolation of these straight lines give $H_1\left(L\right)$ satisfying $cos{\theta }^{\prime }=0$. Finally from the plots of $H_1\left(L\right)$ vs. $1/L$ the thermodynamic limit is reached via extrapolation (cf. Figure 14). The final values of the critical surface fields estimated using this method are $H_{1c}/J=-0.905$ for $J/k_{B}T=0.35$, and $H_{1c}/J=-0.604$ for $J/k_{B}T=0.25$.

4.2. BLK Method for Symmetric Surface Fields Revisited

The original idea of Binder, Landau, and coworkers [13,14] (hereinafter referred to as the BLK method) was to consider the symmetric system, $H_1=H_D$. After setting $H=0$ the location of the critical wetting transition can be estimated out by varying $H_1$ at a constant temperature. The maximum of ${\chi }_1$ as a function of $H_1$ indicates the location of the critical surface field $H_{1c}$ at which the critical wetting transition occurs. First simulations were carried out for symmetric field at $J/k_{B}T=0.35$, in order to directly compare with [15]. We assumed that $D=60$ and $L=$126, 252, 315, 378, and 504. The statistical effort was $10\times {10}^6$ spin flips per site (not counting extra gains due to the preferential sampling technique). A total of 16 systems, each at different $H_1$ were simulated at once using hyper-parallel tempering technique.

Figure 15 shows the results of the simulations. Its clear that the statistics is much improved, when compared to the earlier efforts.

The data for $L=126$ can be compared with that from Ref. [15]. We note that the position of the maximum of $H_{1c}(L=126)/J=-0.889\pm 0.0005$ agrees very well with earlier study ($H_{1c}(L=128)/J=-0.89\pm 0.004$). However, there is a considerable shift for larger systems. This indicates that $H_{1c}(L=\infty )$ is lower. This is confirmed in Figure 16, where we show extrapolation $L\to \infty$, and find, that the estimated value $H_{1c}(L=\infty )/J=-0.905$. In making this extrapolation we did not include data for $L=126$ due to considerable deviations from the rest of the data.

Figure 17 shows the results calculated for $J/k_{B}T=0.25$. Since this temperature is closer to the bulk critical temperature the correlation length is greater. For this reason we carried out calculations for somewhat larger systems with $D=120$ and $L=$252, 315, 378, and 504. The overall result are similar to those obtained at $J/k_{B}T=0.35$. The extrapolation to $L=\infty$ reveals (cf. Figure 18), that that the estimated value $H_{1c}(L=\infty )/J=-0.606$. This is again, noticeably lower than previous estimates calculated using much smaller system sizes ($H_{1c}(L=\infty )/J=-0.555$ [13,14]).

As the careful Reader noticed, there is one vexing feature of the results presented in Figure 15 and Figure 17. Namely the maxima of ${\chi }_1$ get smaller and smaller as the linear system size increases. This means that in the thermodynamic limit this peak disappears and this transition does not exist! We recall that the simulated system is not semi-infinite but forms a slit-like pore. In such systems the only phase transition that remains stable is capillary condensation. It’s clear that the large statistical effort and hyper-parallel tempering technique must give the correct result, i.e., that the critical wetting transition studied in this simulational setup is not a stable transition in the thermodynamic limit. However, we argue that in this particular case one can still trust the positions of the maxima of the plots of ${\chi }_1$ even if their magnitudes decrease with increasing system size. This situation can be interpreted in terms of finite size scaling at the first order transition, as formulated by Binder and Landau [35]. One can approximate the probability distribution of the magnetization of an Ising system by two Gaussians for the two phases that would coexist at the capillary condensation transition. At zero bulk field, the phase with the magnetization oppositely oriented to the surface fields is not the stable one. It still gives a signal in a finite system, which ultimately-in the thermodynamic limit - will be exponentially suppressed. However, for the surface susceptibility, the stable phase (magnetization parallel to the surface field) gives only a small background contribution. Hence one can still detect the developing singular behaviour of the surface layer susceptibility even though the magnitude of the signal decreases with increasing system size. If the sampling is insufficient to reach full equilibrium, one may get the wrong amplitude of the signal, but it will be still possible to detect the location of the anomaly [36]. To conclude this subsection, the BLK method yields the critical surface fields $H_{1c}(L=\infty )/J=-0.905$ for $J/k_{B}T=0.35$ and $H_{1c}(L=\infty )/J=-0.606$ for $J/k_{B}T=0.25$.

5. Discussion

The results showed in the previous Sections are summarized in Table 1.

We note that all three methods considered in this work yield results consistent with the AFSS method. At the same time there is a huge difference between the results of the original BLK method and the rest of the results. This is however to be expected, since at the time of carrying out the calculations in [13,15], it was technically impossible to go beyond the system sizes considered in those papers. Each of the presented here methods has some advantages and drawbacks. The TIN method while very promising proves to be tricky due to the lack of proper finite-size analysis. The BLK method requires calculations for progressively large system sizes. The EIP method requires the bulk fluctuations to be small (relative to the capillary wave fluctuations), hence it is not suitable for determination of the critical wetting transition close to the bulk critical point.

In terms of the computational effort the TIN method is least demanding. This is due to the fact that the interfacial fluctuations associated with critical wetting are absent, when using this route. Hence the statistical effort measured as the number of attempted spin flips per site can be two orders of magnitude smaller. The other methods are computationally more involving. The EIP method seems to be slightly less demanding than the BLK method due to the fact that there is no need to consider very large system sizes.

6. Conclusions

We have studied three methods of determination of the critical wetting transition. Our findings can be summarized as follows:

• The effective interface potential method can be used to determine the location of the critical wetting transition. The limitation of this method is that its accuracy decreases if the bulk fluctuations become important.

• The thermodynamic integration method can be used to estimate the location of the critical wetting transition. Extrapolation to the thermodynamic limit is non-trivial.

• The Binder–Landau–Kroll method of determination of the critical wetting transition also leads to reasonable results if sufficiently big system sizes are considered.

Recently Evans et al. showed [37] that the Nakanishi and Fisher [9] topology of the global surface phase diagram is not complete. Their study unveils novel classes of the surface phase diagram which are not present for lattice models. We hope that our work will be useful in establishing simulational tools that will help in better understanding of the origin of the differences between the atomistic and lattice models. In our opinion the results of the wetting behaviour for the simple lattice models have to be revisited as well, since the old estimates were computed for small system sizes. Some of these issues are being currently considered.

Author Contributions

Conceptualization, methodology, software, investigation P.B.; data analysis and interpretation, writing, P.B. and A.P.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Polish NCN, grant number OPUS 13 UMO-2017/25/B/ST5/00975.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data used in this research can be obtained from the Authors upon a reasonable request.

Acknowledgments

P.B. acknowledges fruitful discussions with K. Binder and W. Rżysko.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. The effective interface potentials calculated for the 2D non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for the surface fields [Forumla omitted. See PDF.], −0.92, −0.925, −0.93, −0.935, −0.94, −0.942, −0.945, −0.948, −0.950, and −0.955. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 2. The effective interface potentials calculated for the 2D non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for the surface fields [Forumla omitted. See PDF.], −0.912, −0.915, −0.918, −0.920, −0.922, −0.925, −0.928, −0.930, −0.932, −0.935, −0.938, and −0.942. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 3. The effective interface potentials calculated for the 2D non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for the surface fields [Forumla omitted. See PDF.], −0.912, −0.915, −0.918, −0.920, −0.925, −0.928, −0.930, −0.932, −0.935, and −0.938. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 4. The estimate of the wetting surface field for the 2D system at [Forumla omitted. See PDF.] Circles correspond to the surface critical fields obtained from simulations. The straight line denotes a regression fit. The black diamond denotes the exact result [7].

Figure 5. The estimate of the wetting surface field for the 2D system at [Forumla omitted. See PDF.] Circles correspond to the surface critical fields obtained from simulations. The straight line denotes a regression fit. The black diamond denotes the exact result [7].

Figure 6. The effective interface potentials calculated for the 3D non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for surface fields from [Forumla omitted. See PDF.], −0.891, −0.892, −0.893, −0.894, −0.895, −0.896, −0.898, −0.90, −0.902, −0.904, −0.906, and −0.908. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 7. The effective interface potentials calculated for the 3D non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for surface fields [Forumla omitted. See PDF.], −0.891, −0.892, −0.893, −0.894, −0.895, −0.896, −0.898, −0.90, −0.902, −0.904, −0.906, and −0.908. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 8. The effective interface potentials calculated for the 3D non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for surface fields [Forumla omitted. See PDF.], −0.891, −0.892, −0.893, −0.894, −0.895, −0.896, −0.898, −0.90, −0.902, −0.904, −0.906, and −0.908. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 9. The effective interface potentials calculated for the non-symmetric system at [Forumla omitted. See PDF.] and for [Forumla omitted. See PDF.]. The potentials are calculated for surface fields from top [Forumla omitted. See PDF.], −0.557, −0.558, −0.559, −0.560, −0.561, −0.562, −0.563, −0.564, −0.565, −0.568, −0.570, −0.575, −0.582, −0.588, −0.592, −0.596, −0.600, −0.604, and −0.608. The thick line corresponding to [Forumla omitted. See PDF.] denotes the effective interface potential with no detectable local minimum.

Figure 10. The integrand of Equation (16) vs. the surface field for [Forumla omitted. See PDF.] and for three system sizes L listed in Figure. The inset shows zoom-out of the main Figure.

Figure 11. The integrand of Equation (16) vs. the surface field for [Forumla omitted. See PDF.] and for three system sizes L listed in Figure. The inset shows zoom-out of the main Figure.

Figure 12. [Forumla omitted. See PDF.] vs. the surface field for [Forumla omitted. See PDF.] and for three system sizes L listed in Figure. The inset shows full dependence of [Forumla omitted. See PDF.] vs. [Forumla omitted. See PDF.].

Figure 13. [Forumla omitted. See PDF.] vs. the surface field for [Forumla omitted. See PDF.] and for three system sizes L listed in Figure. The inset shows full dependence of [Forumla omitted. See PDF.] vs. [Forumla omitted. See PDF.].

Figure 14. Symbols denote the estimate of the surface field for which [Forumla omitted. See PDF.] at a given L. Thick lines denote the linear regression extrapolating [Forumla omitted. See PDF.]. Panel (a) is for [Forumla omitted. See PDF.] while panel (b) is for [Forumla omitted. See PDF.].

Figure 15. [Forumla omitted. See PDF.] for five different system sizes L, evaluated at [Forumla omitted. See PDF.] for the symmetric surface fields, [Forumla omitted. See PDF.]. The big black circle with error bars denotes the simulational result of Binder et al. [15].

Figure 16. Estimation of [Forumla omitted. See PDF.] from simulational data presented in Figure 15, evaluated at [Forumla omitted. See PDF.] for the symmetric surface fields.

Figure 17. [Forumla omitted. See PDF.] for different system sizes, evaluated at [Forumla omitted. See PDF.] for the symmetric surface fields.

Figure 18. Estimation of [Forumla omitted. See PDF.] from simulational data presented in Figure 17, evaluated at [Forumla omitted. See PDF.] for the symmetric surface fields.

Table 1

Critical surface field $$H_{1c}/J$$ obtained from various methods.

$$\mathit{J}/{\mathit{k}}_{\mathit{B}}\mathit{T}$$	AFSS ([18])	EIP (This Work)	TIN (This Work)	BLK (This Work)	BLK ([13,15])
0.25	−0.616	−0.608	−0.604	−0.606	−0.555
0.35	−0.909	−0.908	−0.905	−0.905	−0.89

References

1. Saha, R.; Uppaluri, R.V.S.; Tiwari, P. Effects of interfacial tension, oil layer break time, emulsification and wettability alteration on oil recovery for carbonate reservoirs. Colloid Surf. A; 2018; 559, pp. 92-103. [DOI: https://dx.doi.org/10.1016/j.colsurfa.2018.09.045]

2. Rudolph, M.; Hartmann, R. Specific surface free energy component distributions and flotabilities of mineral microparticles in flotation—An inverse gas chromatography study. Colloid Surf. A; 2017; 513, pp. 380-388. [DOI: https://dx.doi.org/10.1016/j.colsurfa.2016.10.069]

3. Chu, Z.; Seeger, S. Superamphiphobic surfaces. Chem. Soc. Rev.; 2014; 43, pp. 2784-2798. [DOI: https://dx.doi.org/10.1039/C3CS60415B] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24480921]

4. Dietrich, S. Wetting Phenomena. Phase Transitions and Critical Phenomena; Domb, C.; Lebowitz, J.L. Academic: London, UK, 1988; Volume 12, pp. 1-54.

5. Fisher, D.S.; Huse, D.A. Wetting transitions: A functional renormalization-group approach. Phys. Rev. B; 1986; 32, pp. 247-256. [DOI: https://dx.doi.org/10.1103/PhysRevB.32.247] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9936660]

6. Dietrich, S.; Napiórkowski, M. Analytic results for wetting transitions in the presence of van der Waals tails. Phys. Rev. A; 1991; 43, pp. 1861-1885. [DOI: https://dx.doi.org/10.1103/PhysRevA.43.1861]

7. Abraham, D.B. Solvable Model with a Roughening Transition for a Planar Ising Ferromagnet. Phys. Rev. Lett.; 1980; 44, pp. 1165-1168. [DOI: https://dx.doi.org/10.1103/PhysRevLett.44.1165]

8. Brezin, E.; Halperin, B.I.; Leibler, S. Critical Wetting in Three Dimensions. Phys. Rev. Lett.; 1983; 50, pp. 1387-1390. [DOI: https://dx.doi.org/10.1103/PhysRevLett.50.1387]

9. Nakanishi, H. Fisher, M.E. Multicriticality of Wetting, Prewetting, and Surface Transitions. Phys. Rev. Lett.; 1982; 49, pp. 1565-1568. [DOI: https://dx.doi.org/10.1103/PhysRevLett.49.1565]

10. Taborek, P.; Rutledge, J.E. Novel wetting behavior of 4He on cesium. Phys. Rev. Lett.; 1992; 68, pp. 2184-2187. [DOI: https://dx.doi.org/10.1103/PhysRevLett.68.2184]

11. Friedman, S.R.; Khalil, M.; Taborek, P. Wetting Transition in Water. Phys. Rev. Lett.; 2013; 111, 226101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.111.226101]

12. Napiórkowski, M.; Dietrich, S. Wetting Transitions in Terms of Effective Potentials. Phys. Rev. Lett.; 2015; 114, 039601. [DOI: https://dx.doi.org/10.1103/PhysRevLett.114.039601]

13. Binder, K.; Landau, D.P.; Kroll, D.M. Critical Wetting with Short-Range Forces: Is Mean-Field Theory Valid?. Phys. Rev. Lett.; 1986; 56, pp. 2272-2275. [DOI: https://dx.doi.org/10.1103/PhysRevLett.56.2272]

14. Binder, K.; Landau, D.P. Wetting and layering in the nearest-neighbor simple-cubic Ising lattice: A Monte Carlo investigation. Phys. Rev. B; 1988; 37, pp. 1745-1765. [DOI: https://dx.doi.org/10.1103/PhysRevB.37.1745]

15. Binder, K.; Landau, D.P.; Wansleben, S. Wetting transitions near the bulk critical point: Monte Carlo simulations for the Ising model. Phys. Rev. B; 1989; 40, pp. 6971-6979. [DOI: https://dx.doi.org/10.1103/PhysRevB.40.6971] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9991076]

16. Parry, A.O.; Rascon, C.; Bernardino, N.R.; Romero-Enrique, J.M. 3D Short-Range Wetting and Nonlocality. Phys. Rev. Lett.; 2008; 100, 136105. [DOI: https://dx.doi.org/10.1103/PhysRevLett.100.136105] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18517973]

17. Albano, E.V.; Binder, K. Finite-Size Scaling Approach for Critical Wetting: Rationalization in Terms of a Bulk Transition with an Order Parameter Exponent Equal to Zero. Phys. Rev. Lett.; 2012; 109, 036101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.109.036101]

18. Bryk, P.; Binder, K. Non-mean-field behavior of critical wetting transition for short-range forces. Phys. Rev. E; 2013; 88, 030401. [DOI: https://dx.doi.org/10.1103/PhysRevE.88.030401] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24125203]

19. Wansleben, S. Ultrafast vectorized multispin coding algorithm for the Monte Carlo simulation of the 3D Ising model. Comput. Phys. Commun.; 1987; 43, pp. 315-323. [DOI: https://dx.doi.org/10.1016/0010-4655(87)90048-8]

20. Yan, Q.; de Pablo, J.J. Hyper-parallel tempering Monte Carlo: Application to the Lennard-Jones fluid and the restricted primitive model. J. Chem. Phys.; 1999; 111, pp. 9505-9516. [DOI: https://dx.doi.org/10.1063/1.480282]

21. MacDowell, L.G.; Müller, M. Observation of autophobic dewetting on polymer brushes from computer simulation. J. Phys. Condens. Matter; 2005; 17, pp. S3523-S3528. [DOI: https://dx.doi.org/10.1088/0953-8984/17/45/043]

22. MacDowell, L.G.; Müller, M. Adsorption of polymers on a brush: Tuning the order of the wetting phase transition. J. Chem. Phys.; 2006; 124, 084907. [DOI: https://dx.doi.org/10.1063/1.2172597] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16512742]

23. Benet, J.; Palanco, J.G.; Sanz, E.; MacDowell, L.G. Disjoining Pressure, Healing Distance, and Film Height Dependent Surface Tension of Thin Wetting Films. J. Phys. Chem. C; 2014; 118, pp. 22079-22089. [DOI: https://dx.doi.org/10.1021/jp506534b]

24. Evans, R.; Wilding, N.B. Quantifying Density Fluctuations in Water at a Hydrophobic Surface: Evidence for Critical Drying. Phys. Rev. Lett.; 2015; 115, 016103. [DOI: https://dx.doi.org/10.1103/PhysRevLett.115.016103]

25. Rane, K.S.; Kumar, V.; Errington, J.R. Monte Carlo simulation methods for computing the wetting and drying properties of model systems. J. Chem. Phys.; 2011; 135, 234102. [DOI: https://dx.doi.org/10.1063/1.3668137] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22191859]

26. Jain, K.; Schultz, A.J.; Errington, J.R. Application of the interface potential approach for studying wetting behavior within a molecular dynamics framework. J. Chem. Phys.; 2019; 150, 204118. [DOI: https://dx.doi.org/10.1063/1.5096362] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31153165]

27. Jain, K.; Schultz, A.J.; Errington, J.R. Construction of the interface potential froma series of canonical ensemble simulations. J. Chem. Phys.; 2019; 151, 044103. [DOI: https://dx.doi.org/10.1063/1.5110922] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31370560]

28. Jain, K.; Rane, K.S.; Errington, J.R. Using isothermal-isobaric Monte Carlo simulation to study the wetting behavior of model systems. J. Chem. Phys.; 2019; 150, 084110. [DOI: https://dx.doi.org/10.1063/1.5089416]

29. Hughes, A.P.; Thiele, U.; Archer, A.J. Liquid drops on a surface: Using density functional theory to calculate the binding potential and drop profiles and comparing with results from mesoscopic modelling. J. Chem. Phys.; 2015; 142, 074702. [DOI: https://dx.doi.org/10.1063/1.4907732]

30. Hughes, A.P.; Thiele, U.; Archer, A.J. Influence of the fluid structure on the binding potential: Comparing liquid drop profiles from density functional theory with results from mesoscopic theory. J. Chem. Phys.; 2017; 146, 064705. [DOI: https://dx.doi.org/10.1063/1.4974832]

31. Winter, D.; Virnau, P.; Binder, K. Monte Carlo Test of the Classical Theory for Heterogeneous Nucleation Barriers. Phys. Rev. Lett.; 2009; 103, 225703. [DOI: https://dx.doi.org/10.1103/PhysRevLett.103.225703]

32. Winter, D.; Virnau, P.; Binder, K. Heterogeneous nucleation at a wall near a wetting transition: A Monte Carlo test of the classical theory. J. Phys. Condens. Matter; 2009; 21, 464118. [DOI: https://dx.doi.org/10.1088/0953-8984/21/46/464118]

33. Hasenbusch, M.; Pinn, K. Surface tension, surface stiffness, and surface width of the 3-dimensional Ising model on a cubic lattice. Physica A; 1993; 192, pp. 342-374. [DOI: https://dx.doi.org/10.1016/0378-4371(93)90043-4]

34. Hasenbusch, M.; Pinn, K. Comparison of Monte Carlo results for the 3D Ising interface tension and interface energy with (extrapolated) series expansions. Physica A; 1994; 203, pp. 189-213. [DOI: https://dx.doi.org/10.1016/0378-4371(94)90152-X]

35. Binder, K.; Landau, D.P. Finite-size scaling at first-order phase transitions. Phys. Rev. B; 1984; 30, pp. 1477-1485. [DOI: https://dx.doi.org/10.1103/PhysRevB.30.1477]

36. Binder, K. Institut für Physik, Johannes Gutenberg Universität, Mainz, Germany. Private Communication; 2011.

37. Evans, R.; Stewart, M.C.; Wilding, N.B. A unified description of hydrophilic and superhydrophobic surfaces in terms of the wetting and drying transitions of liquids. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 23901-23908. [DOI: https://dx.doi.org/10.1073/pnas.1913587116] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31611388]