Drying and Wetting Transitions of a Lennard-Jones Fluid: Simulations and Density Functional Theory

TL;DR

This study combines simulations and density functional theory to analyze how Lennard-Jones fluids undergo drying and wetting transitions at flat walls, revealing the effects of interaction range on transition order and critical behavior.

Contribution

It provides a detailed comparison of drying and wetting transitions for long-ranged and short-ranged wall-fluid interactions using simulations and theory, highlighting differences in transition order and critical exponents.

Findings

Wetting is first order for long-range interactions; drying is critical and occurs at zero wall attraction.

Drying is critical in the short-range case; wetting transition order depends on potential truncation.

Estimated correlation length exponent exceeds mean field predictions, indicating unexpected critical behavior.

Abstract

We report a theoretical and simulation study of the drying and wetting phase transitions of a truncated Lennard-Jones fluid at a flat structureless wall. Binding potential calculations predict that the nature of these transitions depends on whether the wall-fluid attraction has a long ranged (LR) power law decay, or is instead truncated, rendering it short ranged (SR). Using grand canonical Monte Carlo simulation and classical density functional theory we examine both cases in detail. We find that for the LR case wetting is first order, while drying is continuous (critical) and occurs exactly at zero attractive wall strength, ie. in the limit of a hard wall. In the SR case, drying is also critical but the order of the wetting transition depends on the truncation range of the wall-fluid potential. We characterize the approach to critical drying and wetting in terms of the density and local compressibility profiles and via the finite-size scaling properties of the probability distribution of the overall density. For the LR case, where the drying point is known exactly, this analysis allows us to estimate the exponent $\nu_\parallel$ which controls the parallel correlation length, i.e. the extent of vapor bubbles at the wall. Surprisingly, the value we obtain is over twice that predicted by mean field and renormalization group calculations, despite the fact that our three dimensional system is at the upper critical dimension where mean field theory for critical exponents is expected to hold. We suggest reasons for this discrepancy.

Figures

Equations

$${\gamma }_{vl}\cos (\theta )={\gamma }_{wv}-{\gamma }_{wl},$$

$$\phi_{\rm att}(r)=\left\{\begin{array}[]{ll}4\epsilon_{LJ}\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],&r\leq r_{c}\,,\\ 0,&r>r_{c},\\ \end{array}\right.$$

$$V=(L\sigma {)}^{2}D\sigma .$$

$$W_{\rm SR}(z)=\left\{\begin{array}[]{ll}\infty,\mbox{\hskip 11.38109pt}&z\leq 0\\ -\epsilon,&0<z<\sigma/2\,,\\ 0,&z>\sigma/2,\\ \end{array}\right.$$

$$W_{\rm LR}(z)=\left\{\begin{array}[]{ll}\infty,\mbox{\hskip 11.38109pt}&z\leq 0\\ \epsilon_{w}\epsilon_{LJ}\left[\frac{2}{15}\left(\frac{\sigma}{z}\right)^{9}-\left(\frac{\sigma}{z}\right)^{3}\right],&z>0\,,\\ \end{array}\right.$$

$$\rho (\mathbf{r})=⟨\hat{\rho }(\mathbf{r})⟩\equiv ⟨\sum _{i=1}^N\delta (\mathbf{r}-{\mathbf{r}}_i)⟩$$

$$G({\mathbf{r}}_1,{\mathbf{r}}_2)$$

$$\Gamma \equiv {\int }_0^{\infty }dz(\rho (z)-{\rho }_b)=-\frac{1}{A}{\left(\frac{\partial {\Omega }_{ex}}{\partial \mu }\right)}_T$$

$${\chi }_{ex}\equiv {\left(\frac{\partial \Gamma }{\partial \mu }\right)}_T=\frac{\beta }{A}[(⟨N^2⟩-⟨N{⟩}^2)-(⟨N_b^2⟩-⟨N_b{⟩}^2)]$$

$$\frac{1}{A}d({\Omega }_{ex})=-sdT-\Gamma d\mu -\Xi dϵ,$$

$$\Theta =-\frac{1}{A}{\left(\frac{\partial {\Omega }_{ex}}{\partial ϵ}\right)}_{\mu ,T}=-{\int }_0^{\infty }dz\rho (z)\frac{W_{ϵ}(z)}{ϵ}$$

$${\Gamma }_1\equiv {\left(\frac{\partial \Gamma }{\partial ϵ}\right)}_T={\left(\frac{\partial \Theta }{\partial \mu }\right)}_T\equiv {\chi }_1$$

$${\left(\frac{\partial m_s}{\partial h_1}\right)}_T={\left(\frac{\partial m_1}{\partial h}\right)}_T$$

$$\chi (z)\equiv {\left(\frac{\partial \rho (z)}{\partial \mu }\right)}_T$$

$$S(z_1;q)\equiv {\int }_{-\infty }^{\infty }d{z}_2\int d\mathbf{R}{e}^{i\mathbf{q}\cdot \mathbf{R}}G(z_1,z_2;R)$$

$$\chi (z_1)=\beta S(z_1;0)=\beta {\int }_{-\infty }^{\infty }d{z}_2\int d\mathbf{R}G(z_1,z_2;R)$$

$$\chi (\mathbf{r})={\beta }^{-1}⟨N\hat{\rho }(\mathbf{r})-⟨N⟩⟨\hat{\rho }(\mathbf{r})⟩⟩$$

$${\chi }_{ex}={\int }_0^{\infty }dz[\chi (z)-{\chi }_b]={\int }_0^{\infty }dz[\beta S(z;0)-{\chi }_b]$$

$$|\Gamma |\sim (\delta {ϵ}_w{)}^{-{\beta }_s},\delta {ϵ}_w\to 0$$

$${\xi }_{\parallel }\sim (\delta {ϵ}_w{)}^{-{\nu }_{\parallel }},\delta {ϵ}_w\to 0$$

$$\chi (z)=\beta S(z;0)\sim {\rho }^{\prime }(z){\xi }_{\parallel }^2,z\approx l$$

$${\chi }_{ex}\sim {\xi }_{\parallel }^2\sim (\delta {ϵ}_w{)}^{-2{\nu }_{\parallel }}.$$

$${\chi }_1\sim (\delta {ϵ}_w{)}^{-{\beta }_s-1},\delta {ϵ}_w\to 0$$

$$|{\gamma }^{sing}|\sim (\delta {ϵ}_w{)}^{2-{\alpha }_s}$$

$$2-{\alpha }_s=2({\nu }_{\parallel }-{\beta }_s).$$

$$\delta {\mu }_{ce}\approx \frac{-2{\gamma }_{lv}\cos (\theta )}{D({\rho }_l-{\rho }_v)}$$

$${\omega }^{ex}(l)={\gamma }_{wv}+{\gamma }_{lv}+{\omega }_B(l)+\delta \mu ({\rho }_l-{\rho }_v)l$$

$${\omega }_B(l)=a\exp (-l/{\xi }_b)+b{l}^{-2}+H.O.T.$$

$$b=-({\rho }_l-{\rho }_v){ϵ}_w{ϵ}_{LJ}{\sigma }^3/2$$

$$\frac{-l_{eq}}{{\xi }_b}=\ln {ϵ}_w-3\ln \left(\frac{l_{eq}}{{\xi }_b}\right)+\mathrm{constants};\delta \mu =0^{+}$$

Full text

I Introduction

The behaviour of a liquid drop in equilibrium with its vapor and in contact with a flat substrate (or ‘wall’) is characterised in thermodynamic terms by the contact angle $\theta$ that the drop makes with the substrate de Gennes (1985). The precise value of $\theta$ depends on the surface chemistry of the substrate, but in broad terms, strong wall-fluid attraction is associated with a small contact angle, while weak attraction is associated with a large contact angle, as shown schematically in Fig. 1. On increasing the wall-fluid attraction the contact angle approaches the limit $\theta\to 0^{\circ}$. This corresponds to the wetting transition in which a macroscopic liquid layer intrudes between the wall and the vapor. The drying transition is the counterpart of wetting that occurs as the wall attraction is progressively weakened so that $\theta\to 180^{\circ}$, whereupon a macroscopic vapor layer intrudes between the wall and the liquid. Values of $0^{\circ}<\theta<90^{\circ}$ are often termed partially wet, while values $90^{\circ}<\theta<180^{\circ}$ are termed partially dry. Young’s equation:

[TABLE]

expresses $\theta$ in terms of the wall-vapor ($wv$), wall-liquid ($wl$) and vapor-liquid ($vl$) surface tensions.

Let us consider first the nature of wetting. This has been the subject of enduring experimental, theoretical and simulation interest Bonn et al. (2009). As is well established, the character of the transition can be either discontinuous (first order) or continuous (critical) depending in a subtle fashion on whether the wall-fluid ($wf$) potential is long-ranged (LR) or short ranged (SR) and on whether the fluid-fluid ($ff$) interaction is LR or SR. Nightingale et al. (1983, 1984); Dietrich and Schick (1985); Ebner et al. (1985); Ebner and Saam (1987a, b); Binder et al. (1986); Ross et al. (2001). Most experimental studies of wetting transitions find these to be first order (see Bonn and Ross (2001) for a review and Friedman et al. (2013) for a recent study of water). Critical wetting is much rarer (although potentially more interesting from a fundamental perspective) and experimental reports have, to date, been limited to a few special cases Ragil et al. (1996); Ross et al. (1999, 2001). Consequently, the bulk of progress towards elucidating the character of critical wetting has come from theoretical and simulation approaches, which have revealed a wealth of complex behaviour  Nakanishi and Fisher (1982); Brézin et al. (1983); Nightingale et al. (1983, 1984); Dietrich and Schick (1985); Ebner et al. (1985); Ebner and Saam (1987b); Binder et al. (1986, 1989); Ragil et al. (1996); Ross et al. (1999); Evans and Parry (1990, 1989); Parry et al. (1991, 2008a, 2008b); Albano and Binder (2012); Bryk and Binder (2013). Theoretical treatments typically adopt a mean field (MF) approach based on Landau theory, binding potentials or density functional theory (DFT). These predict that critical wetting occurs when both $wf$ and $ff$ interactions are SR, and that $d=3$ is the upper critical dimension. Hence three dimensional systems are a borderline case for the validity of mean field theory. Renormalization Group (RG) calculations Brézin et al. (1983); Fisher and Huse (1985) predict non-universal behaviour in this instance– a finding that has prompted concerted simulation efforts to clarify the nature of the criticality. Unfortunately simulations are complicated by the finite-size effects that stem from the divergent critical correlations. Early work for a nearest neighbor Ising model with surface fields Binder et al. (1986, 1989) yielded pronounced discrepancies with RG predictions regarding the critical exponents, prompting efforts to extend the theoretical framework to non-local interfacial Hamiltonians Parry et al. (2008b, a); Parry and Rascón (2009) in order to explain the differences. More recent simulation studies have attempted to deploy finite-size scaling (FSS) techniques, using data from a wide range of lattice sizes, in order to pinpoint the transition and elucidate its character Albano and Binder (2012); Bryk and Binder (2013). Significant inaccuracies were identified in the earlier estimates of the transition point, which had in turn skewed the estimates of critical point properties. However, difficulties in applying the FSS methodology to three dimensional systems were also reported Bryk and Binder (2013) and thus the precise relationship between simulation data and theoretical predictions is arguably not settled fully.

As indicated above, simulation studies of critical wetting have been confined exclusively to Ising models on account of their computational tractability. However, it is interesting to ask to what extent realistic fluid models share the properties of the lattice systems. Although fluid-magnet universality suggests that the critical scaling properties should be identical, the particle-hole symmetry inherent in Ising models is expected to engender features that do not occur in real fluids. For example wetting is formally equivalent to drying for SR surface magnetic fields in the Ising model, but there is no reason to think that wetting and drying are equivalent for realistic fluid models. Furthermore it is unclear how to translate some of the model parameters commonly employed in the Ising context, such as enhanced surface layer couplings and surface magnetic fields, to the case of realistic fluids in which the substrate-fluid interactions are fully described by the associated wall potential. What is clear, however, is that simulations of realistic fluids can be expected to be considerably more computationally demanding than for Ising models, and accordingly one should expect that the range of accessible system sizes is correspondingly smaller.

In common with wetting, relatively little is known about the fundamental nature of drying transitions in realistic fluids. Experimental studies are far scarcer than for wetting owing to the challenges of fabricating substrates for which the contact angle of a liquid drop is large. For instance, for water on ‘hydrophobic’ surfaces such as Teflon, wax or self-assembled monolayers, the contact angle does not generally exceed $\theta\approx 130^{\circ}$. While this precludes detailed study of the approach to the drying transition, interesting effects have nevertheless been reported for strongly hydrophobic substrates, such as a depletion region of one or two molecular layers adjacent to the substrate in which the one-body density is considerably reduced compared to its bulk value Mezger et al. (2006); Ocko et al. (2008); Mezger et al. (2010); Chattopadhyay et al. (2010, 2011); Uysal et al. (2013). Encouragingly, the ability to study the drying transition in detail is likely to improve in the future with the advent of novel nano and micro-structured surfaces which exhibit contact angles approaching $170^{\circ}$ Lafuma and Quere (2003); Quéré (2005); Simpson et al. (2015); Li et al. (2007); Ueda and Levkin (2013). Such ‘superhydrophobic’ surfaces are of widespread interest for their potential technological applications, including self cleaning surfaces and chemical separation processes Simpson et al. (2015).

On the theoretical side, it is well established that complete drying occurs for any liquid (that exhibits liquid-vapor coexistence) adsorbed at a planar hard wall, see e.g. Henderson and van Swol (1985); Oettel (2005). However, the situation for attractive $wf$ interactions is less clear. Early work emphasized the central role of the range of the potentials Ebner and Saam (1987a, b) and reported that for a SR lattice-gas system with LR $wf$ interactions, no drying transition can occur, while any wetting transition is first order. This important prediction is specific to the lattice-gas (Ising) model. We reexamine this in the context of a fluid in contact with a hard wall plus attractive LR tail.

On the simulation front, a number of studies have been performed, but no clear consensus regarding the nature of drying has yet emerged. Monte Carlo simulation studies of a Lennard-Jones liquid Oleinikova et al. (2005); Brovchenko et al. (2005, 2004) utilizing LR $wf$ interactions reported no signs of a drying transition, stating this was in accord with the theoretical predictions Ebner and Saam (1987b). Molecular Dynamics (MD) studies by two separate groups, of both truncated Lennard-Jones (LJ) and square-well fluids reported a drying transition at some small but non-zero strength of SR $wf$ attraction, but disagreed regarding its character, with one group claiming that the transition is first order Henderson and van Swol (1990); Henderson et al. (1992); van Swol and Henderson (1989, 1991) and the other Nijmeijer et al. (1992, 1991); Bruin et al. (1995); Bruin (1998) that it is critical. (Note that DFT calculations came out in favor of critical drying Henderson et al. (1992).) More recently Monte Carlo (MC) studies of the contact angle in a LJ fluid Rane et al. (2011) and SPC/E water Kumar and Errington (2013a) at various types of substrate pointed tentatively to a critical drying transition for the LR $wf$ case. Hints of differences in the character of the approach to drying between systems with LR and SR wall-fluid interactions were also noted. Separate MD studies of a model for water at a weakly attractive substrate found evidence for enhanced density fluctuations in the surface region Chandler (2007); Patel et al. (2010, 2012); Willard and Chandler (2014); Mittal and Hummer (2010); Jamadagni et al. (2011). The latter finding was subsequently rationalized by MC simulations for SPC/E water which provided firm evidence that the drying transition in water is critical and that this fact is responsible for the enhanced fluctuations in the surface region Evans and Wilding (2015).

Despite substantial progress made in understanding the physics of wetting and drying in realistic fluids, a number of fundamental and practical questions remain unanswered. Principal among these are: i) What characteristics of the model system are responsible for determining the order of wetting and drying transitions? ii) How can wetting and drying points be located accurately via computer simulation? iii) If a transition is critical, what is the nature of the near-critical fluctuations and the values of the critical exponents and how can these be measured accurately? iv) What is the role of finite-size effects in characterizing surface phase transitions?

In the present contribution we address these issues using a combination of theoretical and computational techniques. For simulational expediency we consider $ff$ interactions that are exclusively SR in nature (a truncated LJ potential), but for the $wf$ interactions we consider both the LR and SR cases. Our focus is mainly on the drying transition in systems with LR $wf$ interactions because: (a) drying turns out to be critical for LR $wf$; (b) the drying point occurs at exactly zero attractive wall strength (ie in the limit of a hard wall) – a feature which allows us to study the transition free from uncertainty regarding its location. However, we also report results for drying in the case when the $wf$ interactions are SR and for wetting in the case of LR and SR $wf$ interactions. The latter case of wetting for SR $wf$ is also found to be critical but the theoretical predictions for the nature of the criticality are different from those for drying with LR $wf$ interactions.

The layout of the paper is as follows 111A short account of some of our results has previously been published elsewhere Evans et al. (2016). Section II describes our model of a Lennard-Jones fluid confined between smooth parallel walls. In Section III we set out some general features of adsorption, surface phase behavior and criticality for a simple fluid. Section. IV provides details of the theoretical methods we have used to study the model, namely MF and RG analysis of a binding potential and classical DFT. The grand canonical Monte Carlo (GCMC) simulation methods are described in Section VI. Our results for the drying and wetting properties are set out in Sections V and VII, and are discussed in Section VIII.

II Model Fluids

Our grand canonical Monte Carlo (GCMC) simulations consider the drying and wetting behaviour of a Lennard-Jones (LJ) fluid in which particles interact via the potential,

[TABLE]

with $\epsilon_{LJ}$ the well-depth of the potential and $\sigma$ the LJ diameter. We choose $r_{c}=2.5\sigma$, for which criticality occurs Wilding (1995) at $k_{B}T_{c}=1.1876(3)\epsilon_{LJ}$. We work at $k_{B}T=0.91954\epsilon_{LJ}=0.775T_{c}$ for which coexistence occurs at $\beta\mu_{co}=-3.865950(20)$, with coexistence densities $\rho_{l}\sigma^{3}=0.704(1)$ and $\rho_{v}\sigma^{3}=0.0286(2)$; and also at $k_{B}T=1.0\epsilon_{LJ}=0.842T_{c}$ for which $\beta\mu_{co}=-3.457131(25)$, $\rho_{l}\sigma^{3}=0.653(1),\rho_{v}\sigma^{3}=0.0504(3)$, $\beta=(k_{B}T)^{-1}$. The choice of cutoff $r_{c}$ is motivated computationally and follows that of most of the LJ community.

The fluid is confined within a slit pore comprising two planar walls of area $(L\sigma)^{2}$ separated by a distance $D\sigma$, so that the volume is

[TABLE]

Periodic boundary conditions are applied in the directions parallel to the planar walls. We employ three types of wall-fluid potential in our GCMC simulations. The SR potential for a single wall is a square-well given by

[TABLE]

where $\epsilon$ is the well-depth.

Two types of LR potential are considered in this work. The first is given by the well known $9$-$3$ form having a decaying attractive part at large $z$ and a steep repulsive part at small $z$:

[TABLE]

where $\epsilon_{w}$ is a dimensionless measure of the strength of the wall-fluid attraction. At the minimum of (5) the value of the wall-fluid potential is $-1.0541\epsilon_{w}\epsilon_{LJ}=-\epsilon$.

The second LR potential, a modification of the first, is obtained by making the replacement $z\to{\tilde{z}}$ in (5) with ${\tilde{z}}=z+(2/5)^{1/6}\sigma$. This shifts the minimum of the $9$-$3$ potential to the hard wall at $z=0$, leading to an infinitely steep repulsive part. The motivation for utilizing the modified form is that when studying a slit geometry (see below) the wall separation and hence the slit volume is unambiguously defined for all $\epsilon_{w}$. This is not the case for the standard form (5) in the regime of interest for drying, namely $\epsilon_{w}\approx 0$. As $\epsilon_{w}\to 0$ the effective wall position (given by the value of $z$ for which the repulsive energy is $\sim k_{B}T$) shifts strongly, leading to artifacts in measurements of the total number density. The GCMC results reported in Sec. VII are for the modified potential. The DFT results in Sec. V are for (5).

III Background to Wetting, Drying and Confinement with Planar Walls

III.1 Statistical Mechanics of Adsorption: Thermodynamics and Correlation Functions

In this subsection we summarize key results in the statistical mechanics of adsorption pertinent to our GCMC and DFT investigations. For convenience we consider the fluid to be adsorbed at a single (planar) wall of (infinite) interfacial area $A$ exerting a potential of the type (4,5). Then the average one-body density $\rho({\bf r})=\rho(z)=0,z<0$. Extension to the fluid confined by two walls is straightforward. As is appropriate for adsorption studies, we work grand canonically with a reservoir at fixed chemical potential $\mu$ and temperature $T$. Thus

[TABLE]

where the brackets $\langle\rangle$ denote a GC average and we introduced the usual particle density operator for $N$ particles with coordinates ${\bf r}_{i}$. Two-body correlations are described by the density-density correlation function:

[TABLE]

where $R=\sqrt{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}}$ is the transverse separation between particles (atoms). Connection with surface thermodynamics is made via the Gibbs adsorption equation

[TABLE]

i.e. $\Gamma$, the excess number of particles per unit area, is the minus of the derivative of the excess grand potential $\Omega_{ex}$ w.r.t. $\mu$. $\rho_{b}\equiv\rho_{b}(\mu,T)$ is the density of the bulk fluid far from the wall and the excess quantity is defined by $\Omega_{ex}\equiv\Omega+pV$ where $\Omega$ is the total grand potential, $V$ is the accessible volume and $p\equiv p(\mu,T)$ is the pressure of the bulk fluid.

The second derivative of $\Omega_{ex}$ w.r.t. $\mu$ yields the surface excess compressibility, and this can be written as a fluctuation formula Evans and Marconi (1987); Evans and Stewart (2015)

[TABLE]

where it is implied that the surface area $A\to\infty$. The first term in (9) is the mean-square fluctuation in the total number of particles, which must be positive to ensure stability. The second is the corresponding quantity for the bulk fluid at the same $(\mu,T)$. The difference can be negative Evans and Marconi (1987). Bratko et al Bratko et al. (2007); Bratko (2010) measured the first term of (9) and, more recently, Kumar and Errington Kumar and Errington (2013b) measured $\chi_{ex}$ in GCMC simulations of SPC/E water at hydrophobic substrates. Whilst $\chi_{ex}$ provides a measure of the overall compressibility of the adsorbed fluid and the strength of fluctuations in the total number of particles, this quantity does not provide information about the spatial location of the important density fluctuations, i.e. at which distances $z$ from the wall these are most pronounced. Below we define the local compressibility $\chi(z)$ which does provide the appropriate measure.

First we introduce a further sum rule relating a thermodynamic quantity to an (integrated) microscopic quantity. Following Henderson (1986); Evans and Parry (1989) we suppose that $W(z)$ is such that $\partial W(z)/\partial\epsilon\equiv W_{\epsilon}(z)/\epsilon$ is independent of the well-depth $\epsilon$. (The wall-fluid potentials we consider here meet this requirement.) The parameter $\epsilon$ acts as a thermodynamic field with a conjugate density $\Theta$. Surface thermodynamics follows from

[TABLE]

where $s$ is the surface excess entropy per unit area and the conjugate density is Evans and Parry (1989)

[TABLE]

Using the Maxwell relation resulting from (10) we obtain the sum rule:

[TABLE]

Eqs. (8,12) are satisfied identically within the DFT approximations that we employ. We determine $\Gamma$ and $\Theta$ in our GCMC simulations and shall use (12) to explore the statistical accuracy of these. Note that the equivalent of (12) is well-known in surface (Ising) magnetism where $\Gamma$ is equivalent to the excess magnetization $m_{s}$, $\epsilon$ plays the role of the locally applied surface magnetic field and $\mu$ plays the role of the bulk magnetic field $h$. For the case of a surface field $h_{1}$ acting in only the first (surface) layer of spins one has the standard result:

[TABLE]

where $m_{1}$ is the magnetization in the surface layer. Clearly this particular magnetic case corresponds to a model fluid in which $W\epsilon)/\epsilon\sim\delta(z)$ so that $\Theta$ reduces to $m_{1}$. In the magnetism literature the quantity $(\partial m_{1}/\partial h)_{T}$ is usually termed $\chi_{1}$, the surface layer susceptibility.

For fluids it is appropriate to define a local compressibility:

[TABLE]

This quantity was introduced in early studies of wetting transitions e.g. Tarazona and Evans (1982); Evans and Parry (1990) and shown recently Evans and Stewart (2015); Evans and Wilding (2015); Evans et al. (2016) to provide a valuable measure of the degree of solvophobicity or hydrophobicity of a substrate, i.e. as the well-depth $\epsilon$ decreases and the macroscopic contact angle increases, there is an accompanying increase in the maximum of $\chi(z)$ located at distances $z$ within $1$-$3$ diameters of the substrate. The connection between $\chi(z)$ and density correlations at the substrate is best made by introducing the local, or transverse, structure factor:

[TABLE]

where $q$ is the transverse wave number. $S(z;q)$ provides a measure of the strength and range of transverse correlations at distance $z$ from the substrate and plays an important role in the theory of wetting Tarazona and Evans (1982); Evans and Parry (1990) and in characterizing the structure of the liquid-vapor interface Parry et al. (2016). It is straightforward to show that $\chi(z)$ is proportional to the $q=0$ limit of (15):

[TABLE]

i.e. the local compressibility $\chi(z_{1})$ is the integral of the density-density correlation function over the transverse coordinate $R$ and over one normal coordinate, $z_{2}$. Generally, $\chi({\bf r})$ can be expressed as a fluctuation formula that follows by differentiating the grand partition function w.r.t. an external potential, to obtain the average one-body density, and then differentiating w.r.t. the chemical potential. One finds Evans and Wilding (2015):

[TABLE]

Clearly $\chi({\bf r})$ is the correlator of the local number density at ${\bf r}$ and the total number $N$ of particles.

Finally if we integrate (16) over $z_{1}$ we obtain the surface compressibility sum rule,e.g. Nicholson and Parsonage (1982); Evans and Parry (1989), relating the surface excess compressibility in (9) to an integral of the surface structure factor, at $q=0$:

[TABLE]

where the bulk contribution $\chi_{b}=\rho_{b}^{2}\kappa_{T}$; $\kappa_{T}$ is the usual isothermal compressibility, proportional to the bulk structure factor at zero wave number.

III.2 Phenomenology of Wetting, Drying and Surface Criticality

Here we remind readers of some of the phenomenology of wetting and drying transitions and describe the critical exponents that characterize such transitions. Wetting and drying are phase transitions that occur strictly in the limit of infinite wall area, $A\to\infty$, and for a single wall, i.e infinite wall separation, $D\to\infty$. Using the adsorption language of the previous subsection, a wetting transition occurs at fixed $T$ when the excess adsorption $\Gamma$ changes from a finite value to an infinite value as the wall-fluid attraction $\epsilon_{w}$ is increased to the transition value $\epsilon_{ww}$; the bulk fluid is a vapor, at fixed chemical potential $\mu=\mu_{co}^{-}(T)$. For $\epsilon_{w}<\epsilon_{ww}$ the (planar) surface tensions satisfy $\gamma_{wv}<\gamma_{wl}+\gamma_{lv}$, where $w$ refers to wall, $l$ to liquid and $v$ to vapor. From Young’s equation (1) it then follows that $\cos(\theta)<1$ and the situation corresponds to partial wetting in Fig. 1.

For $\epsilon_{w}>\epsilon_{ww},\Gamma=\infty,\gamma_{wv}=\gamma_{wl}+\gamma_{lv}$ and $\cos(\theta)=1$, consistent with the wall-vapor interface being wet by a macroscopically thick layer of liquid. If $\Gamma$ jumps abruptly at the transition value the transition is first order. If $\Gamma$ diverges continuously at $\epsilon_{ww}$ the transition is termed continuous or critical. As mentioned in the Introduction, drying is the counterpart of wetting when the bulk fluid is a liquid at $\mu=\mu_{co}^{+}(T)$. For small wall-fluid attraction the local density near the wall can be depleted so that the adsorption, $\Gamma$, as defined by (8), is negative but finite and one finds $\gamma_{wl}<\gamma_{wv}+\gamma_{lv}$. This corresponds to the partial drying situation in Fig. 1, where $\cos(\theta)>-1$. On reducing $\epsilon_{w}$ further, to a value $\epsilon_{wd}$, a drying transition can occur whereby for $\epsilon\leq\epsilon_{wd},\Gamma=-\infty$, and $\gamma_{wl}=\gamma_{wv}+\gamma_{lv}$ , i.e. $\cos(\theta)=-1$, consistent with the wall-liquid interface being wet by a macroscopically thick layer of vapor.

We focus on drying and suppose there is a critical drying transition. The divergence of the adsorption is described by

[TABLE]

with $\delta\epsilon_{w}\equiv\epsilon_{w}-\epsilon_{wd}$, which is accompanied by a divergence of the parallel (transverse) correlation length $\xi_{\parallel}$:

[TABLE]

i.e. density fluctuations parallel to the wall become long-ranged on approaching the transition. This is most easily understood in terms of the surface structure factor introduced in (15). One expects Ornstein-Zernike behaviour: $S(z;q)=S(z;0)/(1+\xi_{\parallel}^{2}q^{2})$, for small wave numbers $q$, when $z$ is located close to $l$, the thickness of the vapor film. The latter is given by $l=-\Gamma/(\rho_{l}-\rho_{g})$. Moreover, considerations of capillary wave fluctuations in the emerging liquid-vapor interface lead to the prediction Evans and Parry (1989):

[TABLE]

Note that $\rho^{\prime}(l)$ is the gradient of the density profile of the emerging gas-liquid interface, as $l\to\infty$. Capillary wave arguments then predict $\rho^{\prime}(l)\sim\xi_{\perp}^{-1}$ as $\delta\epsilon_{w}\to 0$ where $\xi_{\perp}$ is the width of the depinning gas-liquid interface, i.e. the interfacial roughness. In spatial dimension $d=3$ one has $\xi_{\perp}^{2}\sim(2\pi\beta\gamma_{lv})^{-1}\ln(\xi_{\parallel}/\xi_{b}$), where $\xi_{b}$ is the bulk correlation length of the phase that wets; in the present case the vapor. Using the relation (18), i.e. integrating directly (21), we expect the surface excess compressibility to diverge as

[TABLE]

It is also important to consider the quantity $\chi_{1}\equiv(\partial\Theta/\partial\mu)_{T}$ introduced in (12). Clearly this corresponds to an integral of $\chi(z)$ weighted with the wall-fluid potential; see (11). The Maxwell relation (12) with (19) dictates that

[TABLE]

The other key exponent, $\alpha_{s}$, is associated with the singular part of the surface tension $\gamma=\Omega_{ex}/A$. For drying one has $\gamma_{wl}\equiv\gamma_{wv}+\gamma_{lv}+\gamma^{sing}$ with

[TABLE]

Note that from Young’s equation (1) it follows that $\gamma^{sing}=-\gamma_{lv}(1+\cos(\theta)$; this quantity is negative in the partial drying regime. The three critical exponents are not independent. A scaling hypothesis for $\Omega_{ex}$ and a thermodynamic argument employing (10) Evans and Parry (1989) both yield the analogue of the well-known Rushbrooke exponent (in)equality for the corresponding bulk exponents, i.e.

[TABLE]

The critical exponents depend on the dimensionality $d$ and on the ranges of the fluid-fluid ($ff$) and wall-fluid ($wf$) potential, e.g Dietrich (1988). We assume that the hyperscaling relation $2-\alpha_{s}=(d-1)\nu_{\parallel}$ is valid for the $d-1$ dimensional interface of the $d$ dimensional fluid for $d\leq d_{c}$, the upper critical dimension. If $ff$ and $wf$ potentials are both of finite range, as in the case of a truncated LJ fluid (2) near a square-well wall (4), mean-field (MF) analysis, e.g. Dietrich (1988) yields $\alpha_{s}=0,\beta_{s}=0$ (logarithmic divergence) and $\nu_{\parallel}=1$. These MF results for SR potentials are consistent with (25) and when inserted into the hyperscaling relation yield $d_{c}=3$. As outlined in the Introduction, RG calculations, Ising model simulations and concerted theoretical effort to incorporate relevant fluctuation effects conclude that the critical exponents $\nu_{\parallel}$ and $\alpha_{s}$ should depend on the dimensionless parameter $\omega=(4\pi\beta\gamma_{lv}\xi_{b}^{2})^{-1}$ that measures the strength of interfacial fluctuations. Mean field corresponds to an infinitely stiff interface: $\omega=0$.

By contrast, when both $ff$ and $wf$ potentials exhibit appropriate algebraic (power-law) decay critical drying can occur and one finds, using hyperscaling, that the upper critical dimension $d_{c}<3$. Thus MF results for critical exponents should remain valid in $d=3$ when both potentials are LR. An explicit DFT calculation for such a situation is described in Stewart and Evans (2005). In the present paper we focus primarily on the case of a truncated $ff$ potential and a $wf$ potential with attraction that decays algebraically, as $\sim z^{-3}$, see (5). For this special case, a MF binding potential analysis (see Sec.IV.1) yields MF exponents which i) are different from those pertinent to SR potentials described above, ii) satisfy (25) and iii) when inserted into the hyperscaling relation, imply $d_{c}=3$. The same calculation finds that (critical) drying occurs in the limit $\epsilon_{w}\to 0$, i.e. where the $wf$ attraction is vanishing and the (microscopic) DFT calculations presented in Sec. V confirm this result. A simple RG treatment of fluctuation effects, Sec. IV.2, finds that the critical exponents are unchanged from their MF values.

III.3 Phase diagram of a fluid in a slit pore with identical walls

In order to set the scene and assist the reader, we outline here general features of the surface phase behaviour of a simple fluid in a slit pore having wall separation $D$, which is the system that we consider in the present simulations. The natural choice of variables for representing the phase behaviour is the pair of fields: $\delta\mu\equiv\mu-\mu_{co}$, which measures the deviation of the chemical potential from its bulk coexistence value; and $\epsilon_{w}$, which measures the strength of the wall-fluid attraction. Fig. 2 shows a schematic sketch of one possible phase diagram for such a system at some temperature $T<T_{c}$, with $T_{c}$ the bulk critical temperature.

The system exhibits two lines of phase coexistence: the capillary line and the prewetting line. The capillary line is the locus of state points for which gas-liquid coexistence occurs in the slit pore. It is given approximately by the Kelvin equation

[TABLE]

where $\rho_{l}$ and $\rho_{v}$ are the coexisting liquid and vapor densities, respectively. Small values of $\epsilon_{w}$ favor the capillary gas phase while large values favor the capillary liquid. If the transition occurs for $\delta\mu<0$ it is referred to as capillary condensation ($cc$); if the transition occurs for $\delta\mu>0$ it referred to as capillary evaporation ($ce$). The capillary condensation/evaporation line is sketched within a range of $\epsilon_{w}$ between the drying wall strength $\epsilon_{wd}$ at the lower end and the wetting wall strength $\epsilon_{ww}$ at the upper end. The line in Fig. 2 extends smoothly to $\epsilon_{w}<\epsilon_{wd}$ and to $\epsilon_{w}>\epsilon_{ww}$, as shown. The detailed behaviour depends on higher order contributions in (26). At $\epsilon_{w}\geq\epsilon_{ww}$ the contact angle is zero, i.e. $\cos(\theta)=1$, while at $\epsilon\leq\epsilon_{wd}$ it is $180^{\circ}$ i.e. $\cos(\theta)=-1$. The neutral wall, at which $\cos(\theta)=0$, corresponds to the value of $\epsilon_{w}$ at which the capillary line crosses the $\delta\mu=0$ axis. In the sketch in Fig 2 we have drawn the situation that occurs for critical drying and first order wetting, see Sec IV.1. This scenario pertains to the LR wall-fluid potential that we investigate here, though for the particular model that we consider $\epsilon_{wd}=0$. Note however, that other scenarios are possible. Specifically, and as we show below, wetting can also be continuous (critical). The prewetting ($pw$) line shown occurs only for first order wetting. It emerges tangentially from the wetting point $\delta\mu=0^{-},\epsilon_{w}=\epsilon_{ww}$ and extends some distance to $\delta\mu<0$ before terminating at a prewetting critical point whose critical properties correspond to the universality class of the 2d Ising model Nicolaides and Evans (1989). Along the prewetting line a thin layer of liquid on each wall coexists with a thick but finite liquid layer Dietrich (1988); Fan and Monson (1993). However we will not consider prewetting in the present work.

Eq. (26) proves a useful estimate of the range of $\beta\delta\mu$ for which capillary condensation and evaporation can be found in our system. At drying or wetting (26) yields $-\delta\mu_{cc}^{w}=\delta\mu_{ce}^{d}=2\gamma_{lv}/D(\rho_{l}-\rho_{v})$ and the horizontal scale in Fig 2 is determined by $\beta\mu_{ce}^{d}\approx 2\beta\gamma_{lv}\sigma^{2}/[D\sigma^{-1}(\rho_{l}\sigma^{3}-\rho_{v}\sigma^{3})]$. For our simulated LJ fluid at $T/T_{c}=0.775$ the liquid-vapor surface tension is given by $\beta\gamma_{lv}\sigma^{2}=0.404$ and $\beta\mu_{ce}^{d}=1.295/D\sigma^{-1}$. Thus for $D=30\sigma$, the wall separation in most of the simulations, we find $\beta\mu_{ce}^{d}\approx 0.04$. The physics we describe is occurring at small under or over saturations. But these values of $\beta\delta\mu$ are certainly pertinent to experiment.

IV Theory for a Model Fluid with SR Fluid-Fluid and LR Wall-Fluid Potentials

IV.1 Binding potential analysis

We follow the standard treatment, e.g. Dietrich (1988), of wetting/drying transitions and consider $\omega^{ex}(l)$, the excess grand potential per unit surface area, as a function of the thickness $l$ of the wetting/drying layer. For a truncated LJ model adsorbed at a single wall exerting the potential (5) or the modified version, we expect

[TABLE]

with the binding potential

[TABLE]

As previously, $\rho_{l}$ and $\rho_{v}$ are the liquid and vapor densities at coexistence, $\delta\mu=\mu-\mu_{co}\geq 0$ is the deviation of the chemical potential from its value at coexistence and we have specialized now to the case of drying, i.e. $l$ is the thickness of a layer of vapor that can intrude between the weakly attractive wall and the bulk liquid at $z=\infty$. In the limit of complete drying, at $\delta\mu=0^{+}$, $l$ diverges and the $wl$ interface is a composite of the $wv$ and $lv$ interfaces. In this limit $\gamma_{wl}=\gamma_{wv}+\gamma_{lv}$, i.e. $\cos(\theta)=-1$, as mentioned previously in Sec. III.2. The binding potential in (28) has two leading contributions. The exponential term accounts for SR fluid-fluid interactions; $\xi_{b}$ is the true correlation length of the bulk phase that wets, in our case the vapor, and $a$ is a positive coefficient. The term $bl^{-2}$ is associated with the $z^{-3}$ decay of $W_{LR}(z)$ in (5); it arises from dispersion (van der Waals) forces between the substrate and the fluid. The higher order terms in (28) include higher inverse powers such as $cl^{-3}$ as well as more rapidly decaying exponentials. We ignore these in the subsequent analysis. Making a straightforward sharp-kink approximation, or Hamaker type calculation, e.g.  Dietrich (1988), yields

[TABLE]

Since $b<0$ for all $T<T_{c}$, minimizing (27) w.r.t. $l$ at $\delta\mu=0^{+}$, leads to a finite value for the equilibrium thickness:

[TABLE]

A formula equivalent to (30) was derived by Nightingale et al. (see Eq. 6 of Nightingale et al. (1983)) in a study of critical wetting in systems with LR forces. Those authors considered only the case where $\epsilon_{w}>0$ and concluded there was no wetting, critical or first order. Here we focus on the situation where $\epsilon_{w}\to 0^{+}$, $b\to 0^{-}$ and $l_{eq}$ diverges continuously. Note that for $\epsilon_{w}=0$, $W_{LR}(z)$ in (5) reduces to the planar hard-wall potential and minimization of (27) then yields $-l_{eq}/\xi_{b}=\ln(\delta\mu)+{\rm const}$, the mean-field (MF) result appropriate for complete drying from off-coexistence, for all $T<T_{c}$, e.g. Evans (1992); Dietrich (1988); Evans and Parry (1990).

Using (27,28) we can calculate several properties and examine these, within MF, in the approach to critical drying $\epsilon_{w}\to 0^{+}$. The local compressibility, evaluated for $z\approx l_{eq}$, is given by Evans and Stewart (2015); Evans and Parry (1990)

[TABLE]

where the prime denotes differentiation w.r.t. $z$. From (27) it follows that, at leading order,

[TABLE]

Capillary wave arguments predict that in the limit of critical drying $\rho^{\prime}(l_{eq})\sim\xi_{\perp}^{-1}$, where $\xi_{\perp}$ is the interfacial roughness introduced above. Within MF, $\xi_{\perp}^{-1}$ is non-zero, and using (30) we deduce

[TABLE]

The predictions (30) and (33) were tested carefully using the microscopic DFT, as described below. The quantity $\Gamma_{1}$ defined in (12), is proportional to $-\partial l_{eq}/\partial\epsilon_{w}$. It follows that as $\epsilon_{w}\to 0^{+}$,

[TABLE]

We can also extract the correlation length $\xi_{\parallel}$ that describes density-density correlations parallel to the wall. General arguments predict that $\xi_{\parallel}^{2}$ diverges in the same way as the surface excess compressibility, see (22).

Since $\chi_{ex}$ is proportional to $-\partial l_{eq}/\partial\mu$ it follows from (32,30) that $\xi_{\parallel}$ diverges as

[TABLE]

in the limit $\epsilon_{w}\to 0^{+}$. The same result is obtained from standard binding potential considerations Dietrich (1988) where one has $\xi_{\parallel}^{-2}\propto\partial^{2}\omega_{B}(l)/\partial l^{2}$ at $l=l_{eq}$.

The variation of $\cos(\theta)$ close to critical drying is determined by $\omega_{B}(l_{eq})$ at $\delta\mu=0^{+}$, i.e. the singular part of the surface excess free energy $\gamma^{sing}$. Using Young’s equation (1) one finds $1+\cos(\theta)=-\omega_{B}(l_{eq})/\gamma^{lv}$ and for the present binding potential (28) we obtain

[TABLE]

in the limit $\epsilon_{w}\to 0^{+}$. This result is striking. Were the logarithm not present in (36) the theory would predict $1+\cos(\theta)$ vanishing linearly with $\epsilon_{w}$, a signature of a 1st order drying transition. It is only the presence of the logarithm that ensures a continuous (critical) transition. The critical exponent $\alpha_{s}$, defined by the vanishing of the singular part of the surface excess free energy $|\gamma^{sing}|\sim\epsilon_{w}^{2-\alpha_{s}}$, see (24), clearly takes the value $\alpha_{s}=1$, with log corrections, in this particular case. The situation is similar to that for complete drying from off-coexistence where for a planar hard-wall, say, $\gamma^{sing}\sim\delta\mu\ln\delta\mu$, $\delta\mu\to 0^{+}$.

It is important to distinguish the MF scenario presented above from that corresponding to a SR wall-fluid potential such as (4). In the SR case it is well-known, e.g. Brézin et al. (1983); Dietrich (1988), that the second inverse power-law term in (28) must be replaced by a H.O. term proportional to $\exp(-2l/\xi_{b})$ while the coefficient of the leading $\exp{(-l/\xi_{b})}$ term now depends on $\epsilon_{w}$: $a(\epsilon_{w})\sim(\epsilon_{w}-\epsilon^{MF}_{wc})$ where $\epsilon^{MF}_{wc}>0$ is the strength of the wall-fluid attraction at which critical drying occurs in MF. Defining $\delta\epsilon_{w}=\epsilon_{w}-\epsilon^{MF}_{wc}$, MF analysis for the SR case yields, for $\delta\mu=0^{+}$,

[TABLE]

and

[TABLE]

The critical exponents $\beta_{s},\nu_{\parallel}$ and $\alpha_{s}$ take the values mentioned in Sec. III.2. These results are clearly very different from those we obtained above for the LR case.

It is also important to consider the implications of a binding potential of the form (28) for wetting. In this context we recall an argument of Ebner and Saam (ES) Ebner and Saam (1987b) who considered a lattice gas (Ising) model for which the substrate-fluid potential is $W_{n}=-RJn^{-p}$, $p\geq 3$, where $n$ labels the $n$th layer from the substrate and $RJ>0$ is a constant. For the case of a SR $ff$ potential the most slowly decaying term in the ES grand potential, i.e. the binding potential, has the form $(\rho_{\alpha}-\rho_{\beta})RJ/(\rho-1)l^{p-1}$ where $\rho_{\alpha}$ is the density and $l$ is the thickness of the phase $\alpha$ that wets the substrate. In the case of wetting by the denser ‘liquid’ $\rho_{\alpha}>\rho_{\beta}$, the density of the ‘vapor’, and this term $>0$, implying the binding potential can have a relative minimum at $l=\infty$. Suppose that there is incomplete wetting, for some value of $RJ$. Then ES argue the binding potential must have a minimum, lower than at $l=\infty$, for a ‘liquid’ film of finite thickness. On increasing $RJ$, equivalent to increasing $\epsilon_{w}$ in our system, a wetting transition can occur but this can only be first order: there cannot be a continuous evolution from finite $l$ to infinite $l$. Thus if a wetting transition occurs this cannot be critical. For sufficiently large $RJ$, or $\epsilon_{w}$, one expect on physical grounds that wetting should occur. It follows this must be first order.

ES also consider drying where $\alpha$ now corresponds to ‘vapor’ and $\beta$ to ‘liquid’. Now the relevant term in the binding potential is negative, as given by (29), and there is a relative maximum at $l=\infty$. ES then argue that drying cannot occur for any $T<T_{c}$. They do not consider the limit $RJ\to 0$, corresponding to our present limit $\epsilon_{w}\to 0^{+}$.

The numerical work of ES for the lattice gas model confirms that wetting is always first order and ES find no critical or first order drying transitions. In our present DFT and simulation studies we find that for $W_{LR}(z)$ in (5) the wetting transition is first order. In contrast to ES, we do find a critical drying transition. This occurs as the attractive strength $\epsilon_{w}\to 0^{+}$. In this limit our wall-fluid potential reduces to that of a hard-wall, for which drying occurs for all $T<T_{c}$. This hard-wall boundary condition, particular to fluids, drives the drying transition. Both DFT and simulation find critical drying as $\epsilon_{w}\to 0^{+}$.

IV.2 Renormalization Group (RG) treatment of fluctuations

The analysis described in Sec. IV.1 was strictly MF; this omits some of the effects of capillary wave (CW) fluctuations. For example, for infinite surface area, MF predicts a sharp interface with $\xi_{\perp}$ finite in all dimensions $d$ whereas, in reality, we expect $\xi_{\perp}$ to diverge for $d\leq 3$. An important early attempt to incorporate CW fluctuations was that of Brezin et al. Brézin et al. (1983) who introduced a RG treatment for the case of SR forces where the upper critical dimension $d_{c}=3$ for both critical wetting and complete wetting from off-coexistence. We follow their methodology for our binding potential (28).

First we invoke the hyperscaling relation $(2-\alpha_{s})=(d-1)\nu_{\parallel}$, where $\nu_{\parallel}$ is the critical exponent for $\xi_{\parallel}$, insert the MF exponents given in Sec. IV.1, and deduce that the upper critical dimension is $d_{c}=3$ for the present system. Introducing again the standard, dimensionless parameter $\omega$, that measures the strength of CW fluctuations, the RG treatment then implies we should consider an effective binding potential (renormalized) at the scale $\xi_{\parallel}$:

[TABLE]

The exponential term is renormalized but the remaining power-law terms are not; in particular the coefficient $b$ is assumed to be unchanged. Minimization of (42) yields

[TABLE]

as $\epsilon_{w}\to 0^{+}$. The equilibrium thickness still diverges with the MF form (30) but the amplitude is increased by a factor $(1+\omega/2)$. MF is recovered when the interface becomes very stiff so that $\omega\to 0$. The parallel correlation length can be obtained from either $\xi_{\parallel}^{-2}\propto\left(\frac{\partial^{2}\omega_{B}(l)}{\partial l^{2}}\right)$ at $l=l_{eq}$ or from $\xi_{\parallel}^{2}\propto\left(\frac{\partial l_{eq}}{\partial\mu}\right)$, see (22). In both cases we find as $\epsilon_{w}\to 0^{+}$

[TABLE]

The singular part of the surface excess free energy can be calculated from (42) and we obtain

[TABLE]

Once again only the amplitudes are changed from the MF results (35) and (36). Note that (43) is reminiscent of the result for complete drying from off-coexistence for SR forces, e.g. at a planar hard-wall. There the second term in the r.h.s. of (42) is absent but the third remains leading to

[TABLE]

Unlike the case of SR forces considered by Brezin et al. Brézin et al. (1983) and in many subsequent studies, e.g. Binder et al. (1989); Albano and Binder (2012); Parry et al. (2008b, a); Fisher and Huse (1985); Parry and Rascón (2009) where several of the critical exponents for critical wetting are predicted to depend explicitly on the parameter $\omega$, for the binding potential (28) our RG analysis predicts the critical exponents to be unchanged from their MF values and therefore independent of $\omega$ even though the upper critical dimension is also $d_{c}=3$. We note that the conclusions of the MF and RG analyzes are changed little if we consider LR wall-fluid potentials other than the standard $9$-$3$ case (5). Suppose the leading wall-fluid power-law decay is proportional to $-(\sigma/z)^{p}$, with $p>2$. Then the coefficient of the second term in (30) is replaced by $p$, (33) is unchanged, and the power of the logarithm in (35) and (36) is replaced by $p/2$ and $-(p-1)$, respectively. The RG results are changed accordingly.

IV.3 DFT Treatment

The classical DFT that we employ is that used in a previous study of solvophobic substrates but one that did not address critical drying Evans and Stewart (2015). The excess Helmholtz free energy functional is approximated by the sum of a hard-sphere functional, treated by means of Rosenfeld’s fundamental measure theory, and a standard MF treatment of attractive fluid-fluid interactions. Eq. (14) of Ref. Evans and Stewart (2015) displays the grand potential functional. This form of the functional has been used in many studies of fluid interfacial phenomena, including wetting and capillary confinement; Ref Evans and Stewart (2015) provides pertinent references. In the present study the attractive part of the truncated LJ potential is given by

[TABLE]

where $r_{\rm min}=2^{1/6}\sigma$. The potential is truncated at $r_{c}=2.5\sigma$, as in simulation. The critical temperature is given by $k_{B}T_{c}=1.3194\epsilon_{LJ}$ and calculations are performed at $T=0.775T_{c}$. The LR wall-fluid potential is the standard $9$-$3$ model given by $W_{LR}(z)$ in (5). We also investigate the SR case (4). The hard-sphere diameter, entering the hard sphere functional, is $d=\sigma$.

In the DFT calculations we determine equilibrium density profiles $\rho(z)$ and the surface tensions $\gamma_{lv},\gamma_{wl},\gamma_{wv}$, by minimizing the grand potential functional Evans and Stewart (2015). The local compressibility $\chi(z)$ defined in (14) is determined numerically as described in Ref. Evans and Stewart (2015). We have performed calculations for a single wall and for a pair of confining walls, equivalent to the GCMC simulations. In Sec. V.1 we show results for the single wall and in Sec. V.2 for two walls.

V Results from DFT

V.1 Fluid adsorbed at a single planar wall

Our key results are shown in Fig. 3. Here we plot $\rho(z)$ and $\chi(z)$ for very small values of $\epsilon_{w}$ at the LR wall (5). As $\epsilon_{w}$ is reduced towards zero the thickness of the drying film $l_{eq}$ increases (top). We have confirmed in detail, within DFT, that the Gibbs adsorption $\Gamma$ or $l_{eq}$ grows according to (30). This is illustrated in Fig. 4 where we plot $l_{eq}$ vs $\ln\epsilon_{w}-3\ln(\l_{eq})$. The slope yields $\xi_{b}=0.51\sigma$ for the correlation length of the bulk (vapor) phase that wets. This estimate is close to that from a separate DFT calculation of the binding potential. The position of the peak in $\chi(z)$ shifts with the position of the gas-liquid interface and its height increases very rapidly as $\epsilon_{w}\to 0^{+}$ (middle). The bottom panel shows clearly that $\ln\chi(l_{eq})$ increases linearly with $l_{eq}$. The prediction (33), including the correct prefactor, the inverse bulk correlation length, is confirmed by our DFT calculations. We determine the contact angle via Young’s equation and DFT results for $\cos(\theta)$ are shown in Fig. 5. Note that in order to facilitate comparison with subsequent plots of simulation results, the quantity plotted on the abscissa of Fig. 5 is the wall-fluid potential strength in units of $k_{B}T$, i.e. for the LR case (5) this denotes $\epsilon_{w}\epsilon_{LJ}/k_{B}T=0.9779\epsilon_{w}$ and for the square well (4) denotes the well-depth in units of $k_{B}T$. For the LR case (5) we find critical drying at $\epsilon_{w}=0$ and 1st order wetting at a value of $\epsilon_{w}$ that is smaller than in simulation, see later. For the SR case (square-well), both drying and wetting are critical transitions, as found in simulation. However, as we shall see, the separation in $\epsilon_{w}$ between wetting and drying in DFT is smaller than in simulation. The microscopic DFT results for a single LR wall yield the same behavior as those from the simple binding potential treatment, based on (28), i.e., for the LR case, the DFT yields the same MF critical exponents, including any logarithmic ($\ln\epsilon_{w}$) corrections, as those predicted in Sec. IV.1. This is not unexpected: DFT is a MF treatment of fluid interfaces and we expect it to capture the same asymptotic, $|\Gamma|$ or $l_{eq}\to\infty$, behaviour as the binding potential analysis. The key difference between the two approaches lies in the fact that DFT incorporates accurately the short distance behaviour of the density profile. In particular our DFT satisfies exactly the hard-wall sum rule: $k_{B}T\rho(0^{+})=p(\mu)$, where $p$ is the pressure of the bulk fluid, that is important in ensuring that complete drying occurs in the limit $\epsilon_{w}\to 0^{+}$. This is mimicked in the binding potential by the first, repulsive, term in (28).

In our DFT calculations we can compute accurately the excess grand potential $\omega^{ex}(\Gamma)$, for non-equilibrium values of the adsorption $\Gamma$, during the minimization of the functional. This is equivalent to determining numerically the binding potential entering (27). For the LR case $\omega^{ex}(\Gamma)$ exhibits two minima, corresponding to a microscopic liquid film and an infinitely thick liquid film, on approaching the wetting transition. At the transition the two minima are equal but there remains a maximum between these – a clear signature that the transition is first order and in keeping with the lattice gas results of ES Ebner and Saam (1987b). On the other hand, on approaching the drying transition $\omega^{ex}(\Gamma)$ exhibits a single minimum at $|\Gamma|$ corresponding to a thick drying film. The minimum erodes continuously and shifts to larger $|\Gamma|$ as $\epsilon_{w}$ is reduced. In the limit $\epsilon_{w}\to 0^{+}$ the minimum is at $|\Gamma|=\infty$ , i.e. the transition is critical as predicted by the binding potential treatment. For the SR case $\omega^{ex}(\Gamma)$ exhibits a single minimum in the approach to both wetting and drying showing that both are critical transitions. We have not attempted to determine critical exponents numerically for the SR case because of the difficulty of locating accurately the drying and wetting points. However, there is no reason to expect the exponents from DFT to differ from those given by the MF analysis of the binding potential, i.e. (37)-(41).

V.2 Fluid adsorbed between two planar walls

Since the GCMC simulations, to be described in Secs. VI and VII, investigate fluids confined between two planar walls, we also performed some DFT calculations for the confined system. We focused on the LR case, with each wall described by the potential (5). The wall-separation was chosen to be $D=30\sigma$ and the temperature was $T=0.775T_{c}$ which correspond to the system studied in most detail in the simulations. Results for the density profiles and the local compressibility are shown in Fig. 6 for a range of values of $\epsilon_{w}$. Note that the smallest value is well-removed from the drying point $\epsilon_{w}=0^{+}$; for $\epsilon_{w}=0.3$ we find the contact angle $\theta\approx 157^{\circ}$; see Fig. 5. We observe the erosion of oscillations in the density profile and the growth of a depleted region of density at each wall as $\epsilon_{w}$ is reduced in this range. This is accompanied by the smoothing of oscillations in the local compressibility. As $\epsilon_{w}$ is reduced, the height of the maximum in the local compressibility increases strongly and its location shifts to larger distances from the wall, following the position of the maximum gradient of the density profile. These trends are consistent with DFT results in Ref. Evans and Stewart (2015) for a different (SR) wall-fluid potential but pertaining to a similar range of contact angles. It is important to note that the results in Fig. 3 for a single LR wall correspond to tiny values of the wall-fluid attraction where the drying film is very thick; there we test the detailed predictions of the binding potential description in the limit of drying. In Fig. 6 we are examining the overall changes of the density profiles and local compressibility as the substrate becomes more solvophobic and the contact angle becomes very large. We should also note that the results in Fig. 6 are for the liquid at bulk coexistence, i.e. $\delta\mu=0^{+}$. For the confined fluid this state is metastable w.r.t. capillary evaporation. The latter would occur at values of $\beta\delta\mu$ that are typically about $0.04$ –see Sec. III.3. Within DFT there is no difficulty in probing these metastable ‘liquid’ states which correspond to local minima of the excess grand potential. This is illustrated, for smaller wall separations, in Fig. 8 of Ref. Evans and Stewart (2015).

VI Simulation methods

We employ GCMC simulation, which is well suited to studying fluids at vapor-liquid coexistence both in the bulk Wilding (1995) and in confinement Cracknell et al. (1993); Maciolek et al. (2003); Grzelak and Errington (2008); Rane et al. (2011). Within this framework one prescribes the temperature $T$ and chemical potential $\mu$, while the particle number $N$ fluctuates. The relevant observables are the probability function $P(\rho)$ of the total density $\rho=N/V$ and, for a confined system, the density profile $\rho(z)$.

Although $\rho$ fluctuates in our simulations, close to coexistence state points sampling problems can arise due to the free energy cost of traversing the mixed phase (interfacial) states that separate pure vapor from pure liquid. This cost is manifest as a deep valley of low probability in $P(\rho)$ which, on simulation timescales, traps the sampling in one phase. To overcome this problem we have implemented biasing techniques Berg and Neuhaus (1992), utilizing a weight function that is calculated from the transition matrix Smith and Bruce (1995). The role of the weight function is to enhance the sampling of mixed states, ie. to remove the sampling barrier. The effects of the biasing can subsequently be unfolded exactly from distributions of observables.

Sufficiently close to the bulk vapor-liquid critical point, the requisite weights can be taken to be a function of the total density $\rho$. However, at low temperature coexistence points (such as used in the present work) this approach breaks down due to the appearance of finite-size induced first order phase transitions known as droplet transitions Binder et al. (2012). For the purpose of biasing through these transitions, $\rho$ is not a good order parameter and more effective alternatives must be sought. Suitable substitutes have recently been proposed by one of us Wilding (2016), and these were adopted in the present work.

GCMC is most efficient when deployed in conjunction with histogram extrapolation Ferrenberg and Swendsen (1989). This permits the results from a simulation performed at one set of model parameters, e.g. $T,\mu$ and wall strength $\epsilon_{w}$, to be reweighted to provide estimates of observables at nearby parameters, without recourse to further simulation. In the present work we have used histogram extrapolation to measure the local compressibility $\chi(z)$ from the $\mu$ dependence of the density profile $\rho(z)$. We have also used it in conjunction with results for $P(\rho)$ for a fully periodic system to obtain accurate estimates of the coexistence chemical potential at a given temperature: tuning $\mu$ until the equal peak weight criterion is satisfied Wilding (1995).

VII Results from simulations

VII.1 $P(\rho)$ and the contact angle

The contact angle as a function of wall strength $\epsilon_{w}$ can be obtained from the measured form of the density probability function $P(\rho)$ in both a fully periodic system and in the slit. Since $P(\rho)$ can vary over many decades it is convenient to work with its logarithm. Doing so has the additional advantage that the latter links directly to the grand potential which is given by $\beta\Omega(\rho)=-\ln P(\rho)$. Note that in all the simulation results that we present, $\epsilon_{w}$ refers to the wall-fluid potential strength measured in units of $k_{B}T$, i.e. for (5) this quantity denotes $\epsilon_{w}\epsilon_{LJ}/k_{B}T=1.0877\epsilon_{w}$ and for (4) denotes the well-depth in units of $k_{B}T$. Figure 7(a) shows our GCMC results for $\ln P(\rho)$ for the slit system at vapor-liquid coexistence at a range of wall strengths $\epsilon_{w}$ which span the regime from wetting to drying. The data shown is for the modified LR potential, but a similar scenario plays out for the SR potential. For sufficiently large $\epsilon_{w}$, a double peaked structure is evident in $P(\rho)$. The low density peak (of height $P_{\rm vap}$) corresponds to the system in a capillary vapor phase, while the high density peak (of height $P_{\rm liq}$) corresponds to the capillary liquid phase. For large $\epsilon_{w}$ the liquid has the higher peak (i.e. is the stable phase), but as $\epsilon_{w}$ is reduced, the height of the liquid peak diminishes progressively, until it becomes metastable with respect to the vapor. The liquid peak height continues to diminish as $\epsilon_{w}$ is decreased until eventually it disappears into a plateau. At still smaller $\epsilon_{w}$, $P(\rho)$ is a monotonically decreasing function in the region of liquid-like densities, as shown in Fig. 7(b). Also included in Fig. 7 (dashed line) is the coexistence form of $P(\rho)$ for a fully periodic cubic system of side $L$ that matches the linear dimension of the planar walls in the slit system, c.f. eq. (3). For this latter system, $P(\rho)$ exhibits a pair of equal peaks of probability $P_{\rm max}$, corresponding to the respective pure phase states. These are separated by a central plateau (of probability $P_{\rm min}$) corresponding to mixed phase (interfacial) states.

As is well established, the ratio of peak to valley probabilities in the fully periodic case provides an accurate estimate of the vapor-liquid surface tension $\gamma_{lv}$ Binder (1982); Errington (2003):

[TABLE]

where $\beta=(k_{B}T)^{-1}$. Similarly the ratio of peak heights for the distributions in the slit system provides a measure of the surface tension difference Müller and MacDowell (2000):

[TABLE]

Accordingly, one can simply read off these quantities directly from the measured forms of $P(\rho)$ and insert them into Young’s equation (1) in order to obtain an estimate of the contact angle as

[TABLE]

This methodology can be used to estimate $\cos(\theta)$ as a function of wall strength $\epsilon_{w}$ for both the LR and SR wall-fluid potentials. The results are shown in Fig. 8 and should be compared with those of the DFT calculations of Fig. 5. A number of pertinent features are apparent. First we note that for both the LR and SR cases, $\cos(\theta)$ appears to approach $-1$ tangentially both in simulation and DFT. This behaviour is expected for a critical drying transition Dietrich (1988), see Secs. III.2 and IV.1. Second, in both the LR and SR cases, the simulations appear to indicate that drying occurs for a small but non-zero wall strength $\epsilon_{w}$, the value of which is substantially smaller for the LR case than the SR case. This appears to signal a qualitative discrepancy with DFT which unambiguously predicts (in accord with our binding potential calculations-see Sec. IV.1) that for the LR case critical drying occurs at $\epsilon_{w}=0$. We shall return to this discrepancy. Third, there are clear qualitative differences between drying and wetting for the LR case: While for the SR case $\cos(\theta)$ approaches the wetting limit $\cos(\theta)=1$ tangentially, similar to drying, and therefore indicating critical wetting, for the LR case the approach to this limit is with non-zero gradient, indicative of first order wetting Dietrich (1988). The DFT results shown in Fig. 5 are consistent with this finding.

Beyond their utility for determining contact angles, the relations (48)-(50) permit one to forge the link between principal features of Fig. 7 and the surface phase diagram (cf. Fig. 2). The wetting point ($\cos(\theta)=1$) occurs for $P_{\rm liq}/P_{\rm vap}=P_{\rm max}/P_{\rm min}$, which in the LR case occurs for $\epsilon_{w}\approx 2.6$. From fig. 7(a) one sees that while the liquid peak is strongly stable at this point, there is a metastable vapor peak which corresponds to a local free energy minimum. This minimum serves to bind the vapor phase to the wall at the wetting transition and hence –and in accord with the observed behaviour of the contact angle– the wetting transition is first order in this system 222A weak bulge is evident in $\ln P(\rho)$ for $\rho\sigma^{3}\simeq 0.2$ at the wetting point; we speculate that this may be a signature of the incipient prewetting transition..

The partial wetting regime, defined by $0<\cos(\theta)<1$, occurs for $1.8\lesssim\epsilon_{w}\lesssim 2.6$. Its lower boundary is marked by the ‘neutral’ wall for which $P_{\rm vap}=P_{\rm liq}$. This heralds entry into the partial drying regime ($-1\leq\cos(\theta)<0$) within which, for confinement within a slit, the liquid phase is metastable with respect to the vapor. Drying occurs for $\cos(\theta)=-1$ and corresponds to the wall strength for which $P_{\rm liq}/P_{\rm vap}=P_{\rm min}/P_{\rm max}$. Fig. 7(b) shows that the value of $\epsilon_{w}$ at which this equality is satisfied coincides closely with the point at which the liquid peak disappears smoothly into a plateau. We can therefore provisionally identify the drying point as that wall strength for which the liquid peak vanishes. The smoothness with which this occurs is consistent with the arguments presented above advocating that drying is critical.

VII.2 Finite-size scaling: pinning down the drying transition

Contact angle measurements provide an accurate indication of the order of a surface phase transition. However in the case of a critical surface phase transition, they do not yield accurate estimates for its location, nor for the associated critical exponents. The problem goes beyond the inherent difficulty of estimating the wall strength for which $\cos(\theta)=\pm 1$ when the approach to this limit is tangential; the main difficulty is one of finite-size effects. This is demonstrated in Fig. 9 which shows our GCMC estimates of $\cos(\theta)$ for values of $\epsilon_{w}$ in the vicinity of the critical drying transition, for three values of $L$; recall $L^{2}$ is the wall area. The data clearly show that the apparent drying point $\cos(\theta)=-1$ shifts systematically to lower values of $\epsilon_{w}$ as $L$ increases. This observation is important and we return to it later.

To clarify the nature of the near-critical finite-size effects, it is useful to return to the density distribution $P(\rho)$. In examining this quantity we shall exclude the low density region where the vapor peak occurs. This peak corresponds to the capillary evaporation transition that occurs when two vapor-liquid interfaces unbind from the wall and wander to the slit centre where they annihilate. By excluding it from the sampling we can focus on the behaviour at higher (liquid-like) densities, which are the ones relevant for critical drying. A further advantage is that we can switch from a logarithmic to a linear scale which is more revealing as regards exposing the character of the criticality.

Fig. 10 shows the measured forms of $P(\rho)$ for the system with LR wall-fluid interactions at a selection of system sizes $L$ and wall strengths $\epsilon_{w}$. One sees that for sufficiently large $\epsilon_{w}$ and $L$, $P(\rho)$ exhibits a liquid peak. On reducing $\epsilon_{w}$, this peak disappears into a plateau. On further reduction $P(\rho)$ becomes monotonically decreasing with a bulge which gradually diminishes until, at $\epsilon_{w}=0$, the distribution comprises a linear part and a tail. The interesting feature of Fig. 10 is that the presence or otherwise of a liquid peak at a given $\epsilon_{w}$ depends on the value of $L$. Even when $\epsilon_{w}$ is sufficiently small that only a bulge (but no peak) is seen, the bulge grows over the range of accessible $L$, suggesting that a peak would form were sufficiently large values of $L$ attainable.

The range of $\epsilon_{w}$ over which the distribution evolves from having a peak to becoming linear with a tail, decreases with increasing $L$ indicating scaling behaviour. Only for $\epsilon_{w}=0$ is the form of $P(\rho)$ scale invariant, ie. no peak begins to form as $L$ is increased. In view of this behaviour we believe that $\epsilon_{w}=0$ marks the true critical drying point for the system with LR wall-fluid interaction and that in the thermodynamic limit a liquid peak (indicating partial drying) will occur for all $\epsilon_{w}>0$. The finite-size dependence of the value of $\epsilon_{w}$ for which a peak first occurs feeds through to the observed finite-size shift in the apparent drying point as measured via the contact angle calculation (fig. 9). Recall that both DFT (figs. 3-5) and binding potential calculations (Sec. IV.1) also predict critical drying for $\epsilon_{w}=0$. Note also that for $\epsilon_{w}=0$, $W_{LR}(z)$ in (5) reduces to the hard wall potential which we know leads to complete drying Henderson and van Swol (1985); Oettel (2005).

The form of $P(\rho)$ at criticality, namely a linear part with a tail, is essentially a universal finite-size scaling (FSS) function, albeit a trivial one. The simplicity of its form stands in stark contrast to critical point FSS functions for the order parameter found at bulk critical points (see eg. Wilding (1995)). This difference highlights what we believe to be a fundamental difference between bulk and surface criticality in 3d as probed by simulations, namely that one can approach, but never quite reach a surface critical point.

To elaborate, consider the behaviour approaching the critical drying point $\epsilon_{w}\to\epsilon_{wd}^{+}$. As criticality is neared, the parallel correlation length grows like $\xi_{\parallel}\sim(\epsilon_{w}-\epsilon_{wd})^{-\nu_{\parallel}}$ (see eq. 20 and Sec. VII.3 below). Physically, one can regard $\xi_{\parallel}$ as reflecting the lateral size of ‘bubbles’ of the incipient (vapor) phase that form at the wall. The divergence in the bubble size as $\epsilon_{w}\to\epsilon_{wd}^{+}$ implies that $\xi_{\parallel}$ can grow up to the system size $L$. However, the situation is different for the bubble thickness in the perpendicular direction. This thickness is given by $\xi_{\perp}$, the surface roughness which also, in principle, diverges (albeit logarithmically) as $\epsilon_{w}\to\epsilon_{wd}^{+}$. However, unlike $\xi_{\parallel}$, the broken translational symmetry perpendicular to the walls implies that in 3d simulations $\xi_{\perp}$ is strongly dampened by finite-size effects. General capillary wave arguments e.g. Gelfand and Fisher (1990); Dietrich (1988); Evans (1992) for a single unbinding vapor-liquid interface predict that the surface roughness $\xi_{\perp}\simeq\sqrt{(k_{B}T/2\pi\gamma_{lv})\ln(L/\xi_{b})}$. Thus the surface roughness depends on the finite lateral dimension of the system and owing to the strong $\sqrt{\ln L}$ dampening one expects that for currently accessible system sizes the bubble thickness does not become large on the scale of the particle diameter (or indeed the bulk correlation length $\xi_{b}$). A configurational snapshot of the emerging vapor-liquid interface at the wall (Fig. 11) in which particles are colored according to their distance from the wall, confirms qualitatively this picture. Observing the correlated regions of purple shaded particles lying close to the wall and the green shaded particles further from the wall, we note that there is a large but finite $\xi_{\parallel}$ manifest in the large fractal bubbles of ‘vapor’ which almost span the system in the lateral dimension. However, the perpendicular extent of these bubbles is microscopic, extending only a few particle diameters away from the wall. We discuss the accompanying density profile in the next subsection.

As $\xi_{\parallel}\to L$ for some $\epsilon_{w}>\epsilon_{wd}$, a vapor bubble spans the wall allowing the liquid to unbind and form a free ‘slab’, surrounded by vapor. In terms of the form of $P(\rho)$, this happens when the liquid peak –which corresponds to a local free energy minimum that serves to binds the liquid to the wall– disappears. Recall, however, that the vanishing of the liquid peak occurs at the same value of $\epsilon_{w}$ for which the contact angle measurements predict $\cos(\theta)=-1$. Thus the unbinding process can be viewed as premature drying induced by the finite system size. We shall denote the wall strength for which it occurs by $\epsilon_{wd}(L)$. Since the unbound liquid slab can fluctuate freely away from the wall, premature drying marks a spontaneous loss of the near-critical state at $\epsilon_{w}\approx\epsilon_{wd}(L)$. This state of affairs contrasts starkly with the situation for simulations of bulk critical phenomena where correlations diverge isotropically and critical fluctuations can be sampled right up to criticality. It appears not to have been recognized previously.

The existence of the non-critical state allows one to rationalize the form of $P(\rho)$ at critical drying. Owing to the dampening of capillary fluctuations, the surface of the detached liquid slab is rather sharp and localized and hence the slab thickness (in the $z$-direction) is proportional to $\rho$. Accordingly, the linear decrease of $P(\rho|\epsilon_{w}=0)$ seen at low to moderate densities in Fig. 10 arises simply from the ‘entropic repulsion’ of the slab and the wall: the number of positions for the slab center along the $z$ axis that are allowed by the presence of the wall, varies linearly with slab thickness. The high density tail of $P(\rho)$ on the other hand reflects the free energy cost of pushing the liquid up against the wall, the act of which quenches the parallel density fluctuations. Its $L$ dependence arises –as shown in fig S2. of Evans et al. (2016)– from a constant repulsive pressure on the liquid-vapor interface by the wall, giving rise to a force which scales simply with the wall area $L^{2}$.

The critical wall strength, $\epsilon_{wd}$, is determined most accurately as the largest value of $\epsilon_{w}$ for which $P(\rho)$ assumes an $L$-independent form. This value can differ substantially from $\epsilon_{wd}(L)$. For the LR system we find $\epsilon_{wd}=0$, in agreement with the prediction of binding potential and DFT calculations. However, for a system having SR wall-fluid interactions (as described in Sec. II), binding potential calculations, as well as DFT and GCMC estimates of contact angles (cf. Eq.  (41) and Figs. 5,8) predict that a critical drying transition occurs for a non-zero attractive wall strength. The value of the critical drying point is not known a-priori from a binding potential analysis in this case and hence it is interesting to see whether an accurate determination can be made by examining the $\epsilon_{w}$ and $L$ dependence of $P(\rho)$. Fig. 12 shows that the behaviour mirrors qualitatively that found in the LR case (Fig. 10). The main difference is that on decreasing $\epsilon_{w}$ from large values, a non-zero value of $\epsilon_{w}$ is reached below which $P(\rho)$ exhibits a linear part and a tail for all $L$. On the basis of the arguments given above we take this value to be the critical drying point which we estimate to occur for $\epsilon_{wd}=0.52(2)$.

VII.3 Density and compressibility profiles and a Maxwell relation.

The existence of premature drying identified above implies that the critical limit can be accessed only in the thermodynamic limit $L\to\infty$. This begs the practical question: How closely can one approach criticality for a given $L$ such that estimates of observables are representative of the thermodynamic limit in the near-critical region rather than the non-critical state?

To answer this question we have investigated the near-critical behaviour of the local compressibility and density profiles. Fig. 13(a) compares the form of $\chi(z)$ for a selection of values of $L$ for the LR system at the near-critical wall strength $\epsilon_{w}=0.5$. A dramatic finite-size dependence is apparent. Specifically, as $L$ is increased from small values, $\chi(z)$ decreases strongly, before converging for sufficiently large $L$. Such a decrease is at first sight most surprising because in bulk systems the total compressibility generally increases with system size in the vicinity of a critical point. The origin can be traced to a smaller free energy cost for fluctuations to lower densities compared to those to higher densities, as manifest in the finite-size forms of $\ln P(\rho)$ (Fig. 13(b)). In the partial drying regime, the system occupies states whose densities lie under the liquid peak shown in this figure. However, the shape of the peak is strongly asymmetric, with a tail extending to lower densities which runs smoothly into the plateau associated with the non-critical fluctuations of an unbound liquid slab. The tail reflects the relative ‘softness’ of fluctuations that reduce the density near the wall. Fig. 13(b) also shows that the liquid peak height (as measured from the plateau) grows with increasing $L$ 333For a given $\epsilon_{w}$, the height of the liquid peak measured from the plateau corresponds to the free energy cost per unit wall area of forming a vapor-liquid interface by removing the liquid from the wall, ie. it is proportional to $\gamma_{wv}-\gamma_{wl}+\gamma_{vl}$. . Accordingly, the extent to which the density fluctuations can escape the top of the liquid peak and sample the tail and plateau region decreases with increasing $L$. It follows that for small $L$ the sampling includes contributions from slab fluctuations (and/or their precursors in the tail region) resulting in a spurious enhancement of the compressibility.

Given this insight, it is interesting to reassess the role of finite-size effects in previous simulation studies for fluids near weakly attractive substrates. In Ref. Evans and Wilding (2015), $\chi(z)$ was measured for SPC/E water as a function of the attractive wall strength $\epsilon_{w}$. For small $\epsilon_{w}$, the form of $\chi(z)$ that was observed is similar to that shown for $L=15\sigma$ in Fig. 13(a). In view of the smaller values of $L$ that were attainable for the water model compared to the current LJ system, it is now clear that this compressibility profile was affected by finite-size effects. Remaining with SPC/E water, we note reports of asymmetry in the form of the probability function of the fluctuating density within a subvolume located close to a large hydrophobic solute particle Patel et al. (2010, 2012). The findings shown in Fig. 13(b) help rationalize this observation in terms of finite-size effects.

The upshot of our analysis of finite-size effects is that in order to obtain estimates of observables that are representative of the thermodynamic limit, one must ensure that $L$ is sufficiently large for the prescribed $\epsilon_{w}$ that the liquid peak is high. This in turn requires $\xi_{\parallel}\ll L$. Measurements of $\rho(z)$ and $\chi(z)$ that have been found to be $L$-independent are shown in Fig. 14 for $L=50\sigma$. For this rather large value of $L$, finite-size effects are found to be small provided $\epsilon_{w}\gtrsim 0.3$ (recall that criticality is at $\epsilon_{w}=0$ for the LR system). Within this region we observe clear evidence of strong near-critical fluctuations: the local compressibility near the wall exceeds its bulk value by a factor in excess of $200$, cf. Fig. 14(a). Additionally, a growing drying layer is associated with the density profile of Fig. 14(b). Note, however, that owing to its very weak critical divergence – see (30), the drying layer thickness does not attain more than a few particle diameters, for these values of $L$, even when $\xi_{\parallel}$ is an order of magnitude larger.

Comparison with the DFT results in Fig. 6 are revealing. We note first that in the simulations the density and $\chi(z)$ profiles very close to the walls have different shapes from those obtained in DFT. The latter decrease smoothly to zero as $z\to 0$, reflecting the soft $wf$ repulsion in (5). Since the $wf$ potential employed in simulations is the modified form with its minimum at the hard-wall, $z=0$, the density profiles, for larger values of $\epsilon_{w}$, are increasing as $z\to 0$. However, this difference is not important for small values of $\epsilon_{w}$ when the drying layer has developed. The density profiles then have very similar forms and integrated quantities such as the adsorption should not depend on the fine details of the $wf$ potential. $\chi(z$) goes to a non-zero value at $z=0$ in simulation but the overall variation in shape with $\epsilon_{w}$ is similar to that in DFT. What is striking is that for a similar thickness of drying layer the simulation results yield much larger maxima in $\chi(z$).

As criticality is approached, finite-size effects begin to manifest themselves when $\xi_{\parallel}\simeq L$. We find that the form of $\chi(z)$ is considerably more sensitive to changes in $L$ in this regime than is $\rho(z)$ – an effect that is traceable to the much stronger critical divergence of $\xi_{\parallel}$ compared to that of the drying film thickness (or the adsorption $\Gamma$), cf. Eqs. (21),(30),(33). To probe further the relationship between the two profiles, we have examined within simulation the Maxwell relation eq. (12), that links $\Gamma_{1}\equiv\partial\Gamma/\partial\epsilon_{w}$ to the weighted compressibility $\chi_{1}\equiv\partial\Theta/\partial\mu$. The adsorption $\Gamma(\epsilon_{w})$ was obtained from $\rho(z)$ using (8) allowing access to its numerical derivative $\Gamma_{1}$ in (12). A fit to the latter is shown in Fig. 15, where it is compared with our measurements of $\chi_{1}$. One observes that at large $\epsilon_{w}$, far from criticality, there are small discrepancies between the two quantities. These we attribute to the fact that when the drying layer is thin $\rho(z)$ contains more structure and numerical integration is less accurate. On moving to smaller $\epsilon_{w}$, for which $\xi_{\parallel}$ is large (but still small compared to $L$), there is excellent agreement between $\Gamma_{1}$ and $\chi_{1}$ – a finding that verifies our numerics. However, as the wall strength $\epsilon_{w}\approx 0.3$ is reached, a significant discrepancy starts to appear. Given the differing sensitivities of $\chi(z)$ and $\rho(z)$ to finite-size effects noted above, we speculate that this discrepancy serves as an indicator that the limit $\xi\approx L$ has been reached and finite-size effects are significant. The extent to which the Maxwell relation holds therefore appears to serve as a useful tool for diagnosing when finite size effects are significant.

VII.4 Estimates of $\nu_{\parallel}$

In Sec VII.2 the form of $P(\rho)$ was examined for state points near criticality. Precisely at criticality $P(\rho)$ comprises a linear part and a tail. The linear part occurs at low to moderate densities and arises from the entropic repulsion of the wall to the unbound liquid slab. The high density tail arises from the free energy cost of pushing the slab up against the wall. Neither of these phenomena is directly associated with criticality, and thus one cannot expect $P(\rho)$ to exhibit non-trivial finite-size scaling (FSS) behavior as a whole. Rather, the signature of near critical fluctuations is manifest in the density range where the liquid is still (weakly) bound to the wall but exhibits strong parallel density fluctuations. This correspond to the liquid peak in Fig. 10, the height of which depends on $\xi_{\parallel}$ and vanishes when $\xi_{\parallel}\approx L$ allowing the liquid slab to unbind from the wall. Simple FSS dictates that this vanishing occurs not at $\epsilon_{wd}$ but at the larger effective value $\epsilon_{wd}(L)=\epsilon_{wd}+aL^{-1/\nu_{\parallel}}$.

We have determined the value of $\nu_{\parallel}$ via the anticipated FSS $\epsilon_{wd}(L)\sim L^{-1/\nu_{\parallel}}$; $\epsilon_{wd}=0$ for the LR case. For a number of choices of $L$ we measured $\epsilon_{wd}(L)$ accurately (via histogram extrapolation) from the vanishing of the liquid peak of $P(\rho)$ (cf. Fig. 10). As Fig. 16 shows, we do indeed see power law scaling, from which we can extract an estimate of $\nu_{\parallel}$ –see key. Interestingly, however, this estimate exceeds the prediction $\nu_{\parallel}=0.5$ of mean field and RG theories (see Secs. IV.1, IV.2) by more than a factor of two and additionally appears to show a clear dependence on the temperature.

A further independent estimate of $\nu_{\parallel}$ can be obtained from the growth in the maximum of $\chi(z)$ on the approach to critical drying. On theoretical grounds one expects (cf. Eq. (21)) that $\chi_{\rm max}\sim(\epsilon_{w}-\epsilon_{wd})^{-2\nu_{\parallel}}$. The data of Fig. 17, when fitted to a power law, yield an estimate of $\nu_{\parallel}$ (see key of Fig. 17) that is again more than twice the theoretical prediction and (in common with the finding for $\epsilon_{wd}(L)$) seem to demonstrate a clear temperature dependence. We shall return to these findings in Sec. VIII.

VII.5 Critical wetting

So far we have focused exclusively on critical drying. However the contact angle measurements of Figs. 5 and 8 suggest that for a SR wall-fluid potential (4), critical wetting occurs. It is therefore of interest to ask how the associated phenomenology compares with that of critical drying. To this end we have performed simulation studies of the $L$-dependence of $P(\rho)$ in the neighborhood of the wetting point suggested by the contact angle measurements. Our results are shown in Fig. 18 and exhibit closely analogous behaviour to that seen for critical drying in the SR case (Fig. 12). The main difference is that here it is the vapor phase that is metastable, displaying a peak that decays smoothly into a plateau before turning into a bulge, which decays further with increasing $\epsilon_{w}$ until all that remains is a linear part and a tail. As for drying, the effect of increasing $L$ is to increase the height of any peak or the strength of any bulge. However once critical wetting is reached, $P(\rho)$ exhibits a linear part and a tail for all $L$. The smallest value of $\epsilon_{w}$ for which this occurs therefore marks the critical wetting point. For our SR system this appears to occur at $\epsilon_{ww}\approx 4.2(2)$. We note that this is very different from the estimate of $\epsilon_{ww}(L)=3.7(1)$ which emerges from the contact angle measurements for $L=15\sigma$ and which corresponds to the point at which the vapor peak disappears into the plateau. Clearly, therefore, an analysis of finite-size effects is indispensable when seeking to obtain accurate estimates of critical wetting points, as indeed it is for critical drying.

In Fig. 19 we show results for the normalised compressibility profile $\chi(z)/\chi_{b}$ and density profile $\rho(z)/\rho_{b}$ on the approach to critical wetting. For the system size $L=50\sigma$ that we used, the estimates were unaffected by finite-size effects for $\epsilon_{w}\lesssim 3.6$. In common with the results for critical drying (cf. Fig. 14), one sees a very large relative compressibility near the walls. The local compressibility in the liquid-like layers near the walls is up to 500 times that of the bulk vapor at the same coexistence chemical potential. The main difference between critical drying and wetting is that owing to the very strong SR wall attraction in the case of wetting, packing effects occur in the density profile near the wall, and these modulate the compressibility profile as well. An interesting feature of the density profile, Fig. 19(b), is that a very dense single layer of particles is strongly adsorbed on the walls. Given the very short range ($0.5\sigma$) of the square well $wf$ potential, particles beyond this layer do not interact with the wall directly, rather they experience an effective confining potential that stems from the dense first layer, and whose form is that of the truncated LJ $ff$ potential. This implies the system is still governed by SR interactions, allowing critical wetting to occur.

VIII Discussion and conclusions

In summary, we have investigated the properties of a fluid having truncated particle-particle ($ff$) interactions which is confined between smooth planar walls. Two forms of the wall-fluid interaction potential $W(z)$ have been considered: the long ranged (LR) case (5) in which $W(z)$ exhibits power law decay and the short ranged (SR) case (4) of a square well potential. Clear evidence has been obtained that the character of the wetting and drying transitions is sensitive to this range. In particular we find from simulation and DFT that for the LR case wetting is first order, while for the square well wall it is critical 444In a further investigation whose results we do not show here, we have studied the effect of truncating the $9$-$3$ wf potentials at different cut-off distances $z_{cw}$ in (5). As expected, we find that the wetting transition is strongly first order for large values of the cut-off. It remains first order as $z_{cw}$ is decreased to about $2\sigma$. For smaller values the transition appears to be critical, as in the square well case.. By contrast, drying is critical for both the SR and LR cases. For the latter case drying occurs at exactly zero attractive wall strength ie. for a hard wall. Of course it is well known that in this limit complete drying occurs for all $T<T_{c}$ Henderson and van Swol (1985); Oettel (2005). What is remarkable, however, is that the transition is critical and is predicted to occur precisely at $\epsilon_{w}=0$ for all LR (power law) forms of $W(z)$.

Knowledge of the exact location of a surface phase transition in 3d provides a unique opportunity to study a surface critical point free from uncertainty regarding its location (a problem that has previously plagued Ising model studies of critical wetting Bryk and Binder (2013)). We have obtained simulation estimates of the surface critical exponent $\nu_{\parallel}$ in the LR case and compared these with the predictions of mean field theory and a linear renormalization group calculation, finding our estimate to be over twice the predicted value. Furthermore, the exponent estimates exhibit a clear temperature dependence – a feature which is also at variance with the theoretical predictions. Given that $d=3$ is the upper critical dimension for this system, at which mean field theory for the exponents is expected to hold, these findings are unexpected.

A possible reason for the discrepancy is finite-size effects, the nature of which we have sought to elucidate. Finite-size effects for surface criticality have been investigated previously in the context of critical wetting in the Ising model with SR forces (where the special symmetry implies that wetting and drying are equivalent). A finite-size scaling ansatz proposed by Binder and co-workers Albano and Binder (2012); Bryk and Binder (2013) assumes that critical wetting in 3d can be treated on the same footing as a bulk transition in which two independent correlation lengths diverge in orthogonal directions. However, our results call this assumption into question. We find that while $\xi_{\parallel}$ can grow as large as the wall dimension $L$, $\xi_{\perp}$ is, by contrast, heavily damped for finite $L$ and remains microscopic for all accessible system sizes (recall that capillary wave theory predicts $\xi_{\perp}\sim\sqrt{\ln L}$ in $d=3$). Moreover, as $\xi_{\parallel}\to L$ the system spontaneously exits the near-critical state due to the premature unbinding of the metastable phase. This state of affairs differs qualitatively from that for bulk criticality where critical fluctuations can be measured right up to the critical point, together with quantities such as cumulants of the order parameter distribution. Clearly fresh and bespoke FSS approaches are need for dealing with surface criticality in fluid systems.

Given the heavy dampening of $\xi_{\perp}$, the surface critical behaviour observed in simulations appears to be controlled by the single diverging lengthscale $\xi_{\parallel}$. It is therefore tempting to speculate that effective critical surface behaviour is $2d$ Ising like in character. Indeed our measured values of $\nu_{\parallel}$ are much closer to $\nu_{\parallel}=1$ than they are to the theoretical prediction $\nu_{\parallel}=1/2$. Another issue to note is that owing to its slow (logarithmic) divergence as $\epsilon_{w}\to 0$, the thickness of the drying layer does not exceed a few particle diameters at the state point closest to criticality for which we can attain the thermodynamic limit. Of course, for very large system sizes one might expect to observe a crossover to the mean-field and RG predictions of Sec. IV; however no hints of such a crossover are visible in our results, despite the investment of substantial computational resources to study large systems. In our view, if such a crossover exists, there seems little hope of observing it in the foreseeable future.

Our finding that premature drying occurs when $\xi_{\parallel}\simeq L$ helps to explain a longstanding controversy in the literature van Swol and Henderson (1989); Henderson and van Swol (1990); van Swol and Henderson (1991); Nijmeijer et al. (1992, 1991); Henderson et al. (1992); Bruin et al. (1995); Bruin (1998) concerning the order of the drying transition. On the basis of MD simulation, van Swol and Henderson asserted that the drying transition for a square well fluid with SR wall-fluid interactions is first order in character. The evidence for this was that at small attractive wall strength an abrupt change was observed in the density profile $\rho(z)$ from a liquid-like profile to a gas-like one. Given the insights provided by the present work, one can see that rather than being associated with a first order drying transition, this phenomenology arises from the premature unbinding of the liquid slab which then diffuses away from wall.

One of the motivations for the present work was to contribute to ongoing attempts to understand the correlation between the structure of water near a hydrophobic surface and the value of the contact angle. Experimental studies have reported a region of depleted density in the close proximity of hydrophobic surfaces Mezger et al. (2006); Ocko et al. (2008); Mezger et al. (2010); Chattopadhyay et al. (2010, 2011); Uysal et al. (2013), while simulation studies report a growth in density fluctuations near the surface as the contact angle increases Chandler (2007); Patel et al. (2010, 2012); Mittal and Hummer (2010); Willard and Chandler (2014); Jamadagni et al. (2011). As we have shown previously Evans and Wilding (2015), both phenomena in water can be accounted for if hydrophobic substrates can be associated with the approach to the critical drying point. Pertinent is the form of the wall-oxygen potential. In the GCMC simulations of Ref. Kumar and Errington (2013a) this was non-truncated $9$-$3$, equivalent to (5), and the water model was SPC/E which is SR. On the basis of the arguments presented here we would predict critical drying of water in the limit where the attraction strength $\epsilon_{w}$ vanishes. Indeed Kumar and Errington (see Fig. 9 of Kumar and Errington (2013a)) appeared to observe $\cos(\theta)$ approaching $-1$ (tangentially) for very small $\epsilon_{w}$. In our own study of SPC/E Evans and Wilding (2015) the $9$-$3$ wall-oxygen potential was truncated at a large distance $15$Å and similar behaviour was found but we could not identify the drying point accurately. The MD study Willard and Chandler (2014) employed a wall-oxygen potential with a non-truncated $z^{-6}$ tail. Once again we would predict critical drying in the limit of vanishing attraction. It is likely that the strong density fluctuations close to the wall that are observed in Willard and Chandler (2014) are associated with the approach to the critical region; the authors consider very weak $wf$ attraction but their estimates of $\cos(\theta)$ are not sufficiently accurate to address the location of the transition point. The present study of a generic fluid model serves to emphasize that there is nothing special about water with respect to these phenomena. Indeed it seems that critical drying should be expected for all liquids provided the substrate is sufficiently weak and this will be associated with a growing drying layer and enhanced density fluctuations near the substrate. Furthermore it can be expected that the presence of the critical drying point is to be felt throughout the hydrophobic (partial drying) regime i.e. not just in the limit $\cos(\theta)=-1$.

Finally we briefly point out a number of avenues for future work. While the present study has focused on critical drying in the particular case of SR $ff$ with LR $wf$ because here the critical point is known exactly, it would be interesting to see if one can determine critical properties for drying and wetting in a SR system such as a square well wall (4). Our results (Figs. 12,18) strongly suggest that the hallmark of surface criticality is the same as for the LR case, namely a scale invariant limit in which $P(\rho)$ exhibits a linear part with a tail. It remains to be seen whether application of this method for locating criticality provides sufficient accuracy to allow reliable estimates of exponents, which are predicted to depend on the form of the $wf$ and $ff$ potentials-see Sec. III.2. Beyond this, it would also be interesting to investigate the case most pertinent to real systems, where dispersion forces ensure that both $wf$ and $ff$ are LR and critical drying can occur at non-zero wall strength (whose value is known exactly in terms of the coexisting densities Stewart and Evans (2005)) and the critical exponents are predicted to take large values, e.g. $\nu_{\parallel}=5/2$ Dietrich (1988).

Acknowledgements.

We have benefited from helpful discussions with A.J. Archer and A.O. Parry. R.E. acknowledges Leverhulme Trust grant EM-2016-031. This research made use of the Balena High Performance Computing Service at the University of Bath.