Stability and Bifurcation of Dynamic Contact Lines in Two Dimensions

The moving-contact line between a fluid, liquid and a solid is a ubiquitous phenomenon, and determining the maximum speed at which a liquid can wet/dewet a solid is a practically important problem. Using continuum models, previous studies have shown that the maximum speed of wetting/dewetting can be found by calculating steady solutions of the governing equations and locating the critical capillary number, $Ca_{\mathrm{crit}}$, above which no steady-state solution can be found. Below $Ca_{\mathrm{crit}}$, both stable and unstable steady-state solutions exist and if some appropriate measure of these solutions is plotted against $Ca$, a fold bifurcation appears where the stable and unstable branches meet. Interestingly, the significance of this bifurcation structure to the transient dynamics has yet to be explored. This article develops a computational model and uses ideas from dynamical systems theory to show the profound importance of the unstable solutions on the transient behaviour. By perturbing the stable state by the eigenmodes calculated from a linear stability analysis it is shown that the unstable branch is responsible for the eventual dynamical outcomes and that the system can become unstable when $CaCa_{\mathrm{crit}}$, we will show that the trajectories in phase space closely follow the unstable branch.


I. INTRODUCTION
Understanding the shape and evolution of the interface between a fluid, liquid and a solid substrate is a classic problem in fluid mechanics and yet a remarkable number of open questions still remain [1,39]. There are two fundamental cases: an advancing contact line, where a liquid phase advances and 'wets' the solid, see figure 1(a)-(c), and a receding contact line, where a liquid phase recedes and 'dewets' the solid, see figure 1(d)-(f). Both experimental and theoretical studies [3,42] have shown that there is a critical contact line speed relative to the solid, beyond which stability is lost and the system ceases to return to a steady state. In the case of an advancing contact line (see figure 1(c)) this instability is characterised by fluid entrainment (which in many practical cases is air entrainment) and in a receding contact line (see figure 1(f)) a thin liquid film is deposited on the solid. The principle aim of this article is to provide insight into this instability and understand the dynamics of the system near the critical speed.
The critical speed where the instability occurs is associated with a fold bifurcation in the steady solution structure, which divides the steady solutions between a stable branch and an unstable branch (figure 2(a) and see [28] for a detailed mathematical description). For parameter values 'beyond the fold' there are no (known) two-dimensional steady states and the system must becomes transient and/or three-dimensional. In our system the appropriate non-dimensional parameter associated with the speed of the solid is the capillary number, Ca (see next section for a precise definition). Whilst analysis of the unstable branch of solutions (which exists for parameter values 'below the fold') can reveal important information about transient behaviour, the focus of theoretical studies has been mainly to calculate and characterise only the stable steady solutions immediately up to the critical speed [see, for example 5,14,46,55]. However, [5] hypothesised that the set of unstable solutions represents what they termed 'effective dynamics', i.e. the system's time-dependent trajectory for Ca > Ca crit closely matches the unstable branch when plotted using appropriate measures. If so, the unstable branch is not just an insignificant consequence of the fold bifurcation but provides unique insight into the system dynamics. The influence and importance of unstable states in fluid dynamics systems has been investigated in many different contexts, including shear flow [11], droplets [17], finite air bubbles [16,26] and a slide-coating flow [8]. Indeed, as shown in figure 2(b), where the phase-plane is sketched for a generic system with an stable ('attractor') and weakly unstable ('saddle-node') state, the unstable state can act as a separator of dynamical outcomes; the stable manifold is a dividing 'line' and the unstable manifold connects to the stable state. In this article we adapt these dynamical systems ideas to the moving-contact-line problem to reveal the role of the unstable solutions. We calculate the bifurcation structure and stability properties of the steady solutions and relate these to time-dependent calculations in the sub-critical (Ca < Ca crit ) and super-critical (Ca > Ca crit ) regimes.
We now provide some important background on moving contact lines. It is well known that the classical 'moving contact line paradox', as described in [22] can be alleviated if there is slip near the contact point. If this slip occurs in an inner region, as considered by [56] and [9], then bending of the interface occurs in an intermediate region where viscous effects can cause the liquid-fluid interface to curve sharply. In this formulation, it is often assumed that the intermediate region connects to an outer region where the interface retains its static meniscus shape. The possible asymptotic matching of these regions has critical consequences and provides insight into the bifurcation structure of the steady solution space. In a series of remarkable articles, it was shown, by solving a lubrication model for a liquid-vacuum system, how the curvature of the inner and outer regions can be asymptotically matched. For the advancing contact line this can be achieved for all values of Ca, but for the receding contact line, the matching fails when Ca is past some critical threshold, interpreted as Ca crit [12][13][14]. The bifurcation structure of the stable and unstable branches of the receding contact line was then fully described using matched asymptotics and bifurcation theory by [5] for Ca 1, and Ca crit was determined to occur at a fold bifurcation. The aforementioned lubrication analysis has been extended to general liquid-fluid systems, where the viscosity of the fluid phase is considered non-zero [4,6,24] and also for the full Navier-Stokes equations [53][54][55]. A key result from these studies is that for the advancing contact line the presence of viscosity fundamentally alters the bifurcation structure and a fold bifurcation appears at a finite Ca. [55] showed that the fold bifurcation in the advancing contact line problem occurs when the horizontal air pressure gradient matches the strength of capillary-stress gradient near the contact point. It was also demonstrated that using the lubrication model, in both phases, poorly predicts Ca crit when compared to the full Navier-Stokes equations for the advancing contact line [53][54][55]. Other physical effects such as Maragoni flows, inertia and gravity, and shear thinning/thickening were also found to preserve the fold bifurcation [7,[30][31][32][33]55].
In physical terms, in the advancing case, the critical behaviour indicates the threshold at which fluid entrainment occurs where, experimentally, a three-dimensional saw-tooth pattern emerges as observed in a variety of different flow configurations e.g. liquid films [36], drop impact [35,52] and plate penetration in a liquid bath [19]. In the receding case, however, the fold bifurcation marks the onset of thin-film deposition [44,45]. Interestingly, despite the 3D structures of air entrainment [18,19], 2D models appear to accurately predict the transition point, an observation which is yet to be understood [see, for example 31,47,54]. Transversal three-dimensional perturbations have been considered for the receding contact line [43] and the advancing contact line [53], both using a lubrication model, but a stability analysis using the full hydrodynamics equations has not yet been conducted.
In this article we consider the bifurcation structure devoid of any limitations from the lubrication approximation, allowing us to naturally consider both advancing and receding cases simultaneously. Our analysis of two-phase contact-line stability will focus on steadystate solutions using a hybrid model; the liquid phase is modelled using the Navier-Stokes equations and the fluid phase is accurately modelled using a lubrication approximation [see 30-33, 47, 51].
The structure of the article is as follows. In § II we describe the hydrodynamic equations that describe the system. In § III we calculate the steady solution curves to determine the critical parameters associated with the loss of stability of the system. In addition, we perform a numerical linear stability analysis that reveals the significance of the unstable branch to the transient dynamics of the system. By treating the governing equations as a dynamical system we form a generalised eigenproblem that can be solved numerically to determine and quantify the stability of the solution branch. Next, in § IV, by solving a timedependent initial value problem (IVP) numerically we are able to demonstrate that far from having a passive role, the unstable branch represents, in the language of dynamical systems, the 'basin boundary of attraction' of the stable state. Furthermore, by examining the phaseplane of the solution trajectory, we discover that the subsequent unsteady time-evolution is intrinsically linked to the unstable branch and are able to confirm the prediction of [5] that the solution moves quasi-statically along the unstable branch. Viewing the trajectories FIG. 1. The moving contact line problem in a channel geometry in a frame of reference that moves with the liquid. Panels (a)-(c) describe the advancing contact line problem whilst (d)-(f) describe the receding contact line problem. In both cases as the speed of the substrate, U * , increases the system is first stable (panels (a) and (d)) before the system becomes unstable and air entrainment (panel (c)) or a thin-film formation (panel (f)) occurs.   2. (a) A sketch of a typical fold bifurcation structure. A solution measure is plotted against a control parameter to form a 'solution curve'. At a critical value, two branches -one stable (solid line) and one unstable (dashed line -) meet. The location of their intersection is known as a fold bifurcation. Beyond the critical value there are no (known) steady states. In our specific problem the control parameter is Ca and the solution measure is either the interface length or meniscus rise. (b) A generic two-dimensional phase plane for a parameter value less than the critical value with a stable state (an 'attractor' on the stable branch, see (a)) and an weakly unstable state (a 'saddle-node' on the unstable branch, see (a)). The unstable/stable manifolds of the unstable state are dashed/dotted respectively. through the lens of the phase plane will also allow us to understand if, and how, the system becomes unstable when Ca < Ca crit and also provide criteria that could potentially enable suppression of this instability. Finally, in § V, we discuss the implications of these results and some possible future research.

II. GOVERNING EQUATIONS
We now discuss the hydrodynamic model and the assumptions that allow us to derive an accurate simplified hybrid model that is used in the calculations thereafter. The following discussion applies to both the advancing and receding contact lines, although the demonstrative figures only show the advancing contact line.

A. Full Hydrodynamic Model
Motivated by the system used in [54], we consider two-dimensional flow between two parallel plates, as shown in figure 1. In the following discussion we denote dimensional quantities using a * superscript. Two fluids of viscosity µ * 1,2 , and density ρ * 1,2 , fill the channel bounded by two rigid plates which are separated by a fixed height H * , subscripts with 1 indicate the upper fluid (the fluid phase) and 2 indicates the lower fluid (the liquid phase). In our system the left plate moves with constant speed in the y−direction U * and the right plate is stationary. For a receding contact line U * > 0 and an advancing contact line U * < 0. The left wall is moving with a speed U * and the right wall is stationary. The coordinate system is centred on the contact point between the two fluids and the left (moving wall), see panel (a) of figure 3. The fluid flow of each phase is governed by the two-dimensional Navier-Stokes equations. All speeds, lengths, pressures and times are scaled by U * , H * , µ * 2 U * /H * and H * /U * respectively. Finally the viscosity ratio, denoted χ, is defined with respect to the liquid phases, i.e. χ = µ * 1 /µ * 2 . As in previous studies [30-33, 48-50, 54, 55] we apply the Stokes-flow approximation so that the Reynolds number, Re = U * H * ρ * 2 /µ * 2 , is negligible and assumed zero; results in [55] show Re can have an influence at sufficiently high values but it does not qualitatively alter their conclusions. We assume that gravitational effects are negligible throughout. The nondimensional computational domain is shown in figure 3(c). The line corresponding to the left plate is denoted Γ 1 , the right plate Γ 2 , the bottom boundary Γ 3 , the free-surface Γ 4 and the top boundary Γ 5 . The fluid and liquid domains are denoted by Ω 1 and Ω 2 respectively. The Stokes-flow equations, for the fluid velocity, u i = (u i , v i ), and pressure, p i , in each phase can be written as (1) On the left (moving) and right (stationary) walls, Γ 1 and Γ 2 respectively, we implement a Navier-slip condition written as where n and t are the vectors normal and tangential to each wall, U = (0, U ) is the nondimensional speed of the wall and λ is the non-dimensional slip length which, for simplicity, we assume to be the same in each phase (see, [47] for potential extensions). We choose to implement a Navier-slip condition on the stationary wall for consistency with the Navier-slip condition on the moving wall, although we could fix u i = 0 on Γ 2 and get similar results (see, for example [30,55]). The stress tensor in each phase τ i is defined as where I is the identity matrix and δ 1 = χ, δ 2 = 1. On the interface between the two fluids we assume a constant surface tension, γ * , so that the dynamic boundary condition can be written as n is the normal of the interface pointing towards the fluid phase, κ = ∇·n is the curvature of the interface and Ca = µ * 2 |U |/γ * is the capillary number. We denote the unknown position of the interface as r = (x s , y s ), see figure 3(b), so that the kinematic condition on the interface can be written as In addition, we have to specify the angle the interface makes on the lower and upper walls. These angles can be allowed to vary with the capillary number, slip-length or other quantities but we choose the simplest approach and choose constant values, i.e.
It is straightforward to replace the conditions in (9) with an equation involving Ca and other quantities, but this is not the focus of the article. Finally, we implement fully-developed flow conditions on the inflow and outflow boundaries, alongside a pressure drop across the domain so that, The full hydrodynamic system is defined in (1)-(12) with the following high-dimensional state vector (denoted w) of unknowns; It is worth noting that we model the effect of varying the speed of the wall by varying Ca and that the non-dimensional slip-length, λ, can be varied to investigate changes in physical channel width. Finally as we are interested in the contact line of the left plate, we set θ 2 = π/2 on the right plate, in all simulations, for simplicity. Therefore we have a set of control parameters Ca, λ, χ, θ 1 , p out , that need to be specified in order to solve (1)- (12).

B. Hybrid Model
The computational cost of the full model can be drastically reduced by solving the thinfilm equations where they are valid [23,34,38], leading to a hybrid model [see [30][31][32][33]51] which approximately halves the complexity of the problem, as unknowns in the fluid phase are only computed on the interface. The difference of our approach from previous implementations is that our hybrid model takes into account time-dependence so that stability can be probed and IVP calculations can be performed. The key assumption is that typical horizontal length scales are small compared to typical vertical length scales so that the horizontal component of velocity in the fluid phase is small, i.e. u 1 1. The full derivation is discussed in Appendix A and the computational domain is shown in figure 3(b).
The effect of this reduction in the fluid phase is to replace a full two-dimensional description, given in (1), by a one-dimensional equation for the fluid pressure, p 1 , on the interface only. This equation can be stated as where h is the horizontal distance, from the left hand solid, to the interface (see figure 3), Q 1 is the flux, and the constants A and B are functions of χ, λ and u 2 and are given in Appendix A. The fluid phase is coupled to the liquid phase through the applied traction given in (A6). We now have a system of PDE described by equations (2)- (12) and (15) with the high-dimensional state vector of unknowns: We validate the hybrid model by comparing to the full hydrodynamic model in Appendix A.
Finally we note that this approach is strictly only valid for the advancing contact line problem but, as shall be shown later, the receding contact-line problem is effectively a onephase problem, (c.f. figure 8) and implementing the hybrid model for a receding contact line does not significantly change the value of Ca crit (Appendix A).

System Parameters and Integral Measures
We now describe additional system parameters and measures that will be useful in computing and describing the steady and time-dependent solutions. The pressure at the outflow boundary, p out , is determined by an integral volume constraint acting on the liquid phase, i.e.
where V is the volume per unit length of the domain (corresponding to the area of the computational domain). In our numerical calculations the position of the contact points on the moving and stationary plates are both allowed to move so that (17) can be satisfied. For ease of presentation, we post-process and rescale the solution so that the origin is always at the contact point of the moving plate. In the calculations that follow we choose V = 5 which, after careful experimentation, is large enough for fully-developed flow to occur near the outflow boundary, Γ 1 . We find that the solutions are independent of values of V ≥ 5 that we choose. For a fixed set of parameter values, as defined in (14), we calculate steady solution curves by setting the time derivatives in the governing equations to zero and then solving the resulting steady system. As we vary Ca, and then subsequently calculate a solution, a solution curve will be traced and a fold bifurcation will occur at the critical value of the capillary number, denoted Ca crit . Whilst it is possible to trace a solution branch around the fold numerically by a pseudo-arclength continuation method [see, for example 10], we implement an alternative, bespoke approach. We expect the interface length, L, to increase monotonically as the curve is traced out around the fold and therefore is a suitable candidate for a continuation parameter that allows us to calculate solutions smoothly around the fold. To achieve this we let Ca become an unknown parameter that is determined by setting the total length of the interface, i.e.
This approach enables us to trace solution curves around the fold by incrementally increasing L and solving the system of equations with Ca determined by (18). We also emphasise that (18) can only be implemented in steady calculations as a means of tracing the solution curve whereas the volume constraint in (17) is applied in both steady and time-dependent calculations.
Finally, when describing the steady solutions and time-dependent solutions we use the meniscus rise (more specifically, the vertical distance between the two contact lines) defined as as a convenient solution measure (as previously considered, for example, in [24]; see figure 3 (in this paper) for reference).

C. Numerical Method
The governing equations are solved using the finite-element method from within the open-source oomph-lib object-orientated multi-physics library [20]. The structure and implementation of our equations follows that of [50]. Following multiplication of the equations by a test function, ψ, and then an integration over the domain, the boundary integrals that result from integration by parts require the traction to be specified on each of the boundaries. The dynamic condition, (6), and the Navier-slip condition, (4), therefore can be implemented as a natural condition by these boundary integrals.
Special care has to be taken at the contact point. In other studies [30-33, 54, 55] the contact angle is imposed as an essential boundary condition at the expense of solving a component of the momentum equations at the contact point. We adopt the approach of [50] and impose the contact angle as a natural boundary condition on both the intersection of the free-surface with the left plate (Γ 1 ) and the symmetry plate (Γ 2 ). We therefore introduce a field of Lagrange multiplier unknowns on Γ 1 and Γ 2 which are determined from the weak form of the no-penetration condition. We refer the reader to [50] for a detailed description of this implementation (we adopt approach (B) in their nomenclature).
As is standard, the fluid velocities are interpolated using bi-quadratic shape functions and the pressure using linear continuous shape functions with Taylor-Hood triangular elements. We choose to mesh the liquid domain using an unstructured triangular grid; see figures 3(b) and (c). The mesh is considered to be a fictitious pseudo-solid with the position of the nodes coming as part of the solution. The weak form of the kinematic condition, (7), is imposed as an essential condition and determines a field of Lagrange multipliers (not to be confused with the Lagrange multipliers in the previous paragraph) that act on the solid deformation equations which in turn determines the shape of the unknown interface, r, see [37] for more details. We note that this approach results in a large system of equations which is disadvantageous, but it also allows for the interface to become highly deformed and even multi-valued, as well as naturally handling unsteady flow (where the domain could significantly change shape c.f. Figure 1), which is difficult to achieve if the mesh is structured.
To solve the hybrid equation, (15), it is convenient to introduce two fields of unknowns on the fluid interface, the pressure p 1 and flux Q 1 , interpolated using quadratic shape functions.
We solve two equations in their weak form: and S ∂h ∂t Equation (20) projects the flux from the lubrication equation onto the finite element space and then (21) ensures mass is conserved in the fluid phase.
The resulting discretised equations are solved using Newton's method using the SuperLu numerical algebra package [29]. For time-dependent calculations the solution is updated in time using a backwards-difference second-order Euler method (BDF2).
Around the contact line the interface becomes highly deformed due to viscous bending and the pressure and velocity gradients are large, see figures 3(d) and (e). In steady calculations, as Ca → Ca crit , we expect the number of elements required in the vicinity of the contact line to increase to ensure a smooth converged solution. We re-mesh the domain according to a ZZ error estimator [57], which measures the discontinuity of strain rate gradients between adjacent elements and interprets this as a measure for the local error. Typically we set a minimum error as 10 −6 and a maximum error 10 −3 , so that elements with error above this range get refined and those with error below this range get unrefined. We allow element sizes from 10 −12 to 10 −2 to accommodate these error estimates. We do not adapt the mesh at each calculation; rather we adapt the mesh based on the condition that where θ c is the computed angle based on (x 1 , y 1 ) and (x 2 , y 2 ), the coordinates of the nodes on Γ 4 directly at the contact line and immediately adjacent, respectively. The number of elements and their sizes are highly dependent on λ and Ca. As an illustrative example, for steady solutions at λ = 0.1, χ = 0.1 and V = 5, the resulting mesh has ∼ 10 3 triangular elements and ∼ 10 5 discretised unknowns at Ca = Ca crit .

III. LINEAR STABILITY ANALYSIS
We now present the stability algorithm and results. Rather than perform a standard normal modes reduction to the Orr-Somerfeld equations -[see, for example 40] we take a more general approach that determines the modes as part of the solution. The analysis below is independent of the model, and although the results we present are from the hybrid model, these results also follow from the full model.
In both cases the PDE system can be written as where R is a nonlinear operator and w(t) represents a state of the system at time t, given as a vector of all the unknowns (either (13) or (16)). The time derivatives,ẇ, appear in linear combinations in our system so we can decompose R into a linear mass operator, M, that operates on the time-derivatives in the problem and a nonlinear operator, F, that operates on the spatial derivatives in the problem so that (23) becomes To proceed we write the state of the system, w, as a Taylor expansion, i.e.
where w is a base state only dependent on spatial variables, ε 1 is a small parameter, g is an eigenmode that is dependent on spatial variables only and σ is the growth rate of the perturbation. The expansion in (25) represents a general class of perturbations that satisfy the boundary conditions of the problem and are in-plane perturbations; we are not extending to the third dimension, a problem we will discuss later.
Substituting (25) into (24) gives a series of problems that have to be solved at each order of ε. At leading order we have The solution, w , is the steady state of the system. At first order we solve where J (w ) is the functional derivative of the nonlinear operator F applied at the steady state w = w . Equation (27) is a generalised eigenvalue problem that can be solved to find g and σ. The eigenspectrum of σ determines the stability of the steady solutions. If at least one of the spectrum of σ has a positive real part then the steady state is linearly unstable.
Conversely if the entire spectrum lies in the left-half of the complex plane then the solution is stable. In general there will be an infinite number of these eigenmodes and thus we can write the linearised solution as where c.c denotes the complex conjugate and a n are arbitrary constants set by the initial conditions. When the system becomes discretised the operators M and J are represented by the mass-matrix and Jacobian matrix, respectively. The mass matrix representation of M is highly rank deficient as the only time-derivatives occur at the fluid-liquid interface and special care has to be taken to ensure that the solution to (27) has converged. We use the Anazasi linear algebra library which is an iterative eigensolver that can solve highly rank-deficient eigenproblems [21]. As the spectrum has an infinite number eigenvalues, the discretised spectrum will have a finite number of eigenvalues, proportional to the number of unknowns in the problem. We find a small subset of eigenvalues which have the largest real part as these will be the modes visible in the transient dynamics; large negative eigenvalues correspond to eigenmodes that decay very rapidly. We validate the calculations using a simplified lubrication model and present this in Appendix B.

A. Stability of the Solution Branches
We now discuss the bifurcation structure and the corresponding stability results of the advancing and receding dynamic contact line problems. Figures 4 and 5 show the bifurcation structures in a typical advancing case (χ = 0.1, λ = 0.1) and receding case (χ = 0.0, λ = 0.1) respectively. Notably, our focus here is on providing insight into the stability structure, rather than necessarily probing the precise values from experimental analyses, where the slip  R 3 is the solution where the inflection point on the interface first becomes stationary, i.e. when dx/dy = 0. R 4 is the solution just before numerical calculations cease to converge which occurs when the interface becomes sufficiently deformed so that it touches the stationary right plate Sta ble Br anc h U n s t a b l e B r a n c h Un sta ble Br an ch length could be far smaller and therefore typically require more computational resources. Previous works [54] have shown that whilst changes in slip length can have a weak effect on Ca crit , they do not qualitatively alter the physical mechanisms at play (similarly for smaller viscosity ratios, e.g. with a glycerol-air system).
The solution curves are shown in the (Ca, Y ) projection of the solution space (see (19) for a definition of Y ). The markers on the curve indicate specific solutions which are shown in the inset panels labelled A 1 -A 4 for the advancing contact line, and R 1 -R 4 for the receding contact line. The eigenspectra of A 1 ,A 2 and R 1 ,R 2 , and at the fold are shown in inset panels for the advancing and receding contact lines respectively. In both the advancing and receding cases, as Y increases, the solution curve experiences a fold which separates the lower branch and upper branch. The eigenspectra is real and at the fold a single eigenvalue crosses the imaginary axis, as expected. The eigenspectra also indicates that the A 1 /R 1 states are 'attractors' of the system and A 2 /R 2 states are weakly unstable 'saddle-nodes' (figure 2(b)), thus numerically confirming the lower branch is stable (solid curve) and the upper branch is unstable (dashed curve).
The interface has an inflection point near the contact point. We measure the angle at the interface inflection point (to the downwards vertical) and define this as θ app ; the apparent contact angle [33,54]. Notably, as can be seen from solutions A 1 ,A 2 and R 1 ,R 2 in figures 4 and 5, θ app < 180 • not only on the stable branch, but also immediately after the fold on the unstable branch. Further up the unstable branch, as can be seen from solutions A 3 ,A 4 and R 3 ,R 4 in figures 4 and 5, θ app ≥ 180 • and the interface becomes multi-valued (as a function of h). In the advancing/receding cases the solution curve terminates when the interface is sufficiently deformed so that the interface touches the right/left plate, respectively, so that the interface effectively 'pinches' off the fluid domain, as seen in solutions A 4 and R 4 in figures 4 and 5.
The steady-solution curves of the advancing and receding contact line problems, although treated separately in figures 4 and 5, are actually two halves of the same solution space. Figure 6 shows the connection for λ = 0.01, χ = 0.1, θ m = π/2 where the signed meniscus rise, y(1) − y(0), is plotted against Ca × U , where U = ±1 with +/− corresponding to the receding/advancing problem respectively. The advancing and receding curves occupy the second and fourth quadrants in this projection and the location of the respective folds in each quadrant highlights that the receding contact line becomes unstable before the advancing contact line [6]. Finally, we note that in both cases the solution curve does not experience additional bifurcations as Y increases along the unstable branch. For a system where gravity is included it is known that within the lubrication approximation, the solution curve (for the receding contact line, at least) oscillates around a fixed value of Ca = Ca * , [see 5], experiencing multiple saddle-node bifurcations as Y → ∞. Preliminary calculations show that, if gravity is included, the oscillations are also present in the advancing/receding hybrid system, although, for brevity, we do not show the results here.

B. Physical Interpretation of the Bifurcation
As discussed in [55] for the advancing contact line, the fold occurs when the fluid pressure gradients (fluid 1) are comparable to the capillary stress gradients [see 53] near the contact line, i.e. when dp 1 ds This is because as Ca → Ca crit the air pressure gradients near the contact-line will increase as the system seeks to 'pump' air out of the region near the contact line to maintain system state. Eventually these air pressure gradients will exceed the capillary stress gradient and the system will be unable to maintain a stable steady equilibrium. Figure 7 shows the evolution of the quantities on either side of (29) calculated at the inflection point for the advancing case. The pressure and capillary stress gradients balance close to Ca crit , as seen by the intersection of the curves, which confirms the ideas of [55]. For a given χ the bifurcation structure is shown in the inset. We note the qualitative behaviour shown in the inset is independent of the value of λ.

C. Fold-Tracking
We can take advantage of the fact that at the fold bifurcation the leading eigenvalue crosses the imaginary axis to develop an algorithm for finding Ca crit . We augment the system with the additional constraint and let another control parameter come as part of the solution. It is convenient to let the interface length, L, be determined by (30) so we are able to track the evolution of the Ca crit as another parameter, the viscosity ratio χ, for example, is varied. This is a robust way of tracking the fold without having to recalculate the solution curve for every set of parameters, as previously considered in [24] and [54]. If we vary χ and calculate Ca crit we observe that the curve of the loci of Ca crit does not itself experience any bifurcation, (co-dimension 2 bifurcations), as seen in figures 8(a) and (b). In addition we observe that that the bifurcation structure also remains intact when the slip-length, λ, is varied as the different coloured curves in figures 8(a) and (b) indicate. Therefore we expect the dynamics to be qualitatively similar (from a dynamical systems perspective) regardless of the viscosity ratio or slip-length.
An important observation is that the advancing and receding cases differ significantly as χ → 0. For the advancing case, Ca crit → ∞ in this limit, whilst for the receding case it tends to a finite value. This indicates that the viscosity of the fluid phase has to be taken into account for the advancing contact line in order to describe the bifurcation structure. In contrast the receding contact line is essentially a one-phase problem, particularly if the fluid is a gas and qualitative features of the bifurcation structure are the same regardless of the viscosity of the fluid.

D. Eigenmode Perturbations
We now discuss the nature of the eigenmodes resulting from the stability analysis. The modes corresponding to the three leading eigenvalues of the unstable branch are shown in figure 9. These eigenmodes correspond to the base state w = A 2 , R 2 in figures 4 and 5, respectively. In this figure the dotted profile indicates the steady interface shape and the coloured lines indicate the shape of the interface when it is perturbed by a single eigenmode, i.e. w = w ± ρg i , i = 1, 2, 3.
The dashed/solid curves correspond to the +/− sign, respectively. The amplitude of the perturbation, ρ, is constrained so that the meniscus rise of the perturbation is no more than 10% of the rise of the steady solution. Each successive mode intersects the steady interface at precisely one more location in similarity to the form of the eigenmodes in a related lubrication model (see figure 22). Thus, the effect of adding higher-order eigenmodes to the steady state is add extra corrugations to the interface. Concentrating on the leading eigenmode alone, the action of adding a multiple of g 1 to a steady solution, i.e. taking i = 1 in (31), stretches/shrinks the interface according to the ± sign with no additional corrugations. Figure 10 shows the stable A 1 , and unstable A 2 states with solid lines and the dotted curves indicate the perturbed interface profiles from the A 1 state using the leading g 1 eigenmode. This figure demonstrates that if we can continuously 'stretch' the nonlinear stable state by increasing the strength of the perturbation, ρ, and can eventually achieve an interface profile similar to the unstable steady state, A 2 , which will have consequences, as discussed in the next section. We can also continuously deform the unstable branch in the same manner to match the interface of the stable branch. We shall denote perturbations using the leading eigenmode only in (31) as stretch perturbations. In a physical experiment, perturbations will naturally emerge from the presence of 'noise' in the system. This 'noise' will occur from random fluctuations of the interface and contactline, but rather than apply a stochastic perturbation to the system we will mimic this 'noise' by using more eigenmodes in the perturbation, i.e.
where for the purpose of simple illustration we choose the amplitude coefficients, ρ i , to be equal. Figure 11 shows the perturbed interface as the value of N increases from 1 to 10, which demonstrates that by increasing the value of N we are able to perturb the nonlinear steady interface with increasingly more corrugations or 'noise'. Henceforth we shall call perturbations of this form as corrugation perturbations. The two forms of perturbation discussed here, namely the 'stretch' and 'corrugation' perturbations, will be used in the next section to understand the subsequent time-dependent behaviour of the system after systematically applying a perturbation to a steady state.

IV. TRANSIENT DYNAMICS
Now that we have calculated the steady solution branches and quantified their stability, we attempt to answer the two questions of fundamental importance: 1. How do the steady-states, stable and unstable, and the resulting bifurcation structure help us understand the time-dependent behaviour of the system when Ca < Ca crit ?
2. What is the time-dependent behaviour of the system when we choose initial conditions (IC) beyond the fold, i.e. when Ca > Ca crit ?
To address these questions we solve the time-dependent hybrid PDE as an IVP. It is useful to define solution measures that will help visualise and aid the discussion. In our formulation Ca corresponds to the speed of the wall and not the speed of the contact point. It is therefore useful to introduce an 'effective' Ca, based on the contact line speeds relative to the wall (as discussed in [5]), which we denote Ca and is defined as where U = ±1 is the non-dimensional speed of the wall, the ± sign corresponds to the advancing and receding cases, respectively, and U cl = dy cl /dt is the speed of the contact line. We remark that for steady solutions U cl = 0, and hence Ca = Ca and the timedependent phase-plane trajectories can be directly compared to the bifurcation structure in the (Ca, Y ) plane. We also introduce a system measure to quantify the size of the perturbation. Let the meniscus rise of a steady solution be Y , where indicates the base solution the perturbation wil be measured against. We can then define a quantity that measures the deviation of the perturbation from the corresponding steady state: as demonstrated schematically in figure 12. If ∆(t, ) > 0 then the meniscus rise of that current state is larger than that of the steady interface indicated by and vice versa if ∆(t, ) < 0 (see panels (a) and (b) of figure 12 respectively).

A. Perturbations from a steady state IVP: Ca < Ca crit
We now consider the first question and look at the dynamics of the system when Ca < Ca crit . Our methodology is to start at either the stable or unstable state and perturb it using either a 'stretch' or 'corrugation' eigenmode expansion, which we will consider separately. We then run a series of IVPs to examine the transient behaviour and eventual dynamical outcome.

'Stretch' perturbations
For the 'stretch' perturbations we use the leading eigenmode only in the perturbation and set the IC to be w(t = 0) = w + ρg 1 .
Initially we concentrate on the advancing contact line case. Using the IC stated in (35), small perturbations from the A 1 stable state (figure 4) decay and the system (unsurprisingly) relaxes back to its stable configuration. Figure 13 demonstrates this by tracking the value of Y in time for two different perturbations, corresponding to ∆(0, A 1 ) < 0 and ∆(0, A 1 ) > 0, of the stable state near the fold (parameter values quoted in the caption). Furthermore, as seen in the insets on the right of figure 13, the decay rate excellently matches the value of the leading eigenvalue, σ 1 , obtained from the linear stability analysis. For the same parameter values, we can also perturb the A 2 unstable state (figure 4) by its leading eigenmode, which is unstable. This case is more interesting, and we see that if ∆(0, A 2 ) < 0 (i.e. the perturbation 'contracts' the A 2 interface) then the system returns to the stable state, whereas if ∆(0, A 2 ) > 0 (i.e. a 'stretch') then the contact line speed eventually diverges and the calculations fail to converge. Figure 14 shows the time-signal of Y showing these two different dynamical outcomes. As shown by the inset labelled 'Final State' the contact point velocity appears to diverge when θ app is close to 180 • . Similar outcomes occur for the receding contact-line problem. Perturbations, using the leading eigenmode, of the R 1 (figure 5) stable state relax back to the stable equilibrium (results not shown). Figure 15 shows the time evolution of perturbations from the R 2 (figure 5) unstable state. Panel (a) shows the time-signal of Y when ∆(0, R 2 ) < 0 and panel (b) when ∆(0, R 2 ) > 0. In the first case the system relaxes back to the R 1 state but if ∆(0, R 2 ) > 0 then, unlike the advancing contact line where the contact-line speed diverges, a thin-film develops that grows in size at a linear rate as t → ∞ as shown in panel (b).
These IVP calculations show that if we consider the class of perturbations representing 'stretches' using the leading eigenmode only, then the indicator of whether the system returns to the stable state is that meniscus rise of the initial perturbation is smaller than the meniscus rise of the unstable state, i.e. the condition for stability is In the language of dynamical systems the unstable state represents the 'boundary of the basin of attraction' of the stable state, when considering simple stretches of the stable interface corresponding to perturbations using the leading eigenmode. The unstable state is therefore not just a trivial consequence of the steady bifurcation structure but also has an important role in the underlying transient dynamics.
Another important observation is that the calculation in figure 15, for the receding contact line, provides evidence for the claim from [5] that the solution curve represents the effective dynamics of the system "in which the state of the solution moves quasi-statically along the solution curve". In the receding case, and, although to a lesser degree, the advancing case, the trajectory of the system in the phase-space (Ca, Y ) is qualitatively similar to the steady solution curve when plotted in the same diagram. The inset diagrams in figures 14 and 15 labelled 'Phase plane' show the steady solution curve and the trajectory of the system shown with an arrow. In both advancing and receding cases, when ∆(0, A 1 /R 1 ) < 0 the trajectory follows a similar path to the solution curve. In the receding case, when ∆(0, R 2 ) > 0 the trajectory closely follows the upper reaches of the unstable branch. The similarity between the steady solution curve and the time-dependent trajectories will be discussed in more detail in the next section.  efficiency. We set the initial condition of the system therefore as We concentrate solely on the advancing case and, as before, examine perturbations from the stable A 1 and unstable A 2 (figure 4) states using the IC prescribed in (38). Figure 16 shows the time-dependent results for incrementally increasing values of ρ in (38). Panel (a) shows trajectories in the (L, Y ) phase plane projection, where L is the total arclength of the interface. To help orient the trajectories we have added artificial axes which are centred on the unstable A 2 state.
After perturbing the A 1 state, the system experiences initial transient growth in the value of Y , corresponding to the weak attraction of the unstable A 2 state, before either the system settles to the stable state (solid trajectories) or the contact-line velocity diverges (dotted trajectories), indicated by the asterisks in panel (a). The dashed line between the A 1 and A 2 state is the response of the system if only the leading eigenmode is retained in the initial perturbation, as considered in the previous section on 'stretch' modes which can be interpreted as the approximate unstable manifold of the A 2 state.
The nonlinear trajectories that return to the stable state all eventually collapse on the unstable manifold, once the higher-order modes have sufficiently decayed. The combination of higher-order modes in the initial perturbation does however cause initial transient growth in the system, despite the corresponding eigenvalues being highly stable as can be seen by the inset in the main diagram where the linear response, predicted by (28), of a particular trajectory (initial condition is marked with a solid circle) is compared to the full system response.
The unstable A 2 state, as indicated by a large circular symbol, plays a crucial role in the partition of behaviours. We make the important observation that initial conditions with an interface length, L, sufficiently larger than the unstable state will eventually become unstable. In contrast, there are also initial conditions with ∆(0, A 2 ) < 0 (see initial conditions in the lower-half plane of panel(a)) which still become unstable due to the initial transient growth caused by the combination of higher-order modes; in direct contradiction to the simpler perturbation of the 'stretch' modes, where if ∆(0, A 2 ) < 0, then the system would relax to the stable state.
By examining the trajectories that explore the vicinity of the unstable branch in figure 16, it is not unreasonable to hypothesise that by continually refining the value of ρ in the initial perturbation the system would be able to stay in the vicinity of the unstable state for an arbitrary time-period; where the dynamics of the system are dominated by the stable eigenmodes of the unstable branch. This behaviour is typically indicative of interpreting the unstable branch as an edge state, commonly used to describe weakly unstable states in the transition to turbulence and other fluid dynamics problems [see 27, for a description of an edge state]. In these scenarios the weakly unstable edge state acts as an 'edge' between two dynamical outcomes; in our problem it separates the system returning to the stable state or the system diverging.
The role of the unstable state is further emphasised in figure 17 for perturbations from the unstable A 2 branch. In panel (a) different colour trajectories correspond to whether a positive or negative ρ is chosen in the IC given by (38), and panels (b) and (c) show the initial conditions compared to the A 2 steady state. Again, the system outcomes are partitioned by the presence of the unstable state which 'deflects' the system to either relax back to the stable state or the contact line velocity diverges. We note that in this case if ∆(0, A 2 ) < 0 (ICs in the lower half plane of panel (a)) then the system will become unstable which is again in direct contradiction to the 'stretch' mode perturbations. These calculations demonstrate that L is a more robust indicator of dynamical outcome; in this figure, all trajectories in the right side of the unstable state eventually experience instability contact line velocity divergence.

Physical Significance of the perturbations
We now relate these results to the physical system. Firstly, we make the key observation that the system is able to experience instability for Ca < Ca crit due to finite-amplitude perturbations of the stable state. Furthermore, these results suggest that the system is more susceptible to instability caused by perturbations that increase the total arclength rather than increase the meniscus rise; corrugations, representing 'noise', that increase the arclength, but not the meniscus rise, are more dangerous. Given that noise is a ubiquitous phenomena in physical systems we would expect to see this realised in a system with Ca < Ca crit . Secondly, the extent to which we are able to perturb the stable A 1 /R 1 state so that the system remains stable is dictated by the information encoded in the unstable A 2 /R 2 steady state, in particular whether the perturbation causes the interface length to increase beyond the length of the unstable steady state. A consequence of this is that the closer Ca is to Ca crit (from below), the smaller the finite amplitude is required, as the stable and unstable branches are increasingly approaching each other.
As a conclusion, finite disturbances from the steady stable state potentially lead to instability before that interface might become unstable otherwise (due to a lack of steady state). Minimising ambient disturbances and minimising fluctuations in Ca would help maintain a stable interface for higher Ca. The physical mechanisms underlying the observations about stability/instability are not straightforward to address due to the highly nonlinear nature of the problem. Moreover, three-dimensional perturbations will likely be important in practice and introduce additional physical mechanisms. Such 3D calculations are exceedingly challenging and are currently the focus of ongoing research, so we defer a detailed study of physical mechanisms for future work.

B. Dynamics 'beyond' the fold: Ca > Ca crit
We now turn our attention to starting the system from rest with a flat interface beyond the fold, i.e. Ca > Ca crit , so that we are able to answer the second question stated at the start of § IV. The IC is For values of Ca exceeding the critical value there are no (known) steady states which can influence the system, unlike in the previous section. For the advancing case, the system quickly diverges with the contact line speed diverging rapidly causing the calculations to fail to converge, as shown in figure 18 where the time-signal of Y is measured along with time snapshots of the interface at the times indicated for Ca = 1.05, 1.25. We emphasise that the system is not attracted to a different steady state of the system that may exist and we are unable to find any additional steady states beyond the fold bifurcation. Prior to the calculations ceasing to converge, the pressure gradients are particularly strong near the contact point and as explained earlier we would expect the fluid pressure gradient to exceed the capillary stress gradient so that the fluid phase is unable to pump fluid out near the contact line.
For the receding case the contact line velocity does not diverge and a thin-film develops.  The system, although not approaching a steady state, forms a coherent structure whose meniscus rise grows at a constant velocity, similar to the structure seen in figure 15(a). By examining the phase space, shown in figure 19(b), we can see that the system trajectory closely matches the steady solution curve and approaches a fixed value of Ca as t → ∞ which appears to be independent of Ca, indicating the contact-line region and the thin-film region are essentially de-coupled by this point. Unlike the advancing case, the time-dependent trajectories and the steady solution curve for the receding problem, when plotted in (Ca, Y ), are qualitatively similar. Figure 20 shows the comparison of the time-dependent trajectory with the steady solution curve for λ = 0.01, χ = 0.1, Ca = 0.15 > Ca crit as well as comparing the actual interface profiles for the time-dependent calculation and the corresponding steady solution for the same value of Y (inset figures). It is striking how close the two corresponding interface profiles are, although we note that in the time-dependent case the horizontal height of the inflection point is smaller than that of the steady solution in all of the cases. This again shows compelling evidence that in the receding contact line problem the unsteady solution branch represents the effective dynamics of the system.

V. DISCUSSION
We have developed a time-dependent hybrid model and utilised ideas from dynamical systems theory to investigate the advancing and receding contact-line problem. By solving a generalised eigenproblem numerically and performing IVP simulations we have demonstrated that far from being passive, the unstable branch of the bifurcation structure plays a subtle role in the underlying time-dependent behaviour, from a dynamical systems perspective.
By perturbing the stable branch using the eigenmodes we have demonstrated that the unstable branch represents the basin boundary of attraction of the stable solution and it has a profound effect on the eventual evolution of the system. We have demonstrated that perturbations that cause the interface to stretch are robust, in that provided the stretch does not exceed that of the unstable branch the system returns to the stable state. In contrast, perturbations that increase the overall length of the interface by adding 'corrugations' are more dangerous and the system instability is more susceptible to this class of perturbations. This information may be helpful in physical systems as a means of flow control as knowing the structure of the unstable eigenmodes may enable us to stabilise the system using suction/injection techniques. We note that the perturbations we have considered here are purely theoretical eigenmode perturbations and it would be of interest to examine the stability of the contact-line when physically perturbed by, for example, the effect of surface defects on the moving substrate. This is part of the authors' current research.
In addition, by performing time-dependent calculations we have shown that trajectories in phase space qualitatively match the steady bifurcation structure. In both the receding and advancing cases, the solution curve describes what [5] termed the 'effective dynamics' of the system. The bifurcation structure remains structurally stable (i.e. the stable-foldunstable branch structure does not change) as χ and λ is varied, (see figure 8) and hence we predict the overall qualitative behaviour to be similar, i.e. the unstable branch is generically the basin boundary of attraction for systems of this nature.
Gravity and inertial effects can easily be added to the model and it is interesting what effect these would have on the bifurcation structure and stability analysis. We have already performed some preliminary calculations incorporating gravity and have found that the bifurcation structure experiences the same multiple saddle-node bifurcations as predicted by the lubrication model [5]. Our preliminary analysis shows that the entire upper branch past the first saddle-node bifurcation is unstable, but we have not pursued this in detail and leave this for future research. For inertial effects, we could expect other types of bifurcations, including Hopf bifurcations, which would introduce complex eigenvalues to the spectrum of σ and lead to more complex transient behaviour. We also leave this avenue for future research.
In our model we choose the simplest approach and set the dynamic contact angle to be constant. It is hotly debated whether this is indeed physically realistic and whether molecular or hydrodynamic effects are stronger near the contact line. A natural extension to the model would be to investigate the effect of having a contact angle that is dependent on the velocity of the plate, such as that proposed by molecular kinetic theory (MKT), [see, for example 2,15] or the interface formation model [see 41]. Provided the bifurcation structure remains intact (i.e. stable-fold-unstable) then we predict the dynamics will be qualitatively the same. This is the subject of a submitted article where we apply the hybrid model and stability algorithm developed here to account for a Ca-dependent contact angle observed in molecular simulations [25].
For the receding contact line, the instability that occurs manifests itself as the formation of a thin film and a capillary ridge, but for the advancing case the calculations are significantly more difficult to get a converged solution and we are only able to advance in time until air entrainment first starts to occur. The instability in this study is taken in a two-dimensional context but it is well-known that the 'saw-tooth' patterns that emerge in an unstable advancing contact line are intrinsically three-dimensional. A natural extension to our model would be to consider transverse perturbations of the advancing contact line in the third dimension and perform a stability analysis, similar to our approach here. This has been done before in the receding case [43] and the advancing case [53] using a lubrication model, but the time-dependent hybrid model developed here is a perfect testing ground for a three-dimensional calculation as the computational cost is relatively small compared to the full model, and we are also able to fully account for the velocity in the liquid, rather than use a lubrication model. This direction of research is currently being developed.
To validate the hybrid model we performed a series of time-dependent IVP calculations starting the system from rest (u i = 0) with a flat interface. If we choose Ca < Ca Crit we expect the system to eventually relax to a stable state. As χ → 0 the full model and hybrid model should converge. Panel (a) of figure 21 shows the time-signal of the meniscus rise, Y , for different viscosity ratios (different colours) for the full model (dotted lines) and the hybrid model (solid lines). In all cases the interface eventually relaxes to a stable state, as shown by the time signals in the main figure. It is also clear that for the smallest value of χ in the figure the full and hybrid model are virtually indistinguishable. This is because for small χ, the fluid only has an influence on the liquid in thin films in front of the contact line, and this is where the lubrication model is valid. In contrast, at moderate χ the influence of the fluid is felt everywhere, the films are thicker/non-existent and this approximation loses accuracy. In other words, this approximation works best when the fluid is a gas.
Panel (b) of figure 21 shows the comparison of the solution curve in the receding problem for the full and hybrid model. As can be seen the value of Ca Crit changes only minimally for the hybrid model when compared to the full model.

Appendix B: Validation of eigenvalue calculation
We can validate the solutions of (27) against a situation where analytic eigenmodes are known. If we assume the wall is static (U = 0) and thatĤ/L 1, corresponding to a short fat pool of liquid at the bottom of the channel with no-slip beneath it, then we can also apply a lubrication approximation to the liquid domain. If the fluid is treated as a vacuum we have the following equation for the vertical height of the interface: where C is a known constant containing non-dimensional system parameters (in this case we choose C = 3). We choose the boundary conditions y(0) = y(L) = y 0 , y xxx (0) = y xxx (L) = 0 which describe a pinned contact line and no flux through the walls so that the resulting interface is flat in equilibrium. In this case the set of unknowns is 'one-dimensional' in that w = [y] and using the perturbation in (25) yields a set of equations where analytic progress can be made. It can be shown that the eigenvalues and eigenmodes of (B1) can be written as g n = a n sinh(σ These non-trivial analytic expressions can be used to validate our numerical stability calculations. We chose a computational domain with H/L = 40, to reflectĤ/L 1, and calculate the corresponding eigenmodes numerically using the hybrid model described in the previous section. Figure 22 shows the first three eigenmodes corresponding to the largest three eigenmodes, n = 1, 2, 3 as calculated by the hybrid model and the analytic solution (different colour curves). These solutions are stable as the eigenspectra lies exclusively in the left-hand complex plane and the agreement between the numerical calculations and the analytic solution obtained from the lubrication model is excellent, giving us confidence in our computational framework.