A mathematical model and mesh-free numerical method for contact-line motion in lubrication theory

Abstract We introduce a mathematical model with a mesh-free numerical method to describe contact-line motion in lubrication theory. We show how the model resolves the singularity at the contact line, and generates smooth profiles for an evolving, spreading droplet. The model describes well the physics of droplet spreading–including Tanner’s Law for the evolution of the contact line. The model can be configured to describe complete wetting or partial wetting, and we explore both cases numerically. In the case of partial wetting, the model also admits analytical solutions for the droplet profile, which we present here. Article highlights We formulate a mathematical model to regularize the contact-line singularity for droplet spreading. The model can be solved using a fast, accurate mesh-free numerical method. Numerical simulations confirm that the model describes the quantitative aspects of droplet spreading well.


Introduction
When a droplet of fluid (surrounded by a gaseous atmosphere) is deposited on a substrate, it spreads until it reaches an equilibrium configuration. At equilibrium, the angle between the liquid-gas interface and the solid surface is measured (conventionally, through the liquid), to yield the equilibrium contact angle eq . If the angle is less than 90 • , the substrate is deemed hydrophilic, whereas if the angle exceeds 90 • , the substrate is deemed hydrophobic [1]. Droplet spreading then describes the dynamic phase before the attainment of this equilibrium.
In droplet spreading, the point of contact between the substrate and the gas-liquid interface is in motion. And yet this contradicts the classical no-slip assumption in viscous fluid flow, which stipulates that there should be no relative motion between a substrate in contact with a fluid [2,3]. The resolution of this paradox is that there is missing physics, and that on a sufficiently small scale, there is slip, the dynamics of which are governed by the interactions between the fluid molecules and the substrate molecules [4]. These molecular-level interactions can be incorporated into a macroscopic fluid model via a so-called regularization technique. Broadly, there are three regularization techniques in the literature -the slip length [5], the precursor film [1], and the diffuse interface [6,7]. Although methodologically distinct, these yield the same qualitative and quantitative results when used to model droplet dynamics. This consistency between the different approaches gives a solid justification for the general approach of model regularization.
The purpose of this article is to introduce a novel model regularization, albeit one in the spirit of those just described. The focus of the work is on mode regularization for the case of thin-film flows (also called lubrication flows). For simplicity, we focus on two-dimensional configurations (or equivalently, three-dimensional axi-symmetric configurations), however, the generalization to three dimensions is straightforward.
Thin-film flows arise in the context of hydrophilic substrates, where the equilibrium shape of the droplet is such that the typical size of the droplet base R greatly exceeds the maximum droplet height h 0 , leading to a small parameter = h 0 ∕R , with ≪ 1 . In this context, the Navier-Stokes equations reduce down to a single equation for the interface height (the so-called Thin-Film Equation). In this context also, the interface height is conventionally written as z = h(x, t) . Thus, the coordinate x and is in the plane of the substrate, and the z-coordinate is orthogonal to the substrate (e.g. Fig. 1).

Motivation for this work
In the context of droplet spreading, the Thin-Film Equation inherits the contact-line singularity from the full Navier-Stokes equations; the singularity can be regularized by introducing a slip length [5,8], or a precursor film [1,9]. Both these regularizations give accurate and consistent descriptions of contact-line spreading [10], albeit with some drawbacks-the classical Navier slip-length model has a logarithmic stress singularity at the contact line. We note, however, that the stress singularity can be alleviated using the approach in Reference [11]. While as the precursor-film model requires a precursor film to be present that extends indefinitely beyond the droplet core. Although physically, such a precursor film does exist, it has a very small thickness (10-100 nm Fig. 1 Schematic description of droplet spreading on a substrate 1 3 [9]), meaning that such a small scale must be resolved in the model: in particular, the numerical grid size must be at least as small as the precursor-film thickness [12]. Furthermore, the resulting equations are quite stiff numerically [13]. Although this approach is just about feasible for millimetre-scale droplets, it may not be feasible for larger ones. Beyond the millimetre scale, an unphysically large precursor-film thickness can be used in numerical investigations (and the results checked for robustness to changes in the value of the precursor-film thickness). The Diffuse Interface Model has been proposed as a more general regularization of the contact-line singularity problem, valid for the full Navier-Stokes equations beyond the lubrication limit [14]. The Diffuse Interface Model has been implemented for droplet spreading in the thin-film lubrication approximation [6].
Our contribution in this article is to introduce a novel regularization technique similar in spirit to the Diffuse Interface Method-we formulate a theory of droplet spreading involving a smooth interface height h , as well as a sharp interface height h, which interact via a convolution operator and an evolution equation. In a previous work [15], the idea behind this regularization was introduced in the context of thin-film lubrication flows-the so-called Geometric Thin Film Equation.
In the present article, we extend this previous work by introducing a particle method as a novel and highly accurate solution method for the Geometric Thin Film Equation. Also, Reference [15] was for complete wetting-in this work we extend the Geometric Thin-Film Equation to describe partial wetting as well.
The Geometric Thin-Film Equation can be viewed as special case of a mechanical model for energy-dissipation on general configuration spaces-the derivation of the general model involves methods from Geometric Mechanics such as Lie Derivatives and Momentum Maps [16]-hence the name Geometric Thin-Film Equation. The main advantage of this new method so far has been the non-stiff nature of the differential equations in the model, which leads to robust numerical simulation results. A second advantage (the main focus of the present work) is that the Geometric Thin-Film equation admits so-called particle solutions. These give rise to an efficient and accurate numerical method (the particle method) for solving the model equations.
The basis for the particle method lies in the structure of the Geometric Thin-Film Equation: the model includes a 'smoothened' free-surface height h(x, t) and a 'sharp' free-surface height h(x, t), related by convolution, h(x, t) = * h(x, t) , where ≥ 0 is a filter function with a characteristic lengthscale . The model admits a delta-function solution for the sharp free-surface height h( where is the Dirac delta function, and x i (t) is the (time-dependent) centre of the Delta function. The deltafunction centres {x i (t)} N i=1 satisfy a set of ordinary differential equations: where V i is a velocity function which can be derived from the Geometric Thin-Film Equation. Thus, the screened free-surface height admits a regular solution In this context, the weights w i and the delta-function centres x i (with i ∈ {1, 2, … , N} ) can be viewed as pseudo-particles, which satisfy the first-order dynamics (1). We will demonstrate in this article another advantage of this particlemethod: it is a mesh-free method that automatically accumulates particles in regions where | xx h| is large, thereby mimicking the effects of adaptive mesh refinement, with none of the computational overheads associated with that method. A final key advantage of the particle , with w i ≥ 0 and ≥ 0 , the numerical values of h(x, t) are guaranteed never to be negative, hence the numerical method is manifestly positivity-preserving.

This work in the context of environmental fluid mechanics
Before introducing the method, we place the method in the context of Environmental Fluid Mechanics, focusing in particular on droplet spreading on plant leaves. In general, surfaces can be further classified as being (i) super-hydrophobic ( eq > 150 • ), hydrophobic ( 90 • < eq < 150 • ), hydrophilic ( 10 • < eq < 90 • ), or super-hydrophilic ( eq < 10 • ). Plant leaves display this wide range of contact angles. Superhydrophobic surfaces have been frequently found in wetland plants, where the superhydrophobic surface prevents a buildup of water on the leaves, which could otherwise promote the growth of harmful micro-organisms and limit the gas exchange necessary for photosynthesis [17].
A particularly well-studied plant with superhydrophobic a leaf surface is the Lotus plant (Nelumbo nucifera), with a contact angle of about 160 • [17,18]. The leaves of the Lotus plant also demonstrate a low contact-angle hysteresis, such that water droplets roll off the leaf surface when at a low tilt angle ( 4 • ). During rolling, contaminating particles are picked up by the water droplets, and are then removed with the droplets as they roll off [17]. The leaf structure of the Lotus plant has been studied using Scanning Electron Microscopy, and reveals a hierarchical microstructure made up of micron-scale pillars (cell papillae) and a randomly covered by a smaller branch-like nanostructure [19] (the wax crystal, with microstructures ∼ 100 nm in diameter) [20]. Such biological microsturctures are the inspiration for engineered super-hydrophobic surfaces [20].
In contrast, plants with superhydrophilic leaf surfaces are often found in tropical regions [21]. A water droplet that spreads on a superhydrophilic surface will spread and form a wide, flat droplet-effectively a thin film. Evaporation in such films is more efficient than in a spherical droplet, due to the increase of the water-air interface. Thus, water evaporates from a superhydrophilic leaf much faster than that from a hydrophilic or superhydrophobic one, thereby keeping the leaf dry, reducing the accumulation of harmful micro-organisms on the leaf surface, and increasing the gas exchange with the environment, for the purpose of photosynthesis [17]. A well-studied hydrophilic plant is Ruellia devosiana, wherein the superhydrophilic property of the leaf surface is due both to the leaf microstructure, and to a secretion of surfactants by the leaf, which both promote spreading [22].
Other plants which exploit hydrophilicity are the carnivorous plants of the Nepenthes genus (e.g. Fig. 2), the perisotone of which is a fully wettable, water-lubricated anisotropic surface [23]. Insects landing on the peristone or rim of the plant effectively 'aquaplane' down the plant rim [23] before being captured by the viscoelstic fluid inside the pitcher [24].
We emphasize here that both the super-hydrophobic and super-hydrophilic plant surfaces are extreme cases. In a comprehensive investigation of 396 plant species out of 85 families growing in three different continents at various elevations (including 792 leaf surfaces total), only forty leaf surfaces (5.1% out of 792) were found in these extremes, including 24 super-hydrophobic surfaces and 16 super-hydrophilic surfaces [21]. Thus, most leaf surfaces lie inside these extremes. The droplet-model introduced in this paper is most relevant to hydrophilic cases where the equilibrium contact angle is small. Hence, droplet spreading over the leaves of Ruellia devosiana should be well described by our model, but not droplet spreading over lotus leaves. However, we anticipate that future applications of the regularization technique introduced herein will be able to handle more such cases.
Finally, we note that several researchers have further investigated droplet spreading using lubrication theory for non-isothermal configurations, which have important applications in cooling-both environmental and industrial. Readers are referred to standard references [25,26]; furthermore, it can be anticipated that the methods developed herein can be extended to such complex systems as part of future work.

Plan of the paper
In Sect. 2 we present the classical Thin-Film Equation as a model of droplet spreading. Certainly, this is a very well-established topic, however, we briefly summarize this topic here because it enables us to clearly mark out the point of departure of the present work. Thereafter, in Sect. 3 we introduce the Geometric Thin-Film Equation as a regularized model which enables contact-line motion. In Sect. 4 we introduce the particle method for generating numerical solutions of the Geometric Thin-Film Equation. In Sect. 5 we present results for complete wetting. In Sect. 6 we demonstrate how the Geometric Thin-Film Equation can be extended to the case of partial wetting, and we present numerical results for that case also. Concluding remarks are given in Sect. 7.

Review of classical theory
In this section we review the classical theory of the Thin Film Equation, including the problem of the contact-line singularity. The purpose of this review is to put the Geometric Thin-Film Equation into the context of the classical theory; the theoretical formulation of the Geometric Thin-Film Equation is therefore presented subsequently in Sect. 3.

Classical thin-film equation
We review the derivation of the classical Thin Film Equation for a flow in two dimensions, with spatial coordinates x and z (we refer the reader Reference [27] for the details). The starting-point is the kinematic condition valid on the free surface z = h(x, t): The fluid flow is assumed to be incompressible, such that u x + v y = 0 . The incompressiblity condition can be integrated once to give Eqs. (2)-(3) can be combined to give In the lubrication limit, the velocity u(x, z, t) satisfies the equations of Stokes flow, hence where P is the fluid pressure and is the constant dynamic viscosity. We integrate the first equation of the pair in (5) once with respect to z to obtain The standard interfacial condition is that the viscous stress u∕ z should vanish on the free surface z = h(x, t) . Thus, Eq. (6) becomes u∕ z = −1 ( p∕ x)(z − h) . Applying the noslip boundary condition the u-velocity profile becomes

3
The pressure P is identified with the Laplace pressure, P = − xx h , where is the surface tension and h xx is the interfacial curvature in the longwave limit. Hence, Eq. (9) becomes: Substituting Eq. (10) into Eq. (4) gives:

Contact-line singularity
In the context of droplet spreading, it is desirable to propose a similarity solution to Eq. (11), corresponding to a self-similar droplet that retains some overall structural properties even as the base of the droplet spreads out. Dimensional analysis indicates that the similarity solution should be: where h 0 and R are as given in Fig. 1 and t 0 is a timescale to be determined. Substitution of Eq. (12) into Eq. (11) yields: The timescale t 0 is chosen to be the capillary timescale, such that Thus, Eq. (13) becomes: The appropriate droplet-spreading boundary conditions for Eq. (14) are where 0 corresponds to the outermost extent of the droplet. Thus, the position a at which the (microscopic) contact line touches down to zero is described by a(t)∕R = 0 (t∕t 0 ) 1∕7 . Unfortunately, the similarity solution (12) with the boundary conditions (15) fails to exist; instead, f ( ) degenerates into a Dirac delta function centred at = 0 , and the droplet does not spread [28]. The reason for this failure is that the no-slip condition (7) is inconsistent with the phenomenon of droplet spreading. When the model (11) is applied to droplet spreading, the 1 physics which permits slip to occur on sufficiently small scales is missing. The missing physics is then put into the model as part of a regularization. For instance, by allowing for slip on a sufficiently small scale s , Eq. (11) becomes: Equation (16) is the Thin-Film Equation with a slip-length model. Through a procedure of matched asymptotic expansion, the leading-order behaviour of the dynamic contact angle as shown in Reference [29] is given by where Ca = ( ∕ )(dx cl ∕dt) is the capillary number based on the contact-line velocity; here x cl is the position of the contact line, which depends on t. Here also, b is a constant. From this description, it can be inferred that Hence the leading-order behaviour of Eq. (17) is given by x cl ∼ t 1∕7 , which agrees with the so-called Tanner's Law. Similar calculations were done in Reference [30] for radially symmetric droplets in three spatial dimensions, from which the spreading rate x cl ∼ t 1∕10 was found.
The slip-length model therefore provides for a resolution of the contact-line singularity. However, the stress h xx remains singular at the contact line. For these reasons, an alternative regularization of the Thin-Film Equation has been proposed, namely the Precursor-Film model [1,9]. Following Reference [15], in this work we present the Geometric Thin-Film Equation as an alternative regularization, the advantage of this approach as we reveal in subsequent sections is the remarkable simplicity of the numerical solutions produced by this model.

Geometric thin-film equation: theoretical formulation
In the framework of the Geometric Thin-Film Equation, the starting-point is the assumption that there is missing physics on a small scale. Instead of modelling the missing physics, it is parametrized. As such, h(x, t) is used to denote the interface location in a crude model with missing physics -which we call here the 'noisy' interface location. The noisy interface location is to be smoothened by a filtering operation, to produce a smoothened, more accurate, estimate of the interface location, which we denote by h(x, t) . The noisy interface location may be different from the true interface location -for instance, the noisy interface location may be zero, whereas the true interface location may be close to, but different from zero -as would be the case if the noisy interface location was obtained through an incomplete model with missing small-scale physics.

Model A
To take account of the fact that h represents a smoother description of the interface location than h, we propose that h and h be connected via the expression where is a fluctuating quantity with mean zero and variance 2 . Then, h(x, t) can be made into an accurate estimate of the interface location by minimizing the total interfacial energy subject to a fixed-variance-constraint: Here, S is the system size (we take S → ∞ in what follows), and is a positive constant representing the surface tension. The constraint (21) enforces a fixed level of uncertainty between the model with missing physics and the smoothened model. In practice, we minimize the surface area in the long-wave limit: in this limit | x h| is small, and the square root in Eq. (20) can be well approximated by 1 + (1∕2)| x h| 2 . Thus, the energy in Eq. (20) becomes As only energy differences are important, the constant term here can be ignored. Furthermore, as the droplet profile rapidly decreases far from the droplet core at x = 0 , we can take S → ∞ in the foregoing analysis, this then gives the required functional form for the energy in the longwave limit: Equation (22) with constraint (21) (with S → ∞ ) is a constrained minimization problemto solve it, one would introduce an energy functional with a Lagrange multiplier: One would then compute -The evolution of h and h must be such that L tends to a minimum over time; dx must be conserved quantities, reflecting underlying principles of conservation of fluid mass.
Under these conditions, the evolution equation for h must be of a generic conservative-gradient-descent type, hence: where M ≥ 0 is a mobility function to be determined. The evolution equation for h may be similar. However, for simplicity, we may assume that h relaxes instantaneously to a smoothened form of h, hence L∕ h = 0 , hence or Equation (28) establishes a natural smoothing operation and hence, smoothing kernel for the formulation, namely, the Helmholtz kernel.
We now use L∕ h = (h − h) . We compare this to Eq. (27), and get L∕ h = − xx h . Substitution of these results into Eq. (26) yields: The physical model for h is completed by specifying the mobility. This is done by reference to the classical theory in Sect. 2. However, instead of M = (1∕3 )h 2 we take the reason for using h 2 in the mobility becomes apparent when we look at particle-like solutions of the regularized model (Sect. 4). Finally, the value of the Lagrange multiplier is chosen at each point in time to reflect the model uncertainty: We refer to this model with a fixed level of uncertainty as Model A.

Model B
In practice, recomputing the Lagrange multiplier at each time t is a difficult task numerically. However, an equivalent model can be formulated by introducing an unconstrained functional, Here, is a fixed parameter of the model. The dynamical equation is the same as before (Eq. (29)), as is the mobility; however, now h is computed as We refer to this model with uncertainty on a small scale as Model B. Here, we have introduced the standard notation for smoothing kernels [31]: This can in turn be written in several further ways: 1. The inner-product pairing of x h with x h :

Higher-order smoothing
For the purpose of generating particle-like solutions of the regularized model (e.g. Sect. 4), smoothing with the Helmholtz kernel is not sufficient. Therefore, in this paper, we work with a higher-order smoothing. We take as before (specifically, Eq. (36)), with evolution Eq. (39) and smoothing kernel h = K * K * h -this is a straightforward extension of the basic model. We therefore summarize in one place the model studied in this work: The aim of the remainder of the paper is to explore numerically the solutions of Eq. (40)we will use = K * K to denote the double Helmholtz kernel; specifically, Hence, h = * h . The need for such higher-order smoothing will become apparent in Sect. 4.

3
We emphasize that the present higher-order smoothing is introduced to give the required degree of regularity to the numerical solutions (e.g. Sect. 4). Therefore, the higher-order smoothing is introduced here for mathematical convenience. However, the higher-order smoothing can be given a physical basis through the construction of appropriate energy functionals, as in Sect. 3.2, although we do not pursue that approach here.

Discussion
Summarizing our work so far, we have introduced a regularized thin-film equation where the missing small-scale physics is not modelled, but is instead parametrized. The idea of the model is that h(x, t) provides an incomplete description of the droplet evolution (but which nonetheless contains important physical information, such as the problem dependence on the viscosity and time t. A refined description of the interface is then obtained via a smoothened interface profile h = * h . Overall, the model evolves so as to minimize the interfacial energy (area) while keeping the difference between h and h as small as possible.
The model as formulated envisages that h(x, t) is an incomplete description of the interface profile (with missing small-scale physics). The missing physics is encoded either as a fixed level of uncertainty between h and h (model A), or such that the uncertainty in the model description occurs below a lengthscale (model B). These two models are equivalent, although model B is preferred for computational simplicity. The model is concerned with physics on the droplet scale, which is typically on the millimetre or micrometre scale. The model does not include explicitly the intermolecular forces between the droplet and the substrate, which occur typically on the nanometre scale [1]. Therefore, the uncertainty in the model can be said to occur on the nanometre scale, with the value of chosen accordingly. However, our calculations (e.g. Appendix C) show that a much larger value of can be chosen without affecting the large-scale droplet spreading.
The  (42) is then a very simple instance of a general theory of gradient-flow dynamics [16] which uses geometric mechanics (Lie Derivatives, Momentum Maps) to formulate an energy-dissipation mechanical model for general configuration spaces -hence the Geometric Thin-Film equation. Equation (40)-(42) can furthermore be identified as a special case of Darcy's Law [16], where the generalized force is f = − x ( ∕ h) , the Darcy velocity is U = Mf , and the fluxconservative evolution is equation for the conserved scalar quantity h is h t + x (hU) = 0 . These insights will be used in formulating the Geometric Thin-Film equation for the case of partial wetting in Sect. 6, below.
We remark finally here on the intriguing connection between Model A and 'denoising' in Image Processing -in Image Processing one is given a noisy image h and it is desired to produce a smoother image h while keeping the difference between the noisy and the smooth image at a fixed level . This is achieved by minimizing a functional such as Eq. (23) [33,34]. Our evolution Eq. (29) for the free-surface height h(x, t) in a thinfilm flow is equivalent to carrying out denoising on a (one-dimensional) image in Image Processing.

Geometric thin-film equation: particle solutions
The general theory of geometric dissipative mechanics introduced in Reference [16] includes many examples where discrete, particle-like solutions (such as Eq. (1)) are admitted. Motivated by these examples, we seek simplified solutions of Eq. (40) of the following form: where N is a positive integer corresponding to a truncation of an infinite sum, w i ≥ 0 are weights to be computed, and (⋅) is the Dirac delta function. The motivation for seeking out such highly simplified particle solutions is that they make the task of solving the partial differential Eq. (40) numerically very simple: instead of discretizing a fourth-order parabolic-type partial differential equation and solving it numerically, we can solve a set of ordinary differential equations for the delta-function centres x i (t) using standard timemarching algorithms. This simplifies the numerical computations greatly, and is similar to the concept in the Method of Lines. We refer to the weights w i together with the deltafunction centres x i (t) as the 'particles' -thus, we are concerned with a particle-solution of Eq. (40). We describe the construction of such particle-solutions in what follows. The weights w i ≥ 0 are chosen such that where (x) is an arbitrary smooth, integrable test function. This limit can be satisfied by taking where L is a lengthscale such that the set of all points where h 0 (x) is non-zero is contained in the interval [−L, L] ; we must also take Then, Eq. (44) is satisfied automatically. We now multiply both sides of Eq. (40) by the test function (x) , integrate from x = −∞ to x = ∞ , and apply vanishing boundary conditions at these limits. We thereby obtain where ⟨⋅, ⋅⟩ denotes the standard pairing of square-integrable functions: We substitute Eq. (43)  give a so-called singular solution to the Geometric Thin-Film Eq. (40). The centres of the Dirac delta functions x i (t) with associated weights w i are identified as pseudo-particles, and the velocity V i of the i th pseudo-particle is identified with the right-hand side in Eq. (53),

Key properties of the particle evolution equations
We notice that in Eq. (55), evaluation of ′′′ is required -this is the rationale for our choice of the double Helmholtz kernel as the smoothing kernel in Eq. (40). Using the single Helmholtz kernel K would not be sufficient, as K ′′′ is singular at the origin. Furthermore, as the reconstructed interface profile h( ) involves the positive weights w i and a positive kernel ≥ 0 , the particle method is positivity-preserving: if h and h are initially positive, then the stay positive for all time. The numerical particle method is manifestly positivity-preserving, this is a key advantage as a numerical method that led to erroneous negative values of h and h would produce unphysical results.

Numerical solutions using the particle method
In this paper, we solve the Geometric Thin-Film Eq. (40) numerically using the particle method (55). To demonstrate the accuracy of the novel particle method, we compare the particle method to a standard fully-implicit finite-difference method. It can be noted from the particle method (specifically Eq. (55)) that N ordinary differential equations are to be solved; each of the N right-hand-side (RHS) terms requires a summation over all other particles (cf. Eq. (43b) and the second entry in Eq. (55))-this suggests the particle method has computational complexity O(N 2 ) . However, the number of floating-point operations to be performed in evaluating the different RHS terms can be dramatically reduced by symmetry operations, to give an overall computational complexity O(N)-this is the so-called fast-particle method. We give details of the fast particle method and our fully-implicit finite-difference method in Appendix A.
The transformations described in Appendix A rely on the properties of the one-dimensional kernel function (x) . The analogous kernel function in higher dimensions does not possess the same properties, rendering the one-dimensional fast-particle method a somewhat special case. This issue has already been encountered for fast-particle methods for the Camassa-Holm equation [35]-there, the development of a fast-particle method in higher dimensions can be achieved via a fast multipole expansion. The same approach may apply to the Geometric Thin-Film Equation.
Throughout the paper, we use the value = 0.05 in the numerical computations. This is justified in Appendix C, where show that the results of the simulations are insensitive to the choice of the parameter -provided is sufficiently small. This is advantageous, since the numerical solver must resolve the finest scale in the problem -which is precisely . Finally, both the particle method and the fully-implicit finite-difference method are solved in MATLAB. The particle method makes use of MATLAB's built-in time-evolution algorithms for ordinary differential equations, notably, ODE45 and ODE15s.

Complete wetting
In this section we present numerical results for complete wetting for the Geometric Thin-Film Equation. Equation (40) clearly corresponds to the case of complete wetting: the energy functional E = (1∕2) ∫ ∞ −∞ | x h| 2 dx is penalized, meaning that the system evolves to minimize the curved part of the droplet interface, that is, the part of the droplet interface in contact with the surrounding atmosphere.

Non-dimensionalization and initial conditions
To characterize this spreading phenomenon, we solve the Geometric Thin-Film Eq. (40) in dimensionless variables. The interfacial height is made dimensionless on a lengthscale h 0 , which is proportional to the initial maximum droplet height h max . In this section, we take the rather non-standard value h 0 = (8∕3)h max for h 0 , this is done here to compare with Reference [15]. 1 The x-coordinate is then made dimensionless on the initial droplet base L. Finally, time is made dimensionless on the capillary timescale = (3 L∕ )(L∕h 0 ) 3 . The ratio = h 0 ∕L must be small, for inertial effects to be negligible, and hence, for the lubrication theory underlying Eq. (40) to be valid. Henceforth, we assume that all of the relevant variables have been made dimensionless in this way. The initial condition for the droplet height therefore reads: Thus, the droplet area (which is the analogue of droplet volume in two dimensions) is therefore fixed as ∫ h(x, t = 0) dx = 1∕4 . Both the droplet area ∫ h(x, t) dx and ∫ h(x, t) dx are conserved under the evolution Eq. (40). The simulations are carried out in a finite spatial domain x ∈ [−2, 2] with periodic boundary conditions -this condition also establishes the limits of integration on the foregoing integrals.

Results
We solve Eq. (40) with the initial condition (56). The numerical calculations indicate that the particle method and the finite-difference method produce results that are indistinguishable by eye: we therefore show only results for the particle method. A yet more rigorous comparison between the two methods is also presented herein-this analysis also justifies the number of particles used in the calculations as N = 800 , this can be deemed equivalent to a grid spacing of x = 2L∕N = 0.005 in the finite-difference method. A first set of results is shown in Fig. 3. In Fig. 3a we show a space-time plot of the smoothened free-surface height h where the spatial grid is evaluated at the particle positions x i (t)-effectively, a discretization of h on a non-uniform grid. From this plot, the region in space where where h is significantly different from zero increases over time, demonstrating that the droplet is indeed spreading. Fig. 3b shows a snapshot the filtered surface height and its slope at t = 50.
The particle locations are shown explicitly in Fig. 3b-there is a high concentration of particles in the regions of high curvature-this is discussed in more detail in what follows.
In Fig. 4, we plot the particle trajectories x i (t) for the solution of the Geometric Thin-Film Equation via the particle method.
Intriguingly, the particles are seen to accumulate at regions where | xx h| is large, giving a higher spatial resolution in regions of high interfacial curvature. When a finite-difference or finite-volume solver is able to execute local grid refinement in regions of where the spatial derivatives are large in magnitude, this is because of adaptive mesh-refinement. Here, the particle method demonstrates a tendency to mimic the effect of adaptive mesh-refinement, without having to implement adaptive-mesh procedures explicitly. Furthermore, in Fig. 4 we have used the built-in MATLAB ODE solvers to generate the space-time plot, meaning that the adaptive mesh refinement is performed in the temporal as well as the spatial domain.
A further advantage of the particle method is that it provides a numerical description of the contact line by simply following the trajectory of a particle starting near the contact line, say x k (t) such that x k (0) = 0.5 . This is also shown in Fig. 4, where the contact line is found to satisfy Tanner's Law, with x k (t) ∼ t 1∕7 . The scaling persists even at late time, indicating that finite-size effects are not important here. Indeed, the particle method is intrinsically free from finite-size effects: since the kernel function used extends to infinity in both directions, the particle method and the accompanying ODE solver seamlessly track Fig. 4 Evolution of the particle trajectories x i (t) (logarithmic scales on both axes). The colors indicate the weight corresponding to each particle w i . The line t 1∕7 is imposed to show that the trajectories follow a power law at late time the particles as they move away from the origin and hence, as the droplet spreads. This is a key benefit of using the particle method.
In contrast, in the finite-difference method, the position of the contact line is obtained by extrapolation: the value x * (t) = argmax x [−h(x, t)] is computed, the tangent line at the point (x * (t), h(x * , t)) is computed, and the nominal position x ref of the contact line is taken to be the intersection of the tangent line with the x-axis. This procedure was used in Reference [15], and produces a clear demarcation between the droplet core region, and far-field region, where the droplet profile decays to zero exponentially (we take this opportunity here to clarify the procedure used to compute x ref , which was not clear in our previous work). Using either the particle method, or the finite-difference method with x ref chosen is equivalent: both descriptions reproduce the t 1∕7 behaviour for the contact line.
The finding that the Geometric Thin-Film Equation satisfies Tanner's Law of droplet spreading indicates that the regularized model is capturing the large-scale physics in the  Fig. 5 we produce a space-time plot of f ( , t) -this is seen to relax to a constant profile at late times as t → ∞.
Furthermore, the profile of f ( , t) at fixed t (t large) can be compared with a similarity solution of the unregularized problem, f 3 f ��� = f ∕7 (cf. Eq. (14)). This ordinary differential equation is then solved with the shooting method together with appropriate initial conditions [15]. The results are shown in Fig. 6.
This figure therefore shows that the Geometric Thin-Film Equation describes the expected large-scale droplet-spreading physics in the droplet core. Where the classical Thin-Film Equation breaks down at the contact line, the height profile of the Geometric version decays smoothly to zero.

Rigorous error and performance analysis
We analyse the truncation error associated with the standard fully-implicit finite-difference method and the particle method. As such, let h denote the exact solution of Eq. (40), and let h x denote the numerical solution with step size x (finite-difference method), or number of particles N = 2L∕ x (particle method). We assume that the error ‖h − h x ‖ depends smoothly on x , then for some constant C and p. Since h is unknown, we instead compute Using the triangle inequality, it can be shown that We take the natural logarithm on both sides of Eq. (59); this gives:  7 Convergence plot of the finite-difference method and the particle method for the partial wetting case. Both lines have a slope of 2.0 on the log-log plot Thus, the rate of convergence (or the order of accuracy) of the numerical method is p; p can be computed from the numerical simulations as the slope of the log-log plot between the error and the grid spacing x. Figure 7 shows the rate of convergence of the finite-difference method and the particle method -here we use the L 1 norm applied to Eqs. (59)-(60) (our choice of norm is obtained because the particle solutions are expected to converge weakly in an L 1 function space, see Reference [36]). The finite-difference method is implemented with a step size of t = 0.01 , while the particle method uses the ODE45 solver in MAT-LAB, hence an adaptive time step. Both methods use the same numerical parameter of, = 0.05 , with final time T = 1 , and periodic boundary condition on the spatial domain x ∈ [−1, 1] . From this figure, both the particle method and the standard finite-difference method are estimated to be second-order accurate in the spatial domain, albeit that the absolute error is smaller for the finite-difference method.
Furthermore, in Fig. 8 we evaluate the execution time of the different numerical methods to see if any one method outperforms the rest. A comparison of the average execution time over 10 runs between the finite-difference method, the direct implementation of the particle method (computational complexity O(N 2 ) ), and the fast implementation particle method (computational complexity O(N)). The numerical parameters used are the same as the one used in the convergence analysis and the calculations are performed on an Intel i7-9750H with 6 hyper-threaded cores.
Thus, although the fast-particle method is less accurate for fixed N (Fig. 7), the fastparticle method is more efficient (Fig. 8). Therefore, the a fixed execution time (but varying N), the fast-particle-method and the finite-difference method achieve a similar level of accuracy, albeit that the finite-difference method still has the edge. However, beyond this standard performance analysis, the particle method still retains a key benefit, namely that it is free from finite-size effects. Fig. 8 a Performance of the finite-difference method, the direct implementation of the particle method, and the fast implementation of the particle method 1 3

Partial wetting
In this section, we extend the Geometric Thin-Film Equation to the case of partial wetting, where a droplet on a substrate spreads initially before assuming an equilibrium shape. This requires the inclusion of an extra, stabilizing term, to Eq. (40). We derive this additional term. Then, we construct an analytical solution for the equilibrium droplet shape. Finally, we use both the finite-difference method and the particle method to simulate transient droplet spreading, up to the point where the droplet assumes its equilibrium shape.

Theoretical formulation
The starting-point for the theoretical formulation is to consider an unregularized description of the droplet, with h(x) as the droplet profile. Then, the unregularized energy associated with a droplet of radius r is: Here, la is the surface tension between the air and the liquid droplet (previously referred to as ), ls is the surface tension between the liquid and the substrate, and as is the surface tension between the air and the substrate; S is an arbitrary lengthscale denoting the extent of the system in the lateral direction. In the longwave limit, x , and Eq. (61) becomes: The three surface-tension coefficients are related via the Laplace-Young condition, where eq is the equilibrium contact angle. Thus, Eq. (62) can be re-written as Now, inspired by the replacements h → h in Sect. 3, we propose herein a regularized energy, where w is an estimate of the size of the droplet footprint, based on the interfacial profile h, and on the smoothened interfacial profile, h = * h. We estimate the size of the droplet footprint as x dx + 2r la + ls − as + Const.

Non-dimensionalization
We non-dimensionalize Eq. (68) using the lengthscale h 0 = √ A 0 tan eq in the vertical direction and L =

√
A 0 ∕ tan eq in the lateral direction. Thus, h is made dimensionless on h 0 , x is made dimensionless on L, and time is made dimensionless on the capillary timescale (3 L∕ )(tan 3 eq ) . In dimensionless variables, Eq. (68) now reads: Also, in dimensionless variables, ∫ ∞ −∞ h(x, t)dx = 1 . Also in this context, the ratio = h 0 ∕L is precisely tan eq ; strictly speaking therefore, eq should be small, for inertial effects to be negligible, and hence, for the lubrication theory underlying Eq. (69) to be valid.

Equilibrium solution
Equation (69) has an equilibrium solution with h∕ t = 0 . In this limiting case, Eq. (69) reduces to where 2 is a positive constant, Equation (70) has an analytical solution, parametrized by , and by a radius r:

Correspondingly,
Here, B 1 , B 2 , C 1 , and C 2 are constants of integration. These constants are fixed by imposing continuity of h, h x , h xx at x = r , and also by imposing ∫ ∞ −∞ h dx = 1 . These give four conditions in four unknowns. Hence, {B 1 , B 2 , C 1 , C 2 } are fixed in terms of r. The value of r is in turn fixed by imposing continuity of h xxx at x = r , this gives a simple rootfinding condition: The details of this calculation are provided in Appendix B. Equation (74) gives r as a function of . However, is not arbitrary, but is instead fixed by its own rootfinding condition (Eq. (71)). Although this procedure is somewhat involved, the point remains: the Geometric Thin-Film Equation with partial wetting admits an analytical equilibrium solution in terms of elementary functions (via Eqs. (72)-(73)). Furthermore, the elementary solution (72) for h(x) coincides with the expression for a parabolic-cap droplet in the core region x → 0, It is noted however, that this expression is valid only in the core region x → 0 . In particular, the analytical solution h(x) does not converge to a parabolic-cap profile in the limit as → 0 ; rather, h(x) maintains its cosine shape in this limit. Therefore, in order to reproduce the parabolic-cap profile in the limit as → 0 , it would be necessary to make a more judicious choice for estimate of the width of the droplet footprint in Eq. (67). Finding the  optimal choice for the functional form of the width of the droplet footprint will be the subject of future work. The equilibrium contact angle is now computed as: where h x is expressed in dimensionless variables and x ref > 0 is a reference point. Since tan eq = in the chosen dimensionless variables, this requires: Following the procedure described in Sect. 5, we choose the reference point x ref to be the positive value of x which maximizes |h x | , thus x ref = ∕(2 ) . Thus, we require B 1 = 1 . But B 1 and depend parametrically on , hence we require B 1 ( ) ( ) = 1 . This therefore fixes  Table 1.

Results
For the simulations of partial wetting, we use the initial condition with r 0 = 0.5 , and ∫ ∞ −∞ h(x, t = 0) dx = 1 . We also take = 0.05 . We use both the particle method and the finite-difference method: the results are the same in each case. In Fig. 10a, we show a space-time plot of the solution for the partial wetting case up to t = 10 ; panel (b) shows a snapshot of the droplet profile at t = 10 . Figure 11 is based exclusively on the particle method: here we show a log-log plot of the particle trajectories. At intermediate times, the particle trajectories are parallel to the path x = t 1∕7 before attaining a steady state at late times. Thus, in the case of partial wetting, the system obeys Tanner's law at intermediate times, until at late times, the partial wetting stabilizes the droplet and it assumes its equilibrium shape.
The foregoing statement that the droplet spreading obeys Tanner's Law at intermediate times until the onset of equilibrium also applies to the Cox-Voinov Law [9] of droplet spreading, which in equation form is [ (t)] 3 = [ eq ] 3 + ċx cl log(x cl ∕d) . Here, (t) is the dynamic contact angle, which is obtained from the slope of the interface profile h(x, t) at some appropriate location x, and c and d are constants. In order to validate the applicability of the Cox-Voinov law to the droplet spreading within the framework of the Geometric Thin-Film Equation, we operationally define the contact angle as (t) = max x [− x h(x, t)]) . The tangent line to h(x, t) at argmax x [− x h(x, t)] is constructed, and the contact line is then taken to be the intersection of this tangent line with the x-axis. This is the same procedure as in Sect. 5. A plot of [ (t)] 3 constructed in this way is shown in Fig. 12. Shown also is a plot 1 + ċx cl log(x cl ∕d) -here, eq = 1 , and c and d are best-fit constants. Overall, there |x| > r 0 . Fig. 13 Convergence plot of the finite-difference method and the particle method for the complete wetting case. Both lines have a slope of 2.0 on the log-log plot 1 3 is good agreement between the two curves at intermediate times -consistently with the behaviour of the trajectories Fig. 10. There is some disagreement at late times, however, this may be expected, in view of the somewhat imprecise operational definition of (t) and x cl .

Rigorous error and performance analysis
We carry out a rigorous error analysis of both the fully-implicit finite-difference method, and the particle method, for the case of partial wetting. Because there is an analytical, equilibrium solution valid at late times, we analyze the results of the numerical simulations at such late time (specifically, t = 100 ), when the numerical solutions are close to the equilibrium state. In this case, the equation for the rate of convergence of the numerical method is simply where h denotes the analytical equilibrium solution, and h x denotes the numerical equilibrium solution (or what amounts to the same, the numerical solution at t = 100 ). Figure 13 shows the rate of convergence for the finite-difference method and the particle method. Both methods have = 0.05 and were performed on the spatial domain x ∈ [−2, 2] for various spatial resolutions x . We have used the MATLAB solver ODE15s for the particle method. Again, we observed both the finite-difference method and the particle method to be second-order accurate in the spatial domain. Finally, we have looked at the performance of the different numerical methods (fullyimplicit finite-difference method, particle method, and 'fast' particle method): the results are similar to what was observed in the case of complete wetting (Fig. 14). Fig. 14 Performance of the finite-difference method, the direct implementation of the particle method, and the fast implementation of the particle method (partial wetting)

Conclusions
Summarizing, we have introduced a new mathematical model to describe contact-line motion which has the same regime of validity as conventional lubrication theory. The model involves using a 'smooth' interface profile h and a sharp interface profile h. The smooth interface profile h is connected to the sharp interface profile h via convolution, while h is defined through an evolution equation which couples both interface profiles, and which drives the droplet spreading. It is possible to assign a physical interpretation to the model: the equation for h is a model with missing small-scale physics (below a lengthscale ), while h is a complete model of the physics at lengthscales greater than . By minimizing the difference between the two descriptions, a convolution relation between h and h emerges. Furthermore, by formulating the evolution equation so that it involves both descriptions of the interface profile, the contactline singularity is resolved. There is no need to model the missing small-scale physics explicitly: instead, it is parametrized through the lengthscale .
A further advantage of the new formulation is that it naturally suggests a mesh-free numerical method for the purpose of simulating the model numerically. Furthermore, the mathematical model involves equations which are non-stiff, and are straightforward to solve numerically (we include a repository of the code as part of this work; see Reference [37]). We have conducted numerical simulations for both complete wetting and partial wetting using this new mesh-free method (as well as a conventional finite-difference method), and have found that the model describes well the physics of droplet spreading -including Tanner's Law for the evolution of the contact line. Remarkably, in the case of partial wetting, the model also admits a simple analytical solution for the equilibrium profile.
Beyond droplet spreading, the model may find further applications in describing families of droplets (the analytical equilibrium solution already possesses such 'multiple-droplet' solutions), multi-component systems, and problems in droplet evaporation. The model's intrinsically non-stiff equations (as well as the absence of a precursor film extending to infinity) may in the future simplify the description of such complex physical systems.
Eq. (40). We also work in dimensionless variables, such that the prefacor ∕(3 ) which was present in Eq. (40) is rescaled to one. Eq. (76) is discretized in space using a standard finite-difference method. The derivatives are approximated using standard finite-difference stencils (central differencing, first and second-order derivatives second-order accurate in space). The value of h is therefore known only at the discrete grid points; we collect these values into a vector h . Similarly, Eq. (76) is discretized in time such that the values of h are known only at discrete points in time; the fixed interval between each such discrete point in time is the timesetep t . The values of the vector h at a particular discrete point in time are therefore labelled as h n . Thus, in discrete form, Eq. (76) becomes: Here, D 1 and D 3 are sparse matrices corresponding to the centred finite-difference approximation of the first and third spatial derivatives, respectively. A similar notation would hold for D 2 , which corresponds to the centred finite-difference approximation of the second spatial derivative. In this way, K is a convolution matrix, which is constructed as K = (1 − 2 D 2 ) −2 . Furthermore, we use the ⊙ notation in Eq. (77) to denote the pointwise multiplication of vectors, e.g. if a and b are N-dimensional column vectors with entries a = (a 1 , … , a N ) T and b = (b 1 , … , b N ) T , then a ⊙ b is also an N-dimensional row vector with entries a ⊙ b = (a 1 b 1 , … , a N b N ) T . This notation is useful for the exposition of the numerical method, as the numerical code for solving Eq. (76) is vectorized, and makes use of precisely this pointwise vector multiplication.
Crucially, in the discretized Eq. (77) we have treated all terms on the right-hand-side fully implicitly ( h n+1 instead of h n ). This treatment amounts to a backwards-Euler temporal discretization. A fully explicit treatment (Forward Euler, with h n instead of h n+1 ) would suffer from a severe hyperdiffusion-limited CFL constraint t < C x 4 , where x is the grid spacing, and C is an O(1) constant. A suitable combination of h n and h n+1 on the right-hand side of Eq. (77) (amounting to a Crank-Nicolson scheme) is also a possibility, although due to its simplicity, the backwards-Euler scheme is preferred here. Although the choice of Backward-Euler method makes the numerical code unconditionally stable, it means that at each time-step, a set of nonlinear equations for h n+1 must be solved. We therefore describe how these nonlinear equations are solved.
We introduce where D 1 = KD 1 , D 3 = D 3 K , and v = K −1v . Then, solving Eq. (77) for h n+1 given h n becomes a minimization problem on We employ a Newton-Linesearch method to solve for h n+1 . This requires us to compute the Jacobian of F(v) , which is given by We initialize the guess with v 0 =h n . Then, for each iteration I, the descent direction is given by (80) J(v) = I + t (v 2 ⊙ K −1 + 2v ⊙ v) ⊙ D 3v + (v ⊙v 2 ) ⊙ D 3 . and the improved guess is updated using where I is the optimum step size in the direction v I The iterative process is continued until the residual f (v I ) is sufficiently small, and we set the solution for the next time step by letting h n+1 =v I .

A.2 Fast summation algorithm for particle method
Here, we are interested in solving the Geometric Thin-Film Eq. (40) using the particle method. The particle method requires us to compute the evolution equations for the particle trajectories, these are given in Eq. (52). To solve Eq. (52), one is required to evaluate for i = 1, … , N , at each time step. A direct evaluation would require an operational cost of O(N 2 ) . Now, suppose that x i < x j for all i < j , then the contribution from particles of either side of x i becomes a degenerate case of the fast multipole method. We start by defining new variables Since x i − x j > 0 for j = 1, … , i − 1 , we can drop the absolute value in the function , Similarly, we obtain for S 2 , f v I +v I . (87)

B Derivation of the equilibrium solution for partial wetting
In this Appendix, we look at the equilibrium solution for the droplet in the case of partial wetting. The equilibrium solution was outlined in Sect. 6, and involved expressions for h(x) and h(x) in terms of elementary functions. We begin by recalling these expressions: Here, B 1 , B 2 , C 1 , and C 2 (together with r) are constants of integration, which can be obtained via the methods outlined in this Appendix, although we focus herein on deriving the rootfinding condition (74) previously introduced in the main part of the paper.
We start by noticing that the set of boundary conditions  (2), (3) and (4), we get the following equations: Taking (96 − 2(97) + (98)) yields boundary condition (1). Therefore choosing any three out of the four boundary conditions yields the same system of equations.
Next, consider the following set of boundary conditions for expression (95): From these boundary conditions, one can obtain the following two linearly independent equations (99) (101) 2 + 2 2 = C 1 + rC 2 + C 2 C 1 + rC 2 − C 2 , Fig. 15 Effect on the spreading rate of the contact line x cl of varying for a the finite difference method and b the particle method In matrix form, this can be expressed as Requiring the determinant of the matrix to be zero yields Eq. (74). From the rootfinding condition (74), a value for the parameter is found in the main paper, via Eq. (71), recalled here as 2 = 2 ∕⟨h, h⟩ 2 . For this purpose, we lastly provide here an analytical expression for ⟨h,h⟩: where B 1 = B 1 (1 + 2 2 ) 2 .
C Sensitivity analysis with respect to the parameter Į n this appendix, we look at the sensitivity of the numerical results for the case considered in Sect. 5. In Fig. 15, we show the effect on the macroscopic contact line of varying -these effects are small, and diminish as is decreased.