Putting the micro into the macro: A molecularly-augmented hydrodynamic model of dynamic wetting applied to flow instabilities during forced dewetting

We report a molecularly-augmented continuum-based computational model of dynamic wetting and apply it to the displacement of an externally-driven liquid plug between two partially-wetted parallel plates. The results closely follow those obtained in a recent molecular-dynamics (MD) study of the same problem Toledano (2021) which we use as a benchmark. We are able to interpret the maximum speed of dewetting $U^*_{\mathrm{crit}}$ as a fold bifurcation in the steady phase diagram and show that its dependence on the true contact angle $\theta_{\mathrm{cl}}$ is quantitatively similar to that found using MD. A key feature of the model is that the contact angle is dependent on the speed of the contact line, with $\theta_{\mathrm{cl}}$ emerging as part of the solution. The model enables us to study the formation of a thin film at dewetting speeds $U^*>U^*_{\mathrm{crit}}$ across a range of length scales, including those that are computationally prohibitive to MD simulations. We show that the thickness of the film scales linearly with the channel width and is only weakly dependent on the capillary number. This work provides a link between matched asymptotic techniques (valid for larger geometries) and MD simulations (valid for smaller geometries). In addition, we find that the apparent angle, the experimentally visible contact angle at the fold bifurcation, is not zero. This is in contrast to the prediction of conventional treatments based on the lubrication model of flow near the contact line, but consistent with experiment.


I. INTRODUCTION
Dynamic wetting, the process by which a liquid wets a solid surface, is an important phenomenon that underpins a wide range of both industrial and natural processes, including microfluidics [68], liquid coating and printing operations [74], petroleum recovery [33], plant protection [52], ground water hydrology [5] and biological processes [4]. As such, it presents a multiscale problem. Whilst its origin is at the microscopic scale of the moving contact line, it influences outcomes at very much larger scales. However, despite this importance, and consequent research over many decades, there remain fundamental questions about the physics involved and, in particular, the role of solid-liquid interactions at the moving contact line [1,2,58].
In wetting studies, solid-liquid interactions are usually quantified in terms of the angle of contact between the liquid and the solid, and its proper description has attracted much attention [2,8,19,59]. From hydrostatic and hydrodynamic perspectives, this boundary condition is crucial, as it dictates the shape of the liquid volume. The way it changes in response to movement of the contact line across the solid surface is, therefore, fundamental to our ability to predict wetting outcomes. Nevertheless, the description of the true contact angle at a moving contact line remains hotly debated [1,2,58].
In continuum models the true contact angle, measured by the tangent of the interface at the solid (see figure 1), has to be specified in order to solve the governing equations and is usually considered to be constant and equal to the equilibrium value. The observed dynamics of the apparent contact angle (i.e. the one seen experimentally [75]) is attributed to the 'viscous bending' of the interface: this is the so-called 'Hydrodynamic' or Cox-Voinov formulation [17,73]. However, according to the molecularkinetic theory of wetting (MKT) [7,10] and the interface formation model [59], the true contact angle varies and is dependent on the velocity of the contact line.
Here, we will show that viscous bending alone is insufficient to capture the effects seen in molecular simulations, where the velocity dependence of the actual contact angle is observed. Therefore, we develop a new combined approach based on the Navier-Stokes continuum paradigm combined with the MKT, (whose formulation is far simpler than the interface formation model, despite the latter's attractive features), and focus it on the canonical dynamic wetting problem of a liquid plug propagating  [29], with permission from Elsevier. The panels a) to e) show the liquid plug as a the force F * 0 (in the article asterisks were not used to denote dimensional quantities) becomes successively larger and eventually exceeds the critical value (panels d) and e)) where a thin-film begins to develop. In this MD simulation, the external phase is a vacuum. through a channel. In particular, in order to allow for unambiguous comparisons to the results of molecular dynamics on a comparable system, we identify a critical speed at which a flow bifurcation occurs and a thin film is formed.
The study of dynamic wetting using molecular simulations has a long history, see review articles [18] and [42]; but here we focus on a recent paper, [29], that examines both wetting transitions and the behaviour of the contact angle. In this study, large-scale molecular dynamics (MD) is utilised to explore the steady displacement of a water-like liquid plug between two molecularly-smooth solid plates under the influence of an external driving force F * 0 (see figure 1 for the geometry). The study used a coarse-grained model of water and an atomistic Lennard-Jones model for the solid plates. The general behaviour observed as F * 0 was increased, and hence the liquid-plug's speed was raised, is depicted in figure 2. Notably, it was reported that both the 'true', dynamic contact angle at the contact line, θ cl , and a larger-scale 'apparent' angle, θ app , are dependent on the contact-line velocity U * cl for the receding and advancing interfaces. Henceforth, unless otherwise stated, when we refer to a 'contact angle' we mean the 'true' contact angle. We also note that quantities labelled with an asterisk correspond to dimensional physical quantities and those without to dimensionless quantities.
In the MD study, the apparent angle was measured at the system scale by a method that mimics typical measurements of it in macroscopic experiments, where the precise details of the true contact angle's dynamics remain hidden, as they occur on such small length scales [8,22,37]. By varying the solid-liquid affinity (i.e. the solid's wettability), it was possible to investigate the influence of the equilibrium contact angle θ 0 on the results. For all θ 0 , θ cl was found to be velocity dependent in a manner consistent with the MKT of dynamic wetting [7,10]. However, θ app diverged from θ cl as F * 0 was increased, especially at the receding contact line, in a way that closely followed the Voinov equation [73]: where Ca = µ * U * cl /γ * is the capillary number based on the contact-line speed U * cl , dynamic viscosity µ * and surface tension γ * , and L * and L * m are suitably-chosen macroscopic and microscopic length scales. For each θ 0 , there was a critical receding contact-line velocity U * crit and contact angle θ crit at which θ app became small and the receding meniscus deposited a liquid film on the plates. This value could then be used in (1), assuming θ app ≈ 0, to fix L * /L * m and hence reliably predict θ cl at both the advancing and receding contact lines. This result is significant, as θ cl is not usually experimentally accessible and the fact that it varies with U * cl poses questions for hydrodynamic interpretations of dynamic wetting. The result also shows that the critical condition for film deposition encodes crucial information about the hydrodynamics.
The existence of a critical wetting speed has been investigated thoroughly using hydrodynamic models in a range of geometries, including those associated with coating flows [44], and plate withdrawal [63], among others. In [41], both receding and advancing contact line problems were investigated for a coating flow and the stability of the solutions near the critical speed was quantified using a dynamical systems method. Here, our focus will be on the receding contact line, as this is where the first bifurcation will occur. Previous studies have shown that as Ca increases, the receding contact line will attain a steady state provided Ca < Ca crit , where Ca crit is a critical capillary number that is a function, at the very least, of θ cl , the slip length and the viscosity ratio of the liquid and gas phases [17,24,41,61,62], but, if Ca > Ca crit , a thin-film develops with thickness dependent on Ca [41,62]. Using a lubrication model, Ca crit can be approximated when the slip-length is small relative to the film height [25], by considering a small-Ca asymptotic analysis and using the key assumption that θ app = 0 at the critical point. However, in a nano-geometry, as considered here, we will see that this assumption is not valid and the resulting small-Ca asymptotic analysis does not extend to this regime.
In this paper we will develop a hydrodynamic model based on the Navier-Stokes paradigm to calculate steady states and transient behaviour of the liquid plug scenario considered in [29]. An essential aspect of this model is that the true angle, θ cl , has to be specified at the junction of the liquid, gas and solid phases. In many previous studies where a Navier-Stokes model is used (see, for example [40, 47-50, 66, 71, 72]) θ cl is assumed to be constant; but motivated by the results of [29] we relax this assumption and adopt a model that determines θ cl as a function of Ca and the static contact angle θ 0 based on the MKT. Notably, the model remains hydrodynamic throughout, in contrast, for example, to [35], and the molecular-augmentation comes entirely through the contact-angle formula. This approach was also considered for macroscopic flows in [21], where the static contact angle, θ 0 and slip length are independent parameters. However, motivated again by [12] and [28], we will make use of a correlation between slip length and θ 0 that reduces the number of parameters that are required. This correlation is based on an assumption, borne out by MD simulations, that the mechanism of slip between a liquid and a solid is the same across all parts of the solid-liquid interface, including the contact line.
The article is ordered as follows. In § II we describe the system of equations used to model the liquid plug based on the Navier-Stokes (NS) equations. In addition to the NS equations, in § III we discuss asymptotic results, based on a Quasi-Parallel (QP) lubrication approach adapted from [25], that will be relevant here. By calculating numerical solutions of the governing equations using a finite-element framework, we will then show in § IV that the critical speed of wetting for the entire liquid plug is dependent on the receding contact line and not influenced by the advancing contact line. We will also discuss the method for calculating the apparent angle. Next, in § V we will show how augmenting the NS equations with an MKT variable-angle constraint predicts the existence of a critical Ca, and that as the wettability is varied the values of Ca crit match favourably with the MD data in [29], in contrast to the predictions of the fixed-angle model. In addition, we demonstrate how θ cl and θ app vary with the slip length and thus provide an estimate of L * /L * m for the liquid plug system, which shows excellent agreement with the MD simulations. Further, in § VI we examine time-dependent behaviour when Ca > Ca crit so that a thin film develops, whose height obeys a Landau-Levich law. Finally, in § VII, having validated the system in the liquid nano plug geometry, we exploit our computational framework to explore larger-scale systems, which are beyond the scope of MD simulations. By examining systems where the physical size is orders of magnitude larger than the nano-channel studied in [29], we will show that the dimensionless thickness of the film remains constant, for a fixed Ca.

II. MOLECULAR-AUGMENTED HYDRODYNAMIC MODEL
We will now describe the hydrodynamic model. We shall discuss the full system, based on the Navier-Stokes equations and then describe the two different system formulations; the pressure-driven problem and the force-driven problem, as well as the numerical method and the different computational domains.

A. Fully Nonlinear System
To mimic the molecular simulations, we model the liquid-bridge system as a two-dimensional flow between two parallel plates as illustrated in figure 1 and detailed in figure 3(a). A finite liquid region fills the channel bounded by two rigid plates that are separated by a distance H * . We solve in a frame of reference that moves with the plug, with the walls moving with velocity U * wall . The exact formulation depends on the domain and problem that we consider (i.e. pressure-driven or body-force driven), details of which we discuss later. We nondimensionalise all lengths using the half-height, H * /2, all velocities using U * wall , all pressures by µ * U * wall /(H * /2) all timescales by (H * /2)/U * wall and the body force by µ * U * wall /(H * /2) 2 . As in other studies [47-50, 64, 65, 67, 71, 72], we apply the Stokes-flow approximation, (2)-(3), so that the Reynolds number, ρ * U * wall H * /µ * , is assumed to be negligibly small; simple estimates confirm that this is appropriate for the nano-system. We neglect the influence of gravity and assume that that the gas phase can be modelled as a vacuum (as seen in figure 2, there are no molecules in the gas phase). A typical computational domain is shown in figure 3(a). On the moving wall (y = 0) we apply a Navier-slip condition, (4), and, therefore, introduce a dimensionless slip-length, λ . The MD simulations in figure 2 indicate the flow is symmetric around the centreline of the channel and hence we introduce a symmetry wall at y = 1, labelled Γ 2 , where we set the vertical component of velocity to be zero, apply zero tangential stress, and let the horizontal velocity be determined as part of the solution. As well as the fluid velocity field, u(t, x), and pressure, p(t, x), which depend on the dimensionless time, t, and the position, x of the interfaces, denoted R adv = (x a (t, s), y a (t, s)) and R rec = (x r (t, s), y r (t, s)) respectively, are also unknowns in the problem and functions of t and the arclength, s, as measured from the contact point. These are found using dynamic and kinematic conditions on both free surfaces. The governing equations and boundary conditions then become Dynamic Condition (Rec. and Adv.), where n, and t are the vectors normal and tangential, respectively, to the appropriate boundaries denoted Γ i , and κ is the curvature of the corresponding interface. The plate speed U = (U wall , 0) T and λ = λ * /(H * /2) is the dimensionless slip length. The body force is F = (F, 0) T , where F is a set constant. The stress tensor τ is defined as where I is the identity matrix. We shall refer to the system described by (2)-(9) as the 'Half Liquid-Plug' problem.

B. Contact-Angle Models
The system is not well-posed unless a contact angle is specified between the free-surfaces and the horizontal plates. For the symmetry boundary we set θ (s = L) = π/2, but the dynamic contact angle, θ cl , can be freely chosen and depends on the wettability of the solid.
The simplest approach is to specify a constant equilibrium contact angle, i.e. θ cl = Const. Constant Angle Equation (11) However, as is well known, the molecular-kinetic theory predicts θ cl to be dependent on the speed of the contact line. As shown in Appendix A, for the receding contact lines of interest here, this dependence may be written in the linearised form where δ is a dimensionless parameter that corresponds to the width of the three-phase zone, i.e. the contact line viewed at the molecular scale, and Ca is the relative velocity of the contact line to the wall speed, i.e.
where e x is a unit vector in the x direction. Equation (12) is the linearised form of the theory, which will be valid for the system considered in this article. Furthermore, it has long been recognised that a relationship must exist between the slip length and the equilibrium contact angle, e.g., [3,53,69]. Here, motivated by the results of [29] and the theory described in Appendix A,we consider the relationship where a and b are fitting parameters and λ * MD is the physical slip-length derived from the MD data in [29]. We will assume that λ * MD is independent of the physical channel height, H * , and therefore in our non-dimensionalisation λ = 2λ * MD /H * . To investigate the nano-channel used in [29], where H * = 20.2nm, the different values of θ 0 will yield dimensionless slip-lengths in the range λ ∼ 0.02 to 0.2. Alternatively, as we will show in § VII, by varying λ , and keeping θ 0 fixed, we can investigate the effects of varying the physical channel height H * to larger systems. We emphasise that there is a one-to-one correspondence between λ * MD and θ 0 , and hence we are free to prescribe either quantity and use (14) to determine the other. In physical experiments it is more practical to find θ 0 , which can readily be measured, and then determine λ * MD , which is more difficult to measure experimentally. Figure 4 shows the fit of this function to the MD data from [29]. We call the system of equations described in (2)-(9) augmented with the constant angle formula, (11), the Constant-Angle (CA) model, while when augmented with the variable angle formula, (12) and (14), we call it the Variable-Angle (VA) model. FIG. 4. The slip length dependence on the static angle. The markers are the MD data obtained from [29] and the solid line is the curve fit given from (14) with b = −2.342 and a = 5.656 × 10 −9 .

C. Pressure-driven and Force-driven Problems
We now discuss the two different types of problems, i.e. the pressure-driven and force-driven problems. The MD simulations in [29] is a force-driven problem, but pressure-driven problems are relevant in many practical situations, for example, coating flows, [48].
In Pouseille flow, without a free-surface, it is easy to show that a pressure-driven problem can be equivalent to a force-driven one. With a free-surface however, this is not true and each case has to be considered separately. For steady calculations, in both types of problems, the set of equations are ill-posed unless we specify the volume of the liquid-plug. To remove this issue, for pressure-driven flow, we impose a normal stress on Γ 1 : and let the value of p out be determined implicitly by a condition on the overall volume of the liquid-plug, which corresponds to the computational area of the domain. We set F = 0 and solve in a frame of reference that moves such that the walls are non-stationary in the translating frame (U wall = −1). In contrast, for the steady force-driven problem, we set U wall = −1 and p out = 0, but now let Ca be determined implicitly by a volume constraint. In both cases, to overcome the translational invariance, we also have to pin a point on the boundary which depends on the domain of the problem (as discussed below).
For time-dependent problems, whether pressure-driven or force-driven, the volume constraint is unnecessary, as equation (3) ensures that volume is conserved. Instead, we impose a position constraint for the reduced domains (see below) that determines p out or Ca, depending on the problem.

D. Numerical Method
The complete system of equations are discretised and solved using a finite-element method and the open-source oomph-lib package [36], as described in [41]. An unstructured triangular mesh is used which is treated as a pseudo-elastic body, so that changes to the unknown free-surface can be facilitated and the mesh can adapt to capture regions of high velocity or pressure gradients, for example near the contact point. We use a ZZ error estimator, which measures the continuity of the rate of strain in each element, to identify elements that require refinement or unrefinement [77]. As a typical example for the time-dependent calculations, with Ca = 0.5 and λ = 0.01, elemental areas range from ∼ 10 −3 − 10 −1 to accommodate a maximum ZZ error of 10 −3 .

E. 'Half' and 'Quarter' Domains
We can simplify the complexity of the half liquid-plug domain further by solving in two separate 'quarter' domains, each having only one free surface; thus significantly reducing the number of triangular elements required, see figures 3(b) and (c).
To facilitate this, we replace the free-surface (and corresponding dynamic and kinematic boundary conditions) at one end of the computational domain (Γ 1 for the advancing contact line and Γ 3 for the receding contact line) with an imposed normal stress (i.e. equation (15)) and parallel flow condition In these quarter domains, the imposed pressure, p out , is determined implicitly by ensuring the volume per unit length (i.e. the area) of the liquid domain is constant (as in the half liquid-plug domain), so that in each quarter domain we are solving in a frame of reference with a fixed volume. We note that in the quarter domains we impose (15) and (16) in both steady and timedependent calculations. These quarter domain simulations will be referred to as the 'Receding Contact Line' and 'Advancing Contact Line' domains (see figures 3(b) and (c)), respectively, abbreviated to RCL and ACL in the rest of the paper. The origin is different in each of these domains and corresponds to the pinned position of the steady and time-dependent problems. To illustrate the benefit of this reduction, the number of elements required in the computation of figure 5 for the half-liquid plug is ∼ 4000 − 8000, but the number for the RCL is ∼ 400. The reason for the ∼ 90% reduction in elements is because the the pressure gradients are not as severe near the receding contact line, compared to the advancing one. In the next section, we shall show that the dynamics of the whole system and the prediction of a critical Ca are dominated by the receding contact line. Thus, computation of the full half liquid-plug problem, where both the advancing and receding interface are calculated, is not necessary in order to to find the first flow bifurcation and is computationally inefficient when compared to the reduced RCL domain.
We emphasise that the main aim of this study is to investigate the VA model, and not to make a thorough investigation of the differences between pressure-driven and force-driven flow, as the VA model can be applied independently of the problem. Thus, in the results that follow, we mainly consider pressure-driven flow, except when we make a direct comparison with the MD simulations. For the latter, we present force-driven results, as clearly specified. In addition, as we shall show in § IV A, the choice of domain is also independent of the problem (i.e. pressure-driven or force-driven), so we choose the domain that is most important to the flow-bifurcation and which is also the simplest, computationally.

III. REDUCED GOVERNING EQUATIONS (QUASI-PARALLEL SYSTEM)
We shall now discuss a reduced evolution PDE model, the so-called 'Quasi-Parallel' system, before finally obtaining asymptotic results that will help predict the value of Ca crit .
For the receding contact line, the flow near the contact line is approximately parallel, c.f. figure 12, and we can exploit this to reduce the Navier-Stokes equations to a simpler system which requires unknowns only on the fluid interface. As well as the parallel-flow assumption, we assume the horizontal coordinate is approximately the arclength, i.e. x ≈ s, so the full expression for the curvature can be used, not the linearized form as used in conventional lubrication models (see, for example [25]) and long-wave models (see, for example [60]).
Following [39], [56] and [70], we let θ be the angle the interface makes to the horizontal (see figure 1), h = y r (s) be the height of the interface and s be the arclength coordinate measured from the contact line. Using conservation of mass and the kinematic condition on the free-surface, the governing equation for the fluid pressure gradient, ∂ p ∂ s , may be written as [62] ∂ h ∂t where U · e x = −1 for the RCL. The unknown pressure gradient, ∂ p ∂ s , can then be expressed in terms of the exact curvature by differentiating the normal stress balance w.r.t s, i.e.
In order to solve (18), we require two conditions on θ . At the contact line, s = 0, we implement equation (12): and at the symmetry wall, s = L, we set θ = π/2, where L is the overall length of the interface. The shape of the interface can then be recovered by solving Each of these equations requires a single condition, so we set x(s = 0) = y(s = 0) = 0 (choosing the contact line to be at the origin). Finally we note that the length of the interface, L, and hence the size of the domain, s, are not known a priori. To determine L we scale the independent variable s, so that ξ = Ls and ξ = [0, 1]. The total length of the interface, L, can then be determined by the additional constraint that To solve this system of equations, we choose to discretise the spatial derivatives using finite-differences and then the system of equations are solved numerically using Newton's method. We remark that we exclusively concentrate on the steady results of the pressure-driven QP system and do not solve the time-dependent problem, as this is better suited to the full nonlinear system.

A. Asymptotics
We now briefly describe and adapt the analysis of [14] and [25] to find an asymptotic expression for Ca crit . We will not repeat their analysis except for the parts where it differs from the situation we examine here. In both of these previous works gravitational effects are included and the liquid domain is unconfined, whereas we neglect gravity and the system is confined. They also considered only steady solutions, so that time-derivatives in the problem can be ignored.
The matched asymptotics methodology of [14] and [25] is to determine an inner solution, say h = h inner , for small Ca, that is valid close to the contact line, i.e. when s/λ ∼ O(1), and an outer solution, h = h outer say, that is valid far away from the contact line, i.e. when s ∼ O(1). To determine an unknown constant in the outer solution, the inner and outer solutions have to match in a crossover region. This matching procedure yields an equation, for an arbitrary unknown, θ app , which is the angle the outer interface makes with the horizontal. For the ACL domain, θ app is finite for all values of Ca, and thus, in these asymptotic limits at least, there is no critical point for the ACL domain. However, in the RCL domain θ app can only be calculated up to a critical value of Ca, this value being interpreted as Ca crit .
In our problem, for a confined geometry and in the absence of gravity, the inner region analysis near the contact line is identical to the case considered in [14] and [25]. In [25] the effects of gravity are present at first-order in the outer solution. The outer solution, for small Ca, is found by expanding the unknowns as a power series in Ca. We can write a leading-order outer solution of (17) to (20) as where κ s = (π − 2θ app )/2L is the curvature, X rec is the meniscus rise (c.f. figure 5) and θ app is an undetermined constant. The outer solution in (22) describes a sector of a circle with centre (X rec − r, 1) and r = 1/κ s that makes an angle θ app to the horizontal at s = 0. Examining the geometry (see dashed curve in figure 6(a)) gives κ s = sin(π/2 − θ app ), so as θ app → 0, κ s → 1 and hence the interface is a semi-circle of radius 1. In particular we have that When the outer solution, described in (22), is matched to the inner solution, as described in [25], we obtain an expression for θ app where Ai(z) is an Airy function of the first kind. We note this is the same for an unconfined geometry with gravity as considered in [25]. In addition, Ca crit satisfies the same expression as in [25] but has a factor of 2 1/3 in the denominator of the logarithm term, i.e.
Ca crit = θ 3 These results are valid for a constant contact angle model but we can easily extend them to our variable angle formula by expanding (12) in powers of Ca, for Ca ≪ 1, and then θ cl in the above expression is just the static angle θ 0 . However, we will show that the full expression for θ cl in (12) will be required in the formula (25) for it to compare favourably with the numerical results. In the sections that follow, the expressions in (23), (24) and (25) will be compared with the numerical solutions.

IV. SYSTEM MEASUREMENTS AND PARAMETERS
In this section we describe the methods used to determine Ca crit from the numerical calculations of the fully nonlinear and quasi-parallel systems. In addition we also discuss, in detail, the methodology of identifying the value of θ app , and compare two different methods.

A. Finding the Critical Ca
We now describe how we determine the critical Ca computationally. We note that this methodology is valid for both the fully nonlinear and quasi-parallel systems. Initially, we shall assume a constant θ cl , i.e. we impose (11), and for simplicity we shall assume that θ cl = π/2. In the pressure-driven problem, by varying Ca and subsequently solving the steady set of equations, we can trace the state of the system by recording the horizontal distance between where the interface meets the moving wall and the symmetry wall, which we denote X rec and X adv for the receding and advancing contact lines, respectively (the same can be achieved in the force-driven problem by increasing the value of F and finding the maximum value of Ca crit ). Figure 5 shows the resulting solution curves of X rec (upper curve) and X adv (lower curve) plotted against Ca. The solid lines indicate solutions of the half liquid-plug problem, while the broken lines are solutions of the corresponding quarter RCL and ACL domains.
There are a number of important features of these solution curves. The RCL solution curve experiences a limit point (or fold bifurcation) where the curve turns around and the corresponding steady solution becomes unstable. The consequence of this is that the value of Ca where this critical point occurs marks the limiting threshold for which stable (i.e. those which can be experimentally realised) steady solutions exist, and it is, therefore, natural to associate this value with Ca crit . It is important to emphasise that the limit point occurs only for the RCL in this setup (i.e. a liquid-vacuum system). Furthermore, we note that the curves of the half liquid-plug problem completely overlap the curves for the RCL and ACL domains, the only distinction being that the quarter ACL solution curves are able to continue past Ca crit , as no unstable ACL is observed. The critical point of the full nonlinear system coincides with the fold bifurcation of the RCL, and so in order to understand the dynamics of the system before and after criticality, we only need to consider the RCL; thus, significantly reducing the computational demands.

B. Measuring the Apparent Angle
The precision with which the dynamic contact angle can be measured experimentally is limited by the resolution of the method used [22]. This is usually of the order of a few micrometres, and for optical measurements can be no better than the diffraction FIG. 6. The alternative definitions of the apparent angle. In (a) the θ app,circ is defined as the angle a fitted circle makes with the bottom plate. θ app,inf is defined as the minimum angle the interface makes with the horizontal, as measured anti-clockwise, and corresponds to the inflection point of the interface, where the curvature, κ = 0. In (b) we plot the interface angle as a function of y which demonstrates that θ app,inf can be calculated as the minimum value of θ .
limit. Thus, accurate measurement of the true contact angle is not possible, though some progress has been made [15]. A common approach is to fit a curve to the image of the interface and measure its tangent at its point of intersection with the solid surface. Alternatively, the interface can be assumed to have a quasi-equilibrium shape (e.g. a spherical cap) from which the angle may be deduced via an appropriate formula [16,22,37,46]. Neither approach has the ability to resolve significant changes in curvature very close to the contact line, such as those observed in the MD simulations (see figure 2), which occur whenever the true contact angle differs significantly from the apparent angle. Therefore, to a greater or lesser extent, the measured angle inevitably depends on the method used to measure it, and simply represents the slope of the interface at some arbitrary distance from the contact line. This is the reason why these angles are commonly described as 'apparent'.
In the MD study, [29], two methods were investigated to evaluate θ app in a systematic way that was consistent with experiment, despite the very small scale of the system. For both, multiple snapshots were averaged to account for thermal noise. Since the menisci of the liquid-plug are cylindrical at rest, the methods were based on circular fits to the liquid surface. The first approach was to estimate the slope of the interface at the point of its inflection, as shown in figure 6. This was achieved by fitting the arc of a circle to the central 50% of the meniscus (i.e. well away from the inflections) and measuring the slope of the arc at its points of intersection with planes parallel to the solid surfaces and passing through the inflections. The method appealed as being consistent with the asymptotic matching procedure used in hydrodynamic treatments of dynamic wetting [17,73].
The second approach was to mimic experiment more directly by measuring the tangents to a circular arc defined by upper and lower contact lines and passing through the apex of the meniscus at its mid-point, as shown in figure 6. This procedure is commonly used to measure the dynamic contact angle in capillary systems [22]; and because the positions of the three defining points could be measured more accurately from the simulations than the locations of the inflections, this was the method adopted. It gave advancing angles a few degrees smaller than the those found at the inflection points, but the receding angles were indistinguishable within simulation limits. The method was also used in a recent numerical study of microscopic and apparent contact angles [51].
Similarly, in the present paper we calculate θ app in two different ways consistent with those adopted in [29]. In the left panel of figure 6 the liquid-gas interface is shown by a solid line and the arc of a circle that is tangent to the interface at the line of symmetry by a dashed line. We define θ app,circ as the angle the circle makes with the horizontal as shown in figure 6. This definition, used in [29], is useful if the position of the liquid-gas interface is not well defined.
Alternatively, and again as considered in [29], we can define θ app,inf as the angle the interface makes with the horizontal at the inflection point of the curve, see Figure 6. In this definition we measure the angle along the curve using the identity and then find the minimum value that θ takes as a function of x, see right panel of Figure 6. This corresponds to where the curvature is zero and the interface has an inflection point. The value of y at the inflection point is denoted h inf , which will be commented on later in the paper. This approach of measuring θ app is more amenable to FEM calculations, because (26) can be calculated easily as the position of the interface is well-defined and was the approach used in [48], [71], [47,49,50] and [72].

V. STEADY RESULTS
In this section we will describe the steady solution space of the RCL using and comparing the CA model, where θ cl is held constant, and the VA model, where θ cl is determined by (12). First, using parameter values that are representative of the values in [29], we shall compute steady solutions in the pressure-driven VA model and present a bifurcation diagram that demonstrates that the fold bifurcation, which represents the critical speed of dewetting, still exists when using the VA model and predicts the value obtained in [29]. Then we shall vary θ 0 to investigate the effect of wettability on the value of Ca crit and compare this directly with the results of [29] using the pressure-driven and force-driven problem. Finally, we compare the predictions of the continuum model with the results previously obtained by applying the Cox-Voinov law to the data derived from the MD simulations in [29] A

. Bifurcation Diagram and General Features
We now focus our attention to the VA model and discuss steady solutions and the critical point. We emphasise that in this model we require only the static angle, θ 0 , the width of the three-phase zone, δ , and the capillary number, Ca, as specified parameters so that a steady solution can be computed. Figure 7 shows the steady solution space of the RCL domain by plotting X against Ca, as calculated numerically. The solid and dashed curves indicate, respectively, the stable and unstable solution branches of the full system and the circular markers indicate the solution branch of the QP system. The inset diagrams show streamline patterns and the interface position at (A) the critical point, (B) when θ app,circ = 0 and (C) when θ app,inf = 0. There are a number of interesting features that are worth commenting on. First, we note that even though θ cl is now dependent on Ca, the limit point still occurs. We also remark on the close agreement of the QP system and the full nonlinear system; the QP system does remarkably well in approximating the limit point for this particular value of θ 0 , although the curves diverge as X rec increases.
The solution where θ app,circ = 0 is significantly closer on the bifurcation curve to the limit point than where θ app,inf = 0. In fact, as seen from the inset interface profiles, the interface when θ app,inf = 0 (label (C)) is already significantly deformed and approaching a thin-film, whereas the profile when θ app,circ = 0 (label (B)) more closely matches the interface at the limit point (label (A)). This result is interesting, as in many works, e.g. [25], Ca crit is defined as occuring when θ app = 0. This is strictly FIG. 8. a) The measured angles as a function of Ca when λ = 0.2, θ 0 = 105.1 • . In this figure we adopt the convention of [29] where for the RCL domain Ca is negative. valid in the regime λ → 0 and we note that in our geometry the slip length and Ca take moderate values, and therefore we find that θ app = 0 at Ca crit .
The consequence of our findings is that it is not unreasonable to use the definition of θ app,circ = 0 as a lower bound on the critical capillary number when analysing outcomes in experiments and MD simulations where the existence of a smooth bifurcation curve is hidden, including smaller geometries whose dimensions are comparable to those of the slip length. Furthermore, at the limit point, we expect the dynamics of the system to be very slow, as the leading eigenvalue of the linear stability problem will be close to zero [41]. Thus in any experimental/MD setup the time frame may not be large enough to guarantee that a steady state is being approached or whether a thin-film is about to develop. Therefore, we conclude that while the position of the critical point is a well-defined threshold for the critical capillary number in calculations involving a deterministic hydrodynamic model, for experimental and MD results using the definition of Ca crit as the location where θ app,circ = 0 may be operationally acceptable. The difference between θ app,circ and θ app,inf for both the RCL and ACL is shown in figure 8 for values of λ = 0.2 and θ 0 = 105 • . As Ca crit is approached in the RCL domain the apparent angle tends to zero in a way such that θ app,circ < θ app,inf < θ cl . For the ACL domain (the shaded region in the figure), the order of the inequalities is reversed, although there is less of a distinction between θ app,circ and θ app,inf . We should also mention here that in the original MD study [29], the behaviour of the cosine of the advancing contact angle was significantly non-linear over the range of velocities investigated. As a result, the difference between θ cl and θ app,circ was significantly smaller than that depicted in figure 8, where the linear form of the MKT, (12), is used throughout. A critical test of both the CA and VA model is how well it is able to predict Ca crit when compared with the MD. For the CA model we specify λ and θ cl and then Ca crit can be calculated using the method described in [41] to find the fold bifurcation. Figure 9(a) shows the location of Ca crit as θ cl is varied for λ = 0.02 and λ = 0.2; these values roughly corresponding to the lower and upper bounds of the slip-length in the MD calculations. As can be seen from the figure, the comparison with the MD data is poor, with Ca crit underestimated. This provides motivation to implement a VA model where θ cl is a function of Ca and λ .
A much more convincing result is obtained when we apply the VA model. Figure 9(b) shows Ca crit plotted against θ cl . The curves represent the loci of the critical point as λ and, therefore, θ 0 , are varied. Note that λ and θ 0 , are expressly linked by expression (14). Here, we have chosen a range of λ that matches the MD simulations in [29]; the only parameter we have to specify is δ in (12). The solid/dotted lines are for the for the pressure-driven/force-driven problems, respectively, with δ = 0.0525, the value from the MD simulations; see appendix A. The solid markers are the QP data and the dashed line is the asymptotics described by (25).
Remarkably, the QP model replicates the full nonlinear system and the simple asymptotic formula shows excellent agreement with the QP model, despite this being in a regime when the slip-length, and indeed Ca, are not particularly small. We stress that here the asymptotic formula uses the full Ca-dependent formula for θ cl .
Evidently, the VA theory captures the same qualitative behaviour exhibited by the MD simulations and is much better at predicting Ca crit than the CA model and the force-driven problem has a slightly better quantitative fit. For the value of δ obtained from [29], the VA under predicts Ca crit for the same value of θ 0 , but if we make δ larger the comparison becomes more FIG. 9. In each chart, the markers with error bars are the MD data obtained from [29]. (a) The critical Ca as a function of the true angle, θ cl for different values of λ using the constant θ cl model given in (11). (b) The critical Ca as a function of the true angle, θ cl for the pressure-driven (solid line) and force-driven (dotted line) problems, QP system (circular markers) and the asymptotics given by (25) figure 9(c). The uncertainty in selecting the appropriate basis for the measurement of δ from the simulations is discussed in appendix A. Values of δ larger than that used here are certainly compatible with the data, depending on the criteria used to define the three-phase zone. This uncertainty is compounded by the inevitable thermal noise in the MD results, despite averaging the data over long periods relative to the simulation timescale, which means that there is an inherent difficulty in obtaining the exact steady state at the critical point in the MD simulations. This would mean that the critical point from the MD simulations should be treated as a lower bound, rather than a precise value.
As well as measuring θ cl we also measure θ app,inf and θ app,circ at the limit point. Figure 10 shows three curves that represent the Ca crit as a function of 1) θ app,circ , 2) θ app,inf and 3) θ cl . We observe that θ app,circ at Ca crit never exceeds 20 • which is consistent with experimental studies of the apparent angles at receding contact lines [13,32,54,55], where the data shows that dθ app /dCa diverges rapidly as Ca → Ca crit , and is replicated in the VA model as shown in figure 8. In these previous studies apparent angles approaching zero are reported only on surfaces that exhibit little or no contact angle hysteresis [46]. Hysteresis implies the presence of surface imperfections, such as roughness or heterogeneity. These cause local fluctuations in contact-line velocity, which may trigger film deposition prematurely. This might be the reason why [55] report that angles below about 30 degrees were inaccessible. The alternative possibility is that dewetting systems may become intrinsically unstable at some value of θ app > 0, as demonstrated here. The fact that θ app is itself an artificial construct and dependent on the method of observation adds further uncertainty to the interpretation of experimental data. Nevertheless, we comment that the results from the VA model are consistent with these experimental observations. We shall now use the Cox-Voinov law, equation (1), to help further rationalise the MD results. If we interpret θ app as θ app,circ , the approach taken by [29], we can make the approximation that θ app,circ = 0 at the critical capillary number and then (1) reduces to In figure 11, we plot θ 3 cl as a function of Ca crit , based on the solutions at Ca crit for each value of θ 0 . We can estimate the value of L * /L * m by approximating the curve as a straight line and measuring the slope. It is clear from the figure that a straight line is not wholly appropriate, but its slope gives an approximation of L * /L * m = 2.07, which compares favourably with the equivalent approximation from [29] of 2.06 (denoted L/L m in their study). This is further direct evidence that the numerical results of the model closely replicate the MD simulations.

VI. TIME-DEPENDENT RESULTS: THIN-FILM FORMATION
We now discuss time-dependent calculations and the formation/deposition of a thin liquid film. In most of the simulations that follow, we start a pressure-driven system from rest with an initially flat interface and a constant value of Ca. As shown in [41], if we choose Ca < Ca crit then the system will relax to the stable steady solution branch, as seen in the MD simulations when F * 0 < F * crit . However if we choose Ca > Ca crit a thin-film will develop, as also observed in the simulations. In this section we implement (12) and perform time-dependent calculations to understand the effect of the various parameters on the formation of this thin-film. Figure 12 is a visualisation of the velocity field using quivers to represent the strength and direction of the flow once a thin-film has developed. There are three distinct regions; a 'rim' region close to the contact line, a flat, thin-film region of height h film and a static region corresponding to the static meniscus shape. We remark that close to the contact line the flow is approximately parallel, and so a lubrication model would be an appropriate model reduction here. Far away from the contact line, the flow is certainly not parallel and, therefore, to resolve the half liquid-plug a full continuum model is required. Figure 13 shows snapshots of the evolution of the interface at different times, t, for various values of Ca (panels (a) -(d)). Panel (e) shows the time-signal of θ cl and panel (f) compares the final time-snapshot for the different values of Ca chosen. In the super-critical case (i.e. Ca > Ca crit ) the height of the thin-film is approximately constant before an almost circular cap region closes the interface. For macroscopic geometries, it is well known that the film thickness, h film scales according to the Landau-Levich-Derjaguin (LLD) law [20,45]: This value of the film height is shown as dotted lines in figure 13  increasing accuracy as Ca becomes smaller (as expected).
The contact angle at small times rapidly decreases and achieves a minimum value, before gradually increasing to a limiting value, at θ cl ≈ 59 • , as shown in panel (e); the same time-dependent behaviour was observed in the MD simulations. A key observation is that in these time-dependent calculations the limiting relative capillary number Ca is independent of Ca, as seen in figure 14(a) which is consistent with experimental and theoretical studies, for example [62]. This indicates that the flow in the film region becomes increasingly independent of the liquid-plug. In [41] it was shown that, at the RCL especially, the time-dependent trajectories of the system are similar to the steady bifurcation diagram when both are plotted in the (Ca, X) plane. The same phenomenon occurs here; see panel (b) in figure 14 where it is shown the trajectories closely match the steady bifurcation curve. This is consistent with a prediction of [14] that, for plate-withdrawal from a bath flattened in the far-field by gravity, the dynamics closely follow the unstable branch of solutions in a quasi-steady manner. In their case, where gravity plays an important role, the bifurcation curve oscillates around a fixed value of Ca, but in the pressure-driven problem this does not occur and we observe monotonic convergence towards a particular Ca. In the body-force problem we see the exact same phenomena, for both the CA and VA model, (results not shown), but leave a thorough investigation of this as a future research avenue. Finally, we can make a qualitative comparison of the numerical results to the MD results in figure 2. Figure 15 shows the equivalent half liquid plug profiles for the force-driven problem (as in the MD) obtained by computing the receding and advancing interfaces separately and combining them. As can be seen from the profiles, the qualitative comparison is strong.

VII. LARGER SCALE SYSTEMS
Having validated our model in the context of the nano-geometry using MD calculations, we can extend our analysis to investigate thin-film formation in a larger-scale geometry, for which MD simulations are prohibitively computationally expensive. We can achieve this by reducing the value of the dimensionless slip-length λ while keeping θ 0 , and therefore the physical slip-length constant, which has the effect of increasing H * , the physical channel width.
We also investigate the formation of thin-films in the limit as λ → 0 (we note that λ = 0 has no solution [38]). Using the same methods as before, we can track Ca crit as λ is varied. Figure 16(a) shows that Ca crit → 0 as λ → 0, indicating that the system becomes unstable for increasingly slower wall speeds as the scale of the system is increased. The dashed line indicates the asymptotic formula given in (25), with the full expression for θ cl used, and the dotted line shows (25) with θ cl replaced with θ 0 . Because λ ≪ δ , the difference between θ 0 and θ cl is large, see equation (12). As a result, in order to capture the numerics, the full expression for θ cl has to be included in the asymptotic formula, and then the agreement is excellent.
Panel (b) shows how θ app,circ and θ app,inf vary at the critical point. It is clear that for both measures of the apparent angle, as λ decreases, θ app → 0. This is an important observation and provides a link to the work of [24,61,62], where in the lubrication approximation they apply it is perfectly reasonable to apply the Cox-Voinov formula with θ app = 0 as a means of determining Ca crit . In a nano-geometry however, this approximation is not valid, as we have shown that θ app,crit = 0. We also find that as We now turn our attention to time-dependent results in the limit as λ → 0 with θ 0 being kept fixed. Figure 17 shows the thinfilm at t = 14.916 for values of λ = 2 × 10 −2 , 2 × 10 −3 , 2 × 10 −4 , 2 × 10 −5 (panels (a) -(d)) with Ca = 0.05, θ 0 = 64.7 • . The largest value of λ = 0.02 corresponds to the nano-channel considered in [29] while the smallest value of 2 × 10 −5 corresponds to a system 10 3 times larger, i.e. a micro-channel. Panel (e) shows the comparison of the profiles when the x normalised by X. The immediate observation is that the dimensionless film-height, h film is independent of λ once the thin-film has had sufficient time to develop, so that the physical film height will scale linearly with the system size. This is especially evident when comparing the interface profiles for λ = 2 × 10 −3 , 2 × 10 −4 , 2 × 10 −5 . Thus, sufficiently far away from the contact line the structure of the thin-film is independent of the size of the geometry (in physical systems this is measured relative to the physical width of the channel). The 'rim' region is however highly dependent on λ and the details of the contact line angle (see [30]); the smaller the system (larger λ ) the larger the 'rim' near the contact line. Therefore, the physical rim height will increase slower than linearly as system size is increased.

VIII. CONCLUSION
We have developed a novel molecularly-augmented continuum model, based on a variable true contact angle, that describes the dynamics of a liquid bridge between two parallel plates, and, more generally, describes the RCL and ACL physics. By solving the resulting set of equations numerically, we are able to interpret the maximum speed of dewetting as a fold bifurcation in the steady bifurcation diagram. We find that the maximum speed of wetting Ca crit , calculated as a function of θ cl , is qualitatively similar to the MD simulations described in [29] and that the estimate of L * /L * m is in excellent agreement. As well as showing good agreement with the MD simulation, the advantages of this approach is that by replacing the assumption that θ cl is constant with the constraints the issue of deciding what θ cl should be in any hydrodynamic calculation is removed, as it is naturally determined, through (29), as part of the solution. Furthermore, whereas in previous approaches the slip-length and θ cl had to be specified as control parameters, in this model the only hydrodynamic parameter we have to specify is the slip length. With this parameter being difficult to measure, invariably it has been used as an additional fitting parameter that can cover-up for inaccuraces in the constant angle model. We do however have to estimate δ , the width of the TPZ from MD simulations and this provides an additional parameter that has to be known in advance, although this parameter can be far more accurately specified than slip lengths. The comparison between the MD simulations and the VA model is strong, and although some of the physics present in the MD calculations are absent, for example the disjoining pressure, we conclude that the VA model contains the minimum ingredients required to replicate the physics contained in the MD calculations, at least before the thin-film ruptures (see, for example [43,76]). Our results also illuminate the values of θ cl and θ app when a partially wetted substrate is withdawn from a pool of liquid at capillary numbers greater than Ca crit . Experiments have shown that attempts at forced dewetting cause the (three-dimensional) contact line to slant at an angle relative to the direction of withdrawal, such that the capillary number in the direction normal to the contact line remains constant at Ca crit . These observations of avoided critical behaviour led to the postulate of a maximum speed of dewetting [11]. Presumably, θ cl and θ app along the slanted contact line are the smallest possible consistent with a stable flow without film deposition, i.e. those associated with the turning point in the steady phase diagram. For θ cl this is θ cl,crit . For θ app it depends on how the angle is measured.
We are easily able to extend the VA model to larger systems, which are prohibitively computationally expensive for MD calculations, and by examining the thin-film formation in these systems when Ca > Ca crit , we are able to demonstrate that the relative height of the thin-film is independent of size of the system and weakly dependent on Ca. Differences in the interface profile occur close to the contact line, as indicated by the size of the 'rim' that develops, but sufficiently far away from the contact line the relative heights of the thin-film are nearly identical.
Another advantage of the framework is practical, in that the computational time for these calculations is O(minutes) using the open-source oomph-lib framework with state of the art linear algebra solvers, rather than O(days) for the MD simulations. As we are able to obtain the velocity and pressure fields in addition, this unified model has excellent potential for researchers wishing to combine the best aspects of the hydrodynamic and molecular theories in their work. We also remark that viscous and inertial effects can be incorporated in this model by, for example, treating the gas-phase using a lubrication approximation; see [41]. We also remark that the QP model and associated asymptotic results, while not resolving the flow-field, are useful for validation, as demonstrated here.
Nevertheless, there remains a need for more physical experiments with emphasis on the RCL up to the point of film deposition, since, as we have seen, this encodes much valuable information concerning the contact angle on the microscopic scale. While there is a very large body of literature on film deposition, such as that which occurs when a solid surface is withdrawn from a pool of liquid, and much published data on advancing contact angles, comprehensive measurements of dynamic receding angles on partially-wetted surfaces are, unfortunately, rare. A resurgence of interest is overdue.
for net displacement of the contact line at velocity U * cl , work must be done to favour events in the desired direction. This work is provided by the out-of-balance surface tension force that arises when the equilibrium is disturbed: f * = γ * L (cos(θ 0 ) − cos(θ cl )), where γ * L is the surface tension of the liquid. According to the model, as the TPZ moves across the solid surface, this work is expended at n * interaction sites per unit area swept. Application of the Frenkel-Eyring theory of stress-modified activated rate processes [31,34] then leads to the principal equation linking U * cl and θ cl : where k B and T are, respectively, the Boltzmann constant and the absolute temperature. Since its inception, this equation has proved very effective in correlating experimental and MD data for a wide range of systems. For examples see [7,8,23,57]. In the interpretation of experimental data, the interactions sites are usually assumed to be uniformly distributed, so that λ * 0 ≈ 1/ √ n * ; thus, reducing the unknowns to just two: κ * 0 and λ * 0 . For small arguments of sinh, typically when θ cl is not too far from θ 0 (true for θ cl at the RCL in the MD data investigated here) or γ * L is small, this reduces to a linear relationship: which may be written as where ζ * is the coefficient of contact-line friction (per unit length of the contact line): This single coefficient quantifies the localised resistance to the displacement of the contact line.
In previous MD studies [6,12] it has been shown that both contact-line friction and slip between a liquid and a solid depend on the same thermally-activated molecular events. Whereas, at the contact line, the principle driving force comes from the outof-balance surface tension acting across the TPZ, for the latter it is provided by the viscous shear stress acting across the whole solid-liquid interface: µ * L (∂ u * /∂ z * ) = β * λ * , where β * is the slip coefficient and λ * the Navier slip length (i.e. the distance into the solid at which the extrapolated fluid velocity vanishes). Because of the common mechanism, it follows that the two coefficients are directly related; specifically, β * = ζ * /δ * ; (A5) hence, λ * = δ * µ * L /ζ * .
This relationship has been validated by molecular-dynamics simulations, in which both the contact-line friction and the slip length have been measured for the same system over a range of equilibrium contact angles [12,28]. Good agreement has been shown for both Lennard-Jones liquids and atomistically simulated water on molecularly smooth carbon-like surfaces. That said, a precise correlation hinges on the value of δ * . For the Lennard-Jones liquids, the value selected was assessed from the velocity profiles across the TPZ. For the simulated water system, the distance over which the density of the of the liquid in contact with the solid fell to zero was used. See figure 10 in [12] to compare the two approaches. Arguments may be made for both. For the Lennard-Jones system, the difference in the result was in the region of 30%. Slip lengths were smaller if the density profile was used. In addition, the value of δ * appeared to depend weakly on both contact-line velocity and the equilibrium contact angle. Based on the existing data, while (A6) appears to be physically justified, a precise understanding of the subtle influences in play requires more work. The value of δ * found for the coarse-grained water simulations [29] was 0.93 ± 0.14 nm based on the density argument. We use this value in the present paper. If, (A6) is accepted, at least in principle, it allows us to rewrite (A3) in dimensionless variables, as Ca cl = λ δ (cos(θ 0 ) − cos(θ cl )) .
In the stationary frame of the liquid-plug between two solid walls moving at velocity U * wall , this becomes Ca U wall − ∂ x cl ∂t = λ δ (cos(θ 0 ) − cos(θ cl )) , which is (12) in the main body of the paper when we set the nondimensional wall speed to be U wall = −1. Furthermore, and perhaps more significantly, we know that contact-line friction depends strongly on the equilibrium contact angle. This means that the same is true for the slip length. As has been shown [6,7,9,23] the frequency κ * 0 is related to the equilibrium contact angle by where v * L is the molecular flow volume in the Frenkel-Eyring theory. This leads to ζ * ∼ (n * µ * L v * L /λ * 0 ) exp [−γ * L (1 + cos(θ 0 )) /n * k B T ] and, hence, to λ * ∼ δ * (λ * 0 /n * v * L ) exp [−γ * L (1 + cos(θ 0 )) /n * k B T ] .
This suggests the general (dimensionless) form In the present paper we have used this expression, (14), to fit the slip length calculated from the MD data in table 1 of [29].