2.1. Numerics

Our numerical scheme follows others used previously to solve the highly nonlinear thin film equation (Dupont et al. Reference Dupont, Goldstein, Kadanoff and Zhou1993; Eggers & Fontelos Reference Eggers and Fontelos2015): we use a fully implicit finite difference scheme with a staggered grid. To that end, (1.2) is split as

(2.5) \begin{equation} h_t + (h^{n-1} f)_x = 0, \quad f = h h_{xxx}, \end{equation}

and the interval $[0,1]$ is divided into grid points $x_i, i=1\dots k$ , with $x_1 = 0$ and $x_k = 1$ ; $f_i, i=1\dots k-1$ is defined at the midpoints. The profile $h_i$ is used to compute $p_k = h_{xx}$ at the grid points. At the right-hand end, we impose symmetry of $h$ , and antisymmetry of $f$ . The values $f_i$ are calculated from $p_i$ using centred differences, where $p_i$ is calculated with a 5-point scheme, valid for arbitrary grid spacings. The values of $f$ and $f_x$ , needed at the grid points for the first equation, are also calculated with second-order accuracy.

Near the contact line, we impose the expansion (2.2), where the contact angle $\theta$ is a separate variable (which takes the place of $f_1$ ), and which evolves according to (2.3). The value of $f_1$ , between $x_1$ and $x_2$ , is calculated from (2.2). Likewise, for a given $\theta$ , the coefficients $b$ and $c$ are found from equating (2.2) to $h_{2}$ and $h_{3}$ at $x_{2}$ and $x_{3}$ , respectively. Then, if we let $\varDelta = x_2$ , and $x_3 = 2\varDelta$ , we find that

(2.6) \begin{equation} b = \frac {1}{\varDelta \ln 2}\left (\frac {\theta }{2} - \frac {h_3-4h_2}{4\varDelta }\right ), \quad c = -b\ln \varDelta + \frac {1}{\varDelta }\left (\frac {h_2}{\varDelta }-\theta \right ). \end{equation}

In summary, the variables are $h_2, \dots , h_{k-1}$ , $f_2, \dots , f_{k-1}$ and $\theta$ : a total of $2k-4+1 = 2k-3$ . The equations are the first of (2.5), evaluated at $x_2, \dots , x_{k-1}$ , and the second of (2.5), evaluated at $x_{1/2}, \dots , x_{k-1/2}$ , with the extra equation (2.3).

Figure 3.

The change in contact angle $\dot {\theta }$ as a function of time from a numerical solution of (1.2) with $n=3$ and initial condition $h_0(x) = \bar {h}\sin (\pi x/2)$ , $\bar {h} = 0.1$ . Thus $V = 4\bar {h}/\pi$ . For long times, $\dot {\theta }$ decays such as ${\rm{e}}^{-0.57 t}$ to zero (see § 5); this exponential decay is shown as the dashed line. For short times, on the other hand, $\theta$ appears to remain constant for some time ( $\dot {\theta } = 0$ ); we show that, instead, $\dot {\theta }$ behaves according to $\displaystyle {\rm{e}}^{-3/(2t^{1/3})}$ .

We use a strongly graded grid, with the smallest (constant) grid spacing near $x=0$ , and which is slowly increased away from $x=0$ , until a maximum spacing of $10^{-3}$ is reached. Our numerical scheme is fully implicit, and second order in time, using a step-halving method. Comparison between the two steps serves to adjust the time step, in order to maintain sufficient temporal resolution.

3.1. Effective linear dynamics

We assume that the initial condition can be expanded into a power series of the form

(3.1) \begin{equation} h_0 = a_0 x + b_0 x^2 + c_0 x^3 + O(x^4). \end{equation}

Since $x$ is small close to the contact line, we can view $c_0 x^3$ as a small perturbation to the equilibrium profile $h_{eq} = a_0x + b_0 x^2$ , around which we linearise by writing $h(x,t) = h_{eq}(x) + \delta h(x,t)$ . We thus consider the linear equation

(3.2) \begin{equation} \delta h_t + \left (h_{eq}^n\delta h_{xxx}\right )_x = 0, \end{equation}

so we have confirmed numerically that the early time behaviour of (3.2) is indeed the same as that of the full equation.

Figure 4.

The change in contact angle $\dot {\theta }(t)$ as function of time. Solid line: solution of (1.2) with $n = 3$ and initial condition $h_0(x) = x(2 - x) + 0.1 x^2(2-x)^2$ ; contact angle found from (2.3). Boundary conditions are $h(0,t) = 0$ and $h_x(1,t) = h_{xxx}(1,t) = 0$ . Long-dashed line: solution of the linearised leading-order equation (3.2) with $h_{eq} = 2x$ and initial condition $\delta h(x,0)= 0.1 x^2(2-x)^2$ . Boundary conditions for the linearised problem are $\delta h(0,t) = 0$ and $\delta h(1,t) = \delta h(1,0), \delta h_{xx}(1,t) = \delta h_{xx}(1,0)$ . Dot-dashed line: equation (3.2) with $h_{eq} = 2x$ and initial condition $\delta h(x,0) = -0.4 x^3$ , with boundary conditions $\delta h(0,t) = 0$ and $\delta h(1,t) = -0.4, \delta h_{xx}(1,t) = 2.4$ . (a) All three curves; (b) a detail of the dashed and dot-dashed lines, with the amplitude of the latter adjusted.

In a second step, we argue that, close to the contact line, $h_{eq}^n$ is dominated by the linear term $h_{eq} \approx a_0 x$ ; hence, we obtain

(3.3) \begin{align} \delta h_t + \left (a_0^nx^n\delta h_{xxx}\right )_{x} = 0, \end{align}

which is a linear and homogeneous problem. Since $a_0 x + \delta h(x,t)$ is no longer the profile of a droplet, we can no longer impose the same symmetry conditions. Instead, we argue that, as far as representing the far-field conditions are concerned, we can impose that the thickness and curvature of the drop profile changes little when imposed on the scale of the drop radius (which is normalised to unity), and viewed over the short time scale on which the contact line region is evolving. We thus impose

(3.4) \begin{equation} \delta h(1,t) = \delta h(1,0), \quad \delta h_{xx}(1,t) = \delta h_{xx}(1,0) . \end{equation}

This idea is tested in figure 4 for the critical case $n = 3$ . In figure 4(a), using the solid line, we show a simulation of the fully nonlinear equation (1.2) with contact angle condition (2.3) and boundary conditions $h(0,t) = 0$ , $h_x(1,t) = h_{xxx}(1,t) = 0$ ; for the initial profile, we take $h_0(x) = x(2 - x) + 0.1 x^2(2-x)^2$ , where $h_{eq}(x) = x(2-x)$ is the equilibrium profile, and $\delta h(x,0) = 0.1 x^2(2-x)^2$ the initial perturbation. This solution is compared with the linearised equation (3.2), in which we have also replaced $h_{eq}(x)$ by the leading order expansion $h_{eq} = 2x$ of the initial condition. As the initial condition of the linearised problem, we once more take $\delta h(x,0) = 0.1 x^2(2-x)^2$ and apply the boundary conditions (3.4).

The linearised form of the evolution equation for the contact angle (2.3) is

(3.5) \begin{equation} \dot {\theta } = -2(2+\delta )(1+\delta )b, \end{equation}

with $b$ once more defined by (2.2), but with $\theta$ now representing the deviation of the contact angle from its initial value. The resulting solution for $\dot {\theta }$ is seen in figure 4(a) to agree very well with the full nonlinear solution, at least for $t \lesssim 3\times 10^{-5}$ . In a third step, we further approximate the initial condition of the linearised problem as $\delta h(x,0) = -0.4 x^3$ , derived from the cubic term of $\delta h(x,0) = 0.1 x^2(2-x)^2$ . The boundary conditions are once more given by (3.4).

The resulting solution $\dot {\theta }$ is now shown as the dot-dashed line in figure 4(a). To appreciate the agreement with the dashed line, which includes higher-order terms in the initial condition, we show a zoomed-in version of the same curves in figure 4(b). The agreement for earlier times is such that the two curves can almost not be distinguished; however, we have allowed for a small change in amplitude of the dot-dashed line to reach optimal agreement. This is natural, since we are solving a linear equation, and the effective driving amplitude will have changed slightly on account of the higher order terms.

In summary, up to an overall constant, and after appropriate rescaling, we may describe the short-time dynamics of the full non linear contact line problem by

(3.6) \begin{equation} \delta h_t + \left (x^n\delta h_{xxx}\right )_{x} = 0, \end{equation}

with initial condition $\delta h = x^3$ . The boundary conditions are

(3.7) \begin{equation} \delta h(0,t) = 0, \quad \delta h_x(0,t) = \theta , \quad \quad \delta h(1,t) = 1, \quad \delta h_{xx}(1,t) = 6, \end{equation}

where the dynamics of $\theta$ is given by (3.5). To remove unnecessary constants, in our analysis of (3.6) we found it useful to take the second derivative of (3.6), to obtain an equation for the curvature $\kappa = \delta h_{xx}$ of the perturbation

(3.8) \begin{equation} \kappa _t + \left (x^{3-\delta } \kappa _x\right )_{xxx} = 0, \end{equation}

with initial condition $\kappa = 6 x$ . From (2.2), we find that, for small $x$ ,

(3.9) \begin{equation} \kappa = \frac {(2+\delta )(1+\delta )b}{\delta } x^{\delta } + 2\left (c - \frac {b}{\delta }\right ),\; \delta \ne 0,\quad \kappa = 2 b \ln x + 3b + 2c,\; \delta = 0. \end{equation}

Inserting the initial condition into (3.8), we obtain $\kappa _t = -6(3+\delta )(2+\delta )(1+\delta )x^{\delta }$ , so that at short times the solution is

(3.10) \begin{align} \kappa = 6x - 6(3+\delta )(2+\delta )(1+\delta )x^{\delta } t. \end{align}

This means that new terms of the form (3.9) are generated immediately, and the contact angle must change. We can expect a local solution driven by the dynamics of the contact line alone. If, on the other hand, $\delta = 0$ , $\kappa = 6x - 36 t$ solves (3.8) exactly and is compatible with the initial conditions. The dynamics is now driven by the fact that the solution is no longer compatible with the boundary condition at the other end: a non-local process. As a result, the contact angle hardly changes initially, since the driving is extremely weak. We start with the local generic, case $\delta \ne 0$ .

3.2. The contact angle, $\delta \ne 0$

We begin with the case where $n$ does not equal the generic value for a liquid film on a solid substrate. We are looking for similarity solutions of the form $\kappa = t^{\gamma } P(\xi )$ (Bernis et al. Reference Bernis, Hulshof and King2000), where $\xi = x/t^{1/(1+\delta )}$ , so that (3.8) is satisfied. For large $\xi$ , this solution has to match $\kappa = 6x$ ; this implies that $\gamma = 1/(1+\delta )$ , so that

(3.11) \begin{equation} \kappa = t^{1/(1+\delta )} P(\xi ), \quad \xi = x/t^{1/(1+\delta )} \end{equation}

is the form of the similarity solution, with $P(\xi ) = 6\xi$ for large arguments. The similarity equation becomes

(3.12) \begin{equation} \frac {P - \xi P_{\xi }}{1+\delta } + \left (\xi ^{3-\delta } P_{\xi }\right )_{\xi \xi \xi } = 0. \end{equation}

Four linearly independent solutions $P_1,P_2,P_3$ and $P_4$ to (3.12) can be found in terms of generalised hypergeometric functions (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010), as given in Appendix A, cf. (A1)–(A4); they are real for $\xi \gt 0$ . Very similar solutions to an elastic fourth-order problem have been described in (Stone & Duprat Reference Stone and Duprat2016). The solutions $P_2$ and $P_4$ are singular at the origin, and on account of (3.9) have to be excluded. The remaining solutions $P_1$ and $P_3$ grow exponentially at infinity, while we demand $P(\xi ) = 6\xi$ . As shown in detail in Appendix A, we can cancel the exponential growth by superimposing $P_1$ and $P_3$ . Different expressions apply for $\delta \gt 0$ and for $\delta \lt 0$ , respectively. In summary, the solution to (3.12) satisfying all required conditions can be written in the form

(3.13) \begin{equation} P(\xi ) = P_0 \left (P_3(\xi ) - r P_1(\xi )\right ), \end{equation}

where $P_0$ (cf. (A15)) is chosen such that (3.13) matches $P = 6\xi$ ; the constant $r$ is defined in (A11).

Figure 5.

(a) Numerical results for the exponent of $\dot {\theta }$ , based on (3.6), compared with the analytical result (3.15) (solid line). (b) Numerical result for the first term in (3.12) as a function of $\ln \xi$ (black solid line, $\delta = 0.1$ ) compared with the analytical solution (3.13) (red dotted line) .

For small $\xi$ , the solution behaves as

(3.14) \begin{equation} P \approx P_0 \left (\xi ^{\delta } - r\right ). \end{equation}

Comparing with (3.9) and using (3.5), we find that

(3.15) \begin{equation} \dot {\theta } = -2 a_0^3 \delta P_0 t^{\frac {1-\delta }{1+\delta }} \end{equation}

and

(3.16) \begin{equation} c = \frac {P_0}{(2+\delta )(1+\delta )} t^{\frac {1-\delta }{1+\delta }} - \frac {P_0 r}{2} t^{\frac {1}{1+\delta }}. \end{equation}

This means that, as the drop starts to move, the contact angle changes immediately according to (3.15). This is tested by comparison with numerical simulations of the linearised equation (3.6) for $n\neq 3$ . In figure 5(a), we show the exponent of $\dot {\theta }$ . The solid line is $(1-\delta )/(1+\delta )$ is found from (3.15). The symbols are the exponent as found from numerical simulations. In figure 5(b), the left-hand side of (3.12) (as determined numerically, black solid line), is tested against the analytical solution (3.13) (red solid line), and perfect agreement is found.

Strictly speaking, the similarity solution (3.13) is only one of an infinite sequence of solutions of higher order, but which are unstable (Eggers & Fontelos Reference Eggers and Fontelos2015). If the initial condition happens to be such that the coefficient of $x^3$ vanishes exactly, then one must consider the next order, $x^4$ , or an exponent $m\geqslant 3$ in general. However, an arbitrarily small perturbation will generate a term proportional to $x^3$ , rendering solutions of higher order unstable. In that case, the asymptotic behaviour is $P \propto \xi ^{m-2}$ , and the (higher-order) similarity solution is of the form

(3.17) \begin{equation} \kappa = t^{(m-2)/(1+\delta )} P(\xi ), \quad \xi = x/t^{1/(1+\delta )}, \end{equation}

with similarity equation

(3.18) \begin{equation} \frac {(m-2)P - \xi P_{\xi }}{1+\delta } + \left (\xi ^{3-\delta } P_{\xi }\right )_{\xi \xi \xi } = 0. \end{equation}

The two solutions with the correct behaviour at the origin (given in Appendix A, cf. (A16)) can again be superimposed to generate a solution that behaves like $\xi ^{m-2}$ at infinity, as required.

4.1. Similarity solution near $x=1$ : the pulse solution

To understand the origin of the oscillations that are generated at the right-hand end of the interval, we consider the similarity solution

(4.8) \begin{equation} \Delta \kappa = t^{\alpha }\phi (\zeta ), \quad \zeta = \frac {1-x}{t^{1/4}}, \end{equation}

located at $x=1$ . The exponent $1/4$ follows from (4.2) considering that $x\approx 1$ , while $\alpha = 1/2$ follows from $\delta h_t = -\delta h_{xxxx} = -\Delta \kappa _{xx}$ and $\delta h(1,t) = 18t$ , so that $\Delta \kappa _{xx}(1,t) = -18$ . We thus have

(4.9) \begin{equation} \Delta \kappa = t^{1/2}\phi (\zeta ), \quad \zeta = \frac {1-x}{t^{1/4}}, \end{equation}

where $\phi$ satisfies the similarity equation

(4.10) \begin{equation} \frac {\phi }{2} - \frac {\zeta \phi '}{4} + \phi ^{iv} = 0. \end{equation}

Figure 7.

The profile $\phi$ as found from computing $\Delta \kappa$ , rescaled according to the similarity solution (4.9), for $t = 3.2\cdot 10^{-8}$ (black line) and $t = 10^{-7}$ (red line). The dashed line is the solution (4.11).

The four linearly independent solutions $\phi _1,\phi _2,\phi _3$ and $\phi _4$ can once more be written in terms of generalised hypergeometric functions, as given in (B1); of them, $\phi _3$ does not satisfy $\phi (0) = 0$ , which is required to conform with the boundary condition $\kappa (1,t) = 6$ , and drops out. In addition, we want $\phi$ to decay at infinity. As explained in more detail in Appendix B, the ratio between the two remaining hypergeometric functions is fixed to $r = \sqrt {2}\pi /(24\varGamma ^2(3/4))$ , so that the exponential growth cancels. In a second step, we fix the amplitude of $\phi _3$ to $r_1 = \pi /(4\varGamma (3/4))$ to remove the term growing quadratically at infinity. The remaining combination then decays exponentially at infinity. This leaves us with the solution

(4.11) \begin{align} \phi = A\left [\zeta {}_1\textrm {F}_3\left (-\frac {1}{4}; \frac {1}{2}, \frac {3}{4}, \frac {5}{4},\frac {\zeta ^4}{256}\right ) + \frac {\sqrt {2}\pi }{24\varGamma ^2\left (\frac {3}{4}\right )} \zeta ^3 {}_1\textrm {F}_3\left (\frac {1}{4}; \frac {5}{4}, \frac {3}{2}, \frac {7}{4},\frac {\zeta ^4}{256}\right ) - r_1 \zeta ^2\right ], \\[-22pt] \nonumber \end{align}

which behaves according to $\phi \approx A\zeta - r_1\zeta ^2$ for small $\zeta$ . To finally fix $A$ we observe that the boundary condition at $x=1$ is $-18 = \Delta \kappa _{xx}(1,t) = \phi ''(0) = -2 A r_1$ , from which $A=9/r_1=36\varGamma (3/4)/\pi \approx 14.0422$ . This fixes all parameters of the pulse solution (4.11), which is shown as the dashed line in figure 7. Clearly, very good agreement with numerical simulation is found, without an adjustable parameter.

To understand the behaviour for large $\zeta$ (away from the corner), we make the WKB ansatz $\phi \propto {\rm{e}}^{S(\zeta )}$ , as described in chapter 10 of Bender and Orszag (Reference Bender and Orszag1978). Inserting into (4.10), to leading order we get $-\zeta S'/4 + S'^4 =0$ , or

(4.12) \begin{align} S' = \left (\frac {\zeta }{4}\right )^{1/3}\begin{cases} 1, \\ -1/2+\sqrt {3}{\rm{i}}/2, \\ -1/2-\sqrt {3}{\rm{i}}/2, \\ \end{cases} \end{align}

where the last two solutions are the relevant case which decays for large $\zeta$ . Integrating, we obtain

(4.13) \begin{equation} S_{1/2} = S_0 - \frac {3\zeta ^{4/3}}{2\times 4^{4/3}}\left(1\pm \sqrt {3}\rm{i}\right). \end{equation}

To capture algebraic corrections, we need to go to next order: i.e.

(4.14) \begin{align} S' = \left (\frac {\zeta }{4}\right )^{1/3} \left (-\frac {1}{2}\pm \frac {\sqrt {3}\rm{i}}{2}\right ) + \epsilon , \end{align}

for which the leading terms in (4.10) are

(4.15) \begin{align} \frac {1}{2}-\frac {\zeta S^{\prime }}{4} + S^{\prime 4} + 6 S^{\prime \prime } S^{\prime 2} = 0, \end{align}

so that $\epsilon = -4/(3\zeta )$ . The complex roots can be combined to cosine and sine modes, so that the asymptotic solution finally becomes

(4.16) \begin{equation} \phi \approx A_p\zeta ^{-4/3} {\rm{e}}^{{-\frac {3\zeta ^{4/3}}{2\cdot 4^{4/3}}}} \cos \left (\frac {3^{3/2}\zeta ^{4/3}}{2\times 4^{4/3}} - \phi _p\right ). \end{equation}

Here $\phi _p = 0.561$ is a phase factor, and the amplitude is $A_p = 3.6599A \approx 51.39316$ , both found numerically, based on a numerical evaluation of (4.11) using Maple. This means that $\Delta \kappa$ near the right-hand end of the domain is

(4.17) \begin{equation} \Delta \kappa \approx A_p t^{1/2}\zeta ^{-4/3}\exp \left [-\frac {3\zeta ^{4/3}}{2\times 4^{4/3}}\right ] \cos \left (\frac {3^{3/2}\zeta ^{4/3}}{2\times 4^{4/3}} - \phi _p\right ). \end{equation}

4.2. Similarity solution near the contact line

Note that (4.13) implies that

(4.18) \begin{equation} S_{1/2} \propto (1 \pm \sqrt {3}{\rm{i}}) \frac {(1-x)^{4/3}}{t^{1/3}}, \end{equation}

which suggests a similarity solution of the form

(4.19) \begin{equation} \Delta \kappa = {\rm{e}}^{-\frac {\mu }{\alpha t^{\alpha }}} A(\xi ), \quad \xi = \frac {x}{t^{1 + \alpha }}, \end{equation}

where $\mu$ is expected to be complex and $\alpha$ is real. We will find $\alpha =1/3$ corresponding to (4.13), but first continue the calculation for general $\alpha$ ; note that $A$ is also complex.

Inserting (4.19) into (4.2), we obtain

(4.20) \begin{align} \mu A - (1+\alpha ) t^{\alpha }\xi A' + \left (\xi ^3 A'\right )^{\prime \prime \prime} = 0, \end{align}

which, for $\alpha \gt 0$ and $t\rightarrow 0$ , simplifies to

(4.21) \begin{equation} \mu A + \left (\xi ^3 A'\right )^{\prime \prime \prime} = 0. \end{equation}

As given in (C1), solutions can be written in terms of the MeijerG function (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010), which is defined in Appendix C. Of the four solutions to (4.21), $A_3$ and $A_4$ are singular at the origin, and are therefore excluded. Instead, we are looking for solutions that are consistent with (4.3) at the origin.

As shown in Appendix C, the leading-order exponential behaviour of the remaining two solutions is

(4.22) \begin{equation} A_1 \sim {\rm{e}}^{4{\rm{i}}(-\mu )^{1/4}\xi ^{1/4}}, \quad A_2 \sim {\rm{e}}^{4(-\mu )^{1/4}\xi ^{1/4}}. \end{equation}

We show, by matching to the intermediate solution, that only $A_1$ has the correct behaviour at infinity, and thus is the only solution to be considered. Its asymptotic behaviour at the origin can be evaluated very efficiently using Barnes-type integral representations, as detailed in Appendix C, resulting in

(4.23) \begin{equation} A_1 = -\frac {\ln (-\mu \xi )}{2} + \frac {5}{4}-2\gamma + O(\xi ,\xi \ln \xi ), \end{equation}

at leading order; $\gamma$ is Euler’s constant.

4.3. The intermediate solution

We are still missing an intermediate solution that connects the solution (4.19) with $A$ given by (C7) to the pulse solution (4.9) with $\phi$ given by (4.11). In the process, we hope to find the parameters $\alpha ,\mu$ of the contact-line solution. On the one hand, the tail of the pulse solution (4.18) has a space dependence $\displaystyle {\rm{e}}^{(1-x)^{4/3}}$ , while, on the other hand, (4.22) implies that $\displaystyle {\rm{e}}^{x^{1/4}}$ , which clearly is not the same.

The required solution, which fits both asymptotic behaviours, is given by

(4.24) \begin{equation} \Delta \kappa = {\rm{e}}^{S}, \quad S = -\frac {G(x)}{t^{1/3}}, \end{equation}

as we now show; the exponent $t^{1/3}$ in the denominator is motivated by (4.18). Taking $S$ to be of order $t^{-1/3}$ and $x$ to be of order one, the leading order expression is

(4.25) \begin{equation} S_t + x^3 S_x^4 = 0, \end{equation}

and so the equation for $G$ becomes

(4.26) \begin{equation} \frac {G}{3} + x^3 G_x^4 = 0, \end{equation}

whose solution with boundary condition $G(1) = 1$ is

(4.27) \begin{equation} G = 3 (-1)^{1/3} \left (1 - x^{1/4}\right )^{4/3}, \end{equation}

and the roots are taken appropriately.

Near $x=1$ , we put $x=1-s$ , so that, for $s\ll 1$

(4.28) \begin{align} G \approx 3 (-1)^{1/3} \left (\frac {s}{4}\right )^{4/3}, \end{align}

which has to match the leading order behaviour (4.13), which implies that

(4.29) \begin{align} G \approx \frac {3}{2 \times 4^{4/3}}\left (1\pm \sqrt {3}{\rm{i}}\right ) s^{4/3}. \end{align}

Clearly, if the roots $(-1)^{1/3} = (1 \pm \sqrt {3}{\rm{i}})/2$ are selected, this is an exact match, so that now the solution in the intermediate region becomes

(4.30) \begin{equation} G = \frac {3}{2} \left (1\pm \sqrt {3}{\rm{i}}\right )\left (1 - x^{1/4}\right )^{4/3}. \end{equation}

This result is tested in figure 8, by plotting $\Delta \kappa$ multiplied by $t^{1/3}$ . For simplicity, we disregard the oscillations, by plotting maxima of $\Delta \kappa$ only, as function of the maximum position. Thus, the real part of (4.30) is plotted as the dashed line, which agrees progressively for earlier times.

Figure 8.

A solution of (4.1) at $\ln t = -16,-14,-12,-10,-8$ . Plotted are the maxima of $\Delta \kappa$ as function of their locations. For earlier times, the profiles converge toward the amplitude of the intermediate solution (4.30), which is $3(1-x^{1/4})^{4/3}$ .

Now for $x\rightarrow 0$ , (4.30) must match the large- $\xi$ behaviour of (4.22). The former limit leads to

(4.31) \begin{align} \ln \Delta \kappa \sim -\frac {3}{2}\left (1\pm \sqrt {3}{\rm{i}}\right ) \frac {1-4x^{1/4}/3}{t^{1/3}}, \end{align}

while the latter, using (4.22), gives

(4.32) \begin{align} \ln \Delta \kappa \sim -\frac {\mu }{\alpha t^{\alpha }} + 4{\rm{i}}(-\mu )^{1/4} \frac {x^{1/4}}{t^{(1+\alpha )/4}} . \end{align}

It is straightforward to confirm that if we identify

(4.33) \begin{equation} \alpha = \frac {1}{3}, \quad \mu = \frac {1 + \sqrt {3}{\rm{i}}}{2}, \end{equation}

then the two expressions become identical, having used that $(-\mu )^{1/4} = (\sqrt {3} -{\rm{i}})/2$ . To avoid confusion, we have chosen the plus sign for $\mu$ . However, since (4.21) is a real equation, replacing $\mu$ by its complex conjugate yields another solution. However, while $A_1$ grows at the exponential rate ${\rm{e}}^{2\xi ^{1/4}}$ , $A_2$ grows at the faster rate ${\rm{e}}^{2\sqrt {3}\xi ^{1/4}}$ , and therefore does not match.

To summarise, the contact line solution is

(4.34) \begin{equation} \Delta \kappa = {\rm{e}}^{-\frac {3(1 + \sqrt {3}{\rm{i}})}{2t^{1/3}}} A_1(\xi ), \quad \xi = \frac {x}{t^{4/3}}, \end{equation}

with $A_1 = \textrm {MeijerG}([[], []], [[0,0],[-2,-1]], -\mu \xi )$ , together with its complex conjugate. The notation for the MeijerG function is defined in Appendix C. Equation (4.21) is valid under the assumption that $t^{1/3}\xi A_{\xi } \ll A$ . Using (4.22), this implies that $\xi \ll t^{-4/3}$ , or $x$ is small. On the other hand, the validity of the intermediate solution relies on $G'^4/t^{4/3}$ being small in comparison with $G'^2 G''/t$ . In other words, $x$ must be greater than $t^{4/3}$ , so clearly there is an overlap for $t\rightarrow 0$ .

The physical (real-valued) solution based on (4.34) must be a linear superposition of real and imaginary parts. Thus, putting $A_1 = A_r + {\rm{i}} A_i$ , we obtain

(4.35) \begin{align} \Delta \kappa &= {\rm{e}}^{-3/(2 t^{1/3})}\left [\epsilon _1 \left (\cos \left (\frac {3\sqrt {3}}{2 t^{1/3}}\right ) A_r + \sin \left (\frac {3\sqrt {3}}{2 t^{1/3}}\right ) A_i \right )\nonumber\right. \\&\quad \left. + \epsilon _2\left (\cos \left (\frac {3\sqrt {3}}{2 t^{1/3}}\right ) A_i - \sin \left (\frac {3\sqrt {3}}{2 t^{1/3}}\right ) A_r \right )\right ]. \end{align}

Now with (4.33), $\ln (-\mu ) = -2\pi {\rm{i}}/3$ , and so, for small $\xi$ ,

(4.36) \begin{align} A_r = -\frac {\ln x}{2} + \frac {5}{4}-2\gamma + \frac {2}{3}\ln t, \quad A_i = \frac {\pi }{3}. \end{align}

This yields

(4.37) \begin{equation} \dot {\theta } = {\rm{e}}^{-3/(2 t^{1/3})}\left [\epsilon _1 \cos \left (\frac {3\sqrt {3}}{2 t^{1/3}}\right ) - \epsilon _2\sin \left (\frac {3\sqrt {3}}{2 t^{1/3}}\right )\right ], \end{equation}

where $\epsilon _1$ and $\epsilon _2$ remain to be determined. In order to do that, we have to go to next order in the expansion. The reason is that prefactors can be interpreted as logarithmic corrections in the exponential, but which are subdominant, and require a higher-order result.

4.4. The next order

To find the remaining amplitudes $\epsilon _1$ and $\epsilon _2$ , we need to continue to the next order. So far, we have accounted for the exponential terms, but we want to capture the terms from (4.16) and (4.22) that only grow algebraically. In the spirit of (4.24), we now try the generalised ansatz (with $G$ given by (4.30))

(4.38) \begin{equation} S = G(x) t^{-1/3} + G_1(x) + g_1 \ln t + g_2 \ln x + g_3 \ln (1-x). \end{equation}

The constants $g_1,g_2$ and $g_3$ are adjusted to match the required power laws.

Namely, the large- $\xi$ expansion of $A_1$ (cf. (4.22)) amounts to

(4.39) \begin{equation} S_{in} = 2\left (1 \pm \sqrt {3}{\rm{i}}\right ) \frac {x^{1/4}}{t^{1/3}} -\frac {9}{8}\ln x + \frac {3}{2}\ln t + \rm{const.}, \end{equation}

while, for small $s = 1 - x$ , the pulse solution yields

(4.40) \begin{equation} S_{out} = -\frac {3\left (1 \pm \sqrt {3}{\rm{i}}\right )}{2\times 4^{4/3}} \frac {s^{4/3}}{t^{1/3}} - \frac {4}{3}\ln s + \frac {5}{6}\ln t. \end{equation}

Comparing with (4.38), and matching to $\ln t$ and $\ln s$ in the expansion for small $s$ , we find that $g_1 = 5/6$ and $g_3 = -4/3$ ; matching to $\ln x$ in the expansion for small $x$ , we find that $g_2 = -9/8$ . Notice that there is a mismatch in the contribution from $\ln t$ , so (4.19) needs to be multiplied by $t^{-2/3}$ , so that $3/2-2/3 = \alpha$ as required. This change does not alter the leading-order equation (4.21), but changes only the next order.

To make everything consistent, (4.25) has to supplemented with terms of order $t^{-1}$ , leading to

(4.41) \begin{equation} S_t + x^3\left ( S_x^4 + 6 S_x^2 S_{xx}\right ) + 9 x^2 S_x^3 = 0. \end{equation}

Inserting (4.38) into (4.41), at order $t^{-1}$ gives

(4.42) \begin{align} 4(1 + G_1^{\prime})x^{5/4} - G_1^{\prime}x^{9/4} + x(x - 1)G_1^{\prime} - \frac {4x}{3} + \frac {x^{1/4}}{3} = 0, \end{align}

which can be integrated to give

(4.43) \begin{align} G_1 = \frac {4}{3}\ln \left (1 + x^{1/4} + x^{1/2} + x^{3/4}\right ) + \rm{const.} \end{align}

This finally leads to

(4.44) \begin{equation} \Delta \kappa = A_{med} {\rm{e}}^{-\frac {3}{2}\left (1 \pm \sqrt {3}{\rm{i}}\right ) \frac {(1-x^{1/4})^{4/3}}{t^{1/3}}} \frac {1 + x^{1/4} + x^{1/2} + x^{3/4}} {x^{9/8}(1-x)^{4/3}} t^{5/6} \end{equation}

in the intermediate region, where $A_{med}$ is an amplitude to be determined. Recall that (4.44) solves (4.41) and matches inner and outer regions (4.39) and (4.40), respectively, once a correction $-2\ln t/3$ has been added to the inner solution. This means that, to leading order, the inner solution is now $\Delta \kappa = t^{-1/2} {\rm{e}}^{-\frac {\mu }{\alpha t^{\alpha }}} A(\xi )$ instead of (4.19). This will generate an additional term in (4.21), which is proportional to $t^{1/3}$ . This motivates the ansatz to include a term of next order into the inner solution: i.e.

(4.45) \begin{equation} \Delta \kappa = {\rm{e}}^{-3\mu _{\pm }t^{-1/3}}\left [t^{-2/3} A(\xi ) + t^{-1} B(\xi )\right ]. \end{equation}

At leading order, this yields (4.21), as before; at the next order,

(4.46) \begin{equation} -\frac {2 A}{3} - \frac {4\xi A'}{3} + \mu B + \left (\xi ^3 B'\right )^{\prime \prime \prime} = 0, \end{equation}

which needs to be solved for $B$ , if the next order is required.

4.5. Matching

Finally we want to calculate the missing coefficients in the expression (4.37) for $\dot {\theta }$ , also including the correction $t^{-2/3}$ from (4.45). We do that by matching the missing amplitudes and phases of the successive regions, starting from (4.17). First, the limit of the intermediate solution (4.44) for $s\rightarrow 0$ (toward the right-hand end) is

(4.47) \begin{align} \Delta \kappa = 4 A_{med} {\rm{e}}^{-\frac {3}{2}\left (1 \pm \sqrt {3}{\rm{i}}\right ) (\zeta /4)^{4/3}} s^{-4/3} t^{5/6}, \quad \zeta = \frac {1-x}{t^{1/4}}. \end{align}

Comparison with (4.16) yields $A_{med} = A_p / 4$ , and, including the phase factor, for the real version of the intermediate solution, we have

(4.48) \begin{equation} \Delta \kappa = \frac {A_p}{4} \frac {1 + x^{1/4} + x^{1/2} + x^{3/4}}{x^{9/8}(1-x)^{4/3}} t^{5/6} {\rm{e}}^{-\frac {3}{2}\frac {(1-x^{1/4})^{4/3}}{t^{1/3}}} \cos \left (\frac {3\sqrt {3}}{2 t^{1/3}}\left (1-x^{1/4}\right )^{4/3} - \phi _p\right ). \end{equation}

Thus, the expression for the overlap region between (4.48) and (4.17) is

(4.49) \begin{equation} \Delta \kappa = \frac {A_p t^{5/6}}{s^{4/3}} {\rm{e}}^{-\frac {3\zeta ^{4/3}}{2\cdot 4^{4/3}}} \cos \left (\frac {3\sqrt {3}\zeta ^{4/3}}{2\times 4^{4/3}}-\phi _p\right ). \end{equation}

In the opposite limit $x\rightarrow 0$ (near the contact line), the intermediate solution (4.48) becomes (note that $f(x)\approx x^{-9/8}$ )

(4.50) \begin{equation} \Delta \kappa = \frac {A_p}{4} x^{-9/8} t^{5/6} {\rm{e}}^{-\frac {3}{2 t^{1/3}}(1-4x^{1/4}/3)} \cos \left (\frac {3\sqrt {3}}{2t^{1/3}}\left (1-\frac {4}{3}x^{1/4}\right ) - \phi _p\right ). \end{equation}{}

Now (4.50) can be compared with the large- $\xi$ limit of (4.35), but including the factor $t^{-2/3}$ implied by the leading-order contribution to (4.35). This shows that $\epsilon _1 = \epsilon \cos \phi _p$ and $\epsilon _2 = -\epsilon \sin \phi _p$ , with $\epsilon = -\sqrt {\pi } A_p t^{-2/3}/\sqrt {2}$ . Inserting this into (4.35), including the factor of $t^{-2/3}$ , we arrive at the complete solution for the contact line region,

(4.51) \begin{equation} \Delta \kappa = -\sqrt {\frac {\pi }{2}}\frac {A_p}{t^{2/3}} {\rm{e}}^{-3/(2 t^{1/3})}\left [\cos \left (\frac {3\sqrt {3}}{2 t^{1/3}}-\phi _p\right ) A_r + \sin \left (\frac {3\sqrt {3}}{2 t^{1/3}}-\phi _p\right ) A_i\right ], \end{equation}

where

(4.52) \begin{align} A_r + {\rm{i}} A_i = \textrm {MeijerG}\left ([[], []], [[0,0],[-2,-1]], -\frac {(1+\sqrt {3}{\rm{i}})\xi }{2}\right ), \quad \xi = \frac {x}{t^{4/3}}. \end{align}

Figure 9.

The change in contact angle $\dot {\theta }$ as obtained from a numerical simulation of (4.1) (solid line), compared with the prediction (4.54); a best fit (symbols) gives $\displaystyle t^{2/3}{\rm{e}}^{3/(2t^{1/3})}\dot {\theta } = -42.6 \cos ({3\sqrt {3}}/{2 t^{1/3}} - 0.4625)$ .

Now we can deduce the change in contact angle by taking the small- $\xi$ limit of (4.51), keeping the logarithmically diverging terms only, to obtain

(4.53) \begin{align} \Delta \kappa \approx \frac {A_p\sqrt {\pi }}{2\sqrt {2}} t^{-2/3} {\rm{e}}^{-3/(2 t^{1/3})} \cos \left (\frac {3\sqrt {3}}{2 t^{1/3}} - \phi _p\right )\ln x \equiv 2 b \ln x. \end{align}

According to (3.5), $\dot {\theta } = -4 b$ , and thus

(4.54) \begin{equation} \dot {\theta } = -\frac {A_p\sqrt {\pi }}{\sqrt {2}} t^{-2/3} {\rm{e}}^{-3/(2 t^{1/3})} \cos \left (\frac {3\sqrt {3}}{2 t^{1/3}} - \phi _p\right ), \end{equation}

which is the desired final dynamical equation for the contact angle. In figure. 9, the functional form of (4.54) is compared with numerical solutions of the linearised problem (4.1), from which $\dot {\theta }$ is obtained; almost perfect agreement is found. Using that $A_p = 51.39$ as found from the pulse solution, the predicted prefactor is $64.4$ , while a prefactor of $42.6$ is found numerically. The slight disagreement comes from the fact that convergence is slow, so $t \sim 10^{-4}$ is not yet asymptotic. Accessing earlier times numerically is difficult, since this would require an even greater range of spatial scales to be resolved with great accuracy. The observed phase factor of $0.4625$ is also in reasonable agreement with the prediction of $\phi _p = 0.561$ , once more found from the pulse solution.

5.1. Spectral gap

Again, the drop is considered to be symmetrical, and boundary conditions imposed at $x=0$ also apply at $x=2$ . Using (2.1), we supplement the eigenvalue equation (5.2) with the boundary conditions

(5.3) \begin{equation} G(0)=G(2) = 0,\ \lim _{x\rightarrow 0}h_{eq}^{3}G_{xxx}= \lim _{x\rightarrow 2}h_{eq}^{3}G_{xxx}=0. \end{equation}

Integrating (5.2) from $0$ to $2$ , the second condition guarantees that $\int \limits _{0}^{2}G(x){\rm d}x=0$ , ensuring conservation of mass. To transform the eigenvalue problem into one that is self-adjoint, we put $u=G_{x}$ , so that (5.2) becomes

(5.4) \begin{equation} \lambda u=(h_{eq}^{3}u_{xx})_{xx} \equiv A u, \end{equation}

defining an operator $A$ . Using (5.2), the boundary conditions (5.3) become

(5.5) \begin{equation} \lim _{x\rightarrow 0}(h_{eq}^{3}u_{xx})_{x} = 0,\ \lim _{x\rightarrow 0}h_{eq}^{3}u_{xx}=0, \end{equation}

and correspondingly at $x=2$ .

Integrating by parts and using (5.5), we find

(5.6) \begin{equation} \int \limits _{0}^{2}v(h_{eq}^{3}u_{xx})_{xx}{\rm d}x=\int \limits _{0}^{2}v_{xx}(h_{eq}^{3}u_{xx}){\rm d}x=\int \limits _{0}^{2} \left (h_{eq}^{3}v_{xx}\right )_{xx} u{\rm d}x, \end{equation}

as required for $A$ to be self-adjoint. As explained in more detail in Appendix D, we can conclude that the spectrum $\{\lambda _i\}$ of $A$ is real and discrete. It also holds (see Appendix D), that

(5.7) \begin{align} \int \limits _{0}^{2}h_{eq}^{3}u_{xx}^{2}{\rm d}x \geqslant C\int \limits _{0}^{2}u^{2}{\rm d}x \end{align}

for some positive constant $C$ . But since $A$ is self-adjoint, we also have

(5.8) \begin{align} \lambda \int \limits _{0}^{2}u^{2}{\rm d}x = \int \limits _{0}^{2}h_{eq}^{3}u_{xx}^{2}{\rm d}x, \end{align}

and thus $\lambda _i \geqslant C$ . This implies that relaxation towards equilibrium is exponential, at a strictly positive rate greater than $C$ .

Notice that the above arguments do not depend on the particular form of $h_{eq}$ , or the type of boundary condition, provided the velocity of the contact line is equal to zero and the domain is bounded, so that the volume is conserved. In these cases, we would still expect a strictly positive discrete spectrum. If, however, the physical domain is infinite, the problem will be associated with a non-compact operator and the spectrum is no longer expected to be discrete and possibly not strictly positive. To actually calculate the lowest $\lambda _i$ , we turn to a numerical solution.

5.2. Numerical solution

We need to find a solution to (5.2) that satisfies the boundary condition $G(0) = 0$ , once more placing the contact line at $x=0$ . Solving (5.2) by shooting, we use the initial conditions $G(1) = 1$ (a normalisation), $G'(1) = G'''(1) = 0$ (symmetry) and $G''(1) = \psi$ , where $\psi$ is a shooting parameter. For a generic value of $\psi$ , $G$ behaves according to

(5.9) \begin{equation} G \approx d\ln x + c - {\rm d}\lambda x\ln ^2x/2 + (d-c)\lambda x \ln (x) + \dots , \end{equation}

near $x = 0$ , as found by expanding $G$ in $x$ and $\ln x$ .

Figure 10.

(a) The constant $c$ , as defined in (5.9), as a function of $\lambda$ . Zeroes determine the eigenvalues – the first two are shown, with $\lambda _1 = 8.19$ and $\lambda _2 = 38.4$ , based on $V=2/3$ . (b) The first eigenfunction $G_1(x)$ , corresponding to $\lambda _1$ .

To find the eigenfunction, we first adjust $\psi$ such that the term ${\rm d}\ln x$ disappears, satisfying the flux condition (second equation of (5.3)). Since $G' = d/x + O(\ln ^2x)$ , this is achieved by demanding that $G'x = {\rm d} x / x + O(x\ln (x)) \rightarrow 0$ for $x\rightarrow 0$ . Next, we plot $G(0)$ or $c$ as a function of $\lambda$ (see figure 10(a); the condition $c = G(0) = 0$ determines the eigenvalue, since this is the boundary condition that $G$ has to satisfy at the contact line. We have chosen $V=2/3$ ; the result for arbitrary volume is found from rescaling. In figure 10(a), the range of $\lambda$ shown includes the first two zeroes of $G(0)$ , corresponding to the first two eigenvalues, which are strictly positive, in agreement with the above analysis. In figure 10(b), we show the first eigenfunction $G_1(x)$ , corresponding to $\lambda _1 = 8.19$ . Clearly, the eigenvalue scales with $V^3$ , so that the decay exponent in figure 3 is

(5.10) \begin{align} \lambda = (3V/2)^3\lambda _1 = (6\bar {h}/\pi )^3\lambda _1 = 0.057. \end{align}

The resulting decay law $\dot {\theta } \propto {\rm{e}}^{-\lambda t}$ (with the prefactor adjusted) is shown as the dashed line in figure 3, and demonstrates excellent agreement with the fully nonlinear simulation of the original partial differential equation, shown as the solid line.

From (5.1) it follows that $\dot {\theta } = -\epsilon \lambda G'(0) {\rm{e}}^{-\lambda t}$ . Since $G_1'(0) \approx -14.05$ is finite, it follows that the contact angle is changing as $t\rightarrow \infty$ . Alternatively, $\dot {\theta }$ can be calculated from (2.3). With the coefficients $c,d$ in (5.9) having been made to vanish, the local expansion of $G_1$ becomes

(5.11) \begin{align} G_1(x) = a x + \frac {a\lambda _1}{4} x^2\ln x + \dots . \end{align}

As shown on the left of figure 10, the next eigenvalue $\lambda _2$ is significantly larger than $\lambda _1$ , so the first eigenvalue will dominate for times of order unity.