Chasing the Critical Wetting Transition. An Effective Interface Potential Method

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.


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 ω = k B T 4πΣξ 2 (1) where T is the temperature, ξ is the correlation length of the bulk wetting phase and Σ is the interfacial stiffness. Capillary-wave-like fluctuations give rise to a diverging transverse correlation length ξ and the relevant exponent is believed to behave for ω > 2 . (2) 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 ξ NL = lξ ∝ ln ξ . This gives rise to an effective wetting parameter, ω e f f , of the form where κ = ξ −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.

Materials and Methods
Our model consists of cubic lattice of dimensions L × L × D with two free boundary layers L × 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 = ±1. The Hamiltonian for the system is 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 = 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 where 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 total magnetization and "mixed" surface layer susceptibility

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(m p ) of finding a system with magnetization m p defined as where D r denotes the range of interest, D/2 ≤ D r ≤ 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 In our work, we used P(m p ), while MacDowell and Müller [21,22] calculated the adsorption distribution P(Γ). 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 e f f (m p ) differs from that used traditionally in the literature [4,5,8]. Those authors considered V e f f (l) where l is a distance of the interface from the wall, a local quantity. V e f f (m p ) 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]. [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]

MacDowell and Müller (and others
These calculations provide a clear-cut, stringent test of the proposed approach. Figures 1-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 → ∞. We note, that the simulational value of H 1c (∞)/J = −0.935 agrees very well with the exact result, −0.927. Another important observation is that the critical surface field H 1c (L) 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 + √ 2)/2 ≈ 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(m p ) 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. Figures 6-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 ± 0.002. Interestingly, for the larger systems (cf. Figures 7 and 8) the transition is located also at H 1 /J = −0.908 ± 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 ± 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].

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 where γ 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, where α s is the surface critical exponent for specific heat, t = 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 ∆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. Figures 10 and 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 Figures 12 and 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 Figures 12 and 13). The derivative should be zero at the critical surface field, H 1c (L). In order arrive at an estimate of H 1c (L) we plotted two parallel lines cos θ = c 1 and cos θ = c 2 and constructed tangential lines via finite difference. Extrapolation of these straight lines give H 1 (L) satisfying cos θ = 0. Finally from the plots of H 1 (L) vs. 1/L the thermodynamic limit is reached via extrapolation (cf. Figure 14).

.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 χ 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 × 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 ± 0.0005 agrees very well with earlier study (H 1c (L = 128)/J = −0.89 ± 0.004). However, there is a considerable shift for larger systems. This indicates that H 1c (L = ∞) is lower. This is confirmed in Figure 16, where we show extrapolation L → ∞, and find, that the estimated value H 1c (L = ∞)/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 = ∞ reveals (cf. Figure 18), that that the estimated value H 1c (L = ∞)/J = −0.606. This is again, noticeably lower than previous estimates calculated using much smaller system sizes (H 1c (L = ∞)/J = −0.555 [13,14]).
As the careful Reader noticed, there is one vexing feature of the results presented in Figures 15 and 17. Namely the maxima of χ 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 χ 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 = ∞)/J = −0.905 for J/k B T = 0. 35

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.

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.