Analytical and Numerical Results on the Positivity of Steady State Solutions of a Thin Film Equation

We consider an equation for a thin-film of fluid on a rotating cylinder and present several new analytical and numerical results on steady state solutions. First, we provide an elementary proof that both weak and classical steady states must be strictly positive so long as the speed of rotation is nonzero. Next, we formulate an iterative spectral algorithm for computing these steady states. Finally, we explore a non-existence inequality for steady state solutions from the recent work of Chugunova, Pugh,&Taranets.


Introduction
Thin liquid films appear in a wide range of natural and industrial applications, such as coating processes, and are generically characterized by a small ratio between the thickness of the fluid and the characteristic length scale in the transverse direction. The physics of such processes are reviewed in Oron, Davis & Bankoff, [16] and more recently in Craster & Matar, [13]. If one includes surface tension in the model, the fluid is governed by a degenerate fourth order parabolic equation. Such equations, studied by Bernis & Friedman, [8], Bertozzi & Pugh, [9], and Beretta, Bertsch & dal Passo [7] still hold many challenges. See, for example, [12,23,24] for progress on models which include convection, and Becker & Grün, [5], for a thorough review of analytical progress.
There has been recent interest in equations governing the evolution of a thin, viscous film on the outer (or inner) surface of a cylinder of radius R rotating with constant angular velocity ω, as in Figure 1. Subject to non-dimensionalization, this can be modeled by (1) u t + (u n (u θθθ + u θ − sin θ) + ωu) θ = 0, Of particular interest is the case n = 3, in which one assumes there is no slip at the fluid-cylinder interface. The solution, u, measures the (dimensionless) thickness of the fluid, parameterized by the angular variable, θ ∈ Ω ≡ [−π, π). To derive (1), one begins with the Navier-Stokes equations, together with a kinematic boundary condition for the free surface surface tension. Subject to appropriate physical scalings and geometric symmetries, we are left with a dimensionless equation with just the ω parameter. We refer the reader to, for instance, [20] and [19] for a full derivation. Though Buchard, Chugunova, & Stephens, [10], gave a detailed examination of (1) when ω = 0, our understanding of the case ω = 0 remains incomplete. Some g ω R R + u(θ, t) θ O Figure 1. The grey region represents the thin film of liquid on the surface of the cylinder, which is rotating with constant angular velocity ω. We assume symmetry in the z coordinate, so that u is a function of t and θ only.
partial results on this case have made assumptions that solutions of (1) are strictly positive, such as in [22], where Pukhnachov proved, under a smallness assumption, positive solutions of (1) exist and are unique.
In this work, we further study steady state solutionsnon-negative, time independent solutions u(θ) of (1), satisfying (2) (u n (u θθθ + u θ − sin θ) + ωu) θ = 0 This can be integrated once, to where q, the constant of integration, is the dimensionless flux of the fluid. In particular, we give an elementary proof that when ω = 0, both weak and classical solutions must be strictly positive; see Section 2. For a more extended study of the steady states, it is helpful to explore the problem numerically. This can provide information on how, and if, u and its derivatives develop singularities as the rotation parameter, ω, goes to zero. For such a task, one needs a robust algorithm for computing solutions of (3). In Section 3, we formulate and benchmark an iterative spectral algorithm that is easy to implement and rapidly converges.
Our last result, appearing in Section 4, is motivated by Chugunova, Pugh, & Taranets, [12], who proved that for n = 3 there can be no positive steady state solutions when the rotation and flux parameters satisfy the inequality This was refinement of an earlier bound due to Pukhnachov, [21]. We explore the (ω, q) phase space and find evidence that this bound may be sharp. Furthermore, it appears that this curve is related to a transition between two different regimes of steady state solutions.

Notation.
In what follows, we shall use the following notational conventions. C k per (Ω) is the set of continuous 2π-periodic functions on Ω = [−π, π) with k continuous derivatives. H k per (Ω) is the set of square integrable 2π-periodic functions on Ω = [−π, π) with k square integrable weak derivatives. For a function φ in one of these spaces, we shall denote its mean byφ. For brevity, we will also sometimes write a b to indicate that there exists a constant C > 0 such that a ≤ Cb.

Positivity of Steady States
Here, we prove that both classical and weak steady state solutions of (1) are strictly positive provided ω = 0. Buchard, Chugunova, & Stephens, [11], demonstrated that for ω = 0, such solutions need not be positive.
2.1. Classical Steady State Solutions. As noted, the case n = 3 corresponds to a no-slip condition at the surface of the cylinder. With this in mind, it is physically obvious that, as long as ω = 0, any steady state solution of (1) with n = 3 will coat the entire cylinder -that is, there will be no "dry patches" on the surface of the cylinder. The following elementary argument shows that this physical intuition is correct, and is in fact valid for a broad range of exponents n.
On the support of u, Since u θθθ is bounded as we send θ towards the root, ω = 0.
Though a classical solution of (1) needs to be C 4 , the above proof succeeds without change if u is merely C 2 with bounded third derivative.
We first establish some pointwise properties of a solution u, in the neighborhood of a hypothetical touchdown point, θ 0 . Lemma 1. Let u ∈ H 2 per (Ω) be non-negative. If there exists θ 0 such that u(θ 0 ) = 0, then Proof. By virtue of a Sobolev embedding theorem in dimension one, H 2 → C 1,1/2 , [2,15]. Therefore, both u and u are Hölder continuous with exponent 1/2; Because u is C 1 and non-negative, we must have that u (θ 0 ) = 0 too, otherwise there would be a zero crossing.
Refining (9), we can apply the mean value theorem and the Hölder continuity of u to get where θ lies between θ and θ.
To prove the weak form of Proposition 1, we begin by deriving a weak analog of (3): Lemma 2. If u is a weak solution in the sense of (7), then there is a constant q, depending only on u, such that Proof. First, observe that for any φ ∈ C ∞ per (Ω), φ −φ has mean zero and therefore has an antiderivative in C ∞ per (Ω). Given a test function φ, let ϕ solve ϕ θ = φ −φ. Plugging in ϕ in (7) in place of φ, we have We may now take We next derive a weak analog of equation (5) in the event that u has a root.
Proof. Given ε > 0, let φ ∈ C ∞ per (Ω) be a compactly supported test function satisfying: The constant C depends on the L 2 norms of u and its second derivative, which are both finite by assumption.
Using estimates (8), we have that within I ε , Substituting this pointwise analysis into (14), Since this holds for all ε > 0, we conclude that q = 0 for n > 1.
Using Lemmas 2 and 3, we now know that a weak solution with a touchdown satisfies Contrast this with its classical counterpart, equation (5).
We now proceed to our main result in the weak case, under the additional restriction that n ≥ 2: Proposition 2 (Positivity of Weak Solutions). Let u be a weak solution satisfying (7) with n ≥ 2 for all ϕ ∈ C ∞ per (Ω). If S = {θ : u(θ) > 0} = Ω, then ω = 0. Proof. The proof is by contradiction. Assuming that S = Ω, let θ 0 be a point on the boundary of S. Given ε > 0 sufficiently small, take y ∈ Ω such that the interval is contained in S and |θ 0 − y| = ε. Now, let ψ be a compactly supported test function as in (12), and let This function then has support in I ε/2 as in Figure 2. Finally, define φ ε = ψ ε /u n , which is as smooth as u, by the conditions on ε and y.
Plugging φ ε into (17), we see: Since u ε 3 2 on the interval and n > 1, we can bound the integral on the left from below: By the construction of ψ ε , ψ ε θ L 2 (I ε/2 ) ε −3/2 . Therefore, (20) |ω| ε For n > 2, we immediately observe that since ε > 0 is arbitrary, ω = 0. For n = 2, we know that since (u θθ ) 2 and u 2 are measurable, In the proof of proposition 2, we choose a smooth, compactly supported function ψ with support contained in the support of u, but close to a zero of u.

A Computational Method for Steady State Solutions
While we have given a simple proof that when ω = 0,steady state solutions must be positive, we have not established their existence. For that, we rely on [22], where it was shown that in the case n = 3 "small solutions" existed and were unique. Let us motivate this with some heuristic asymptotics. First, assume that 0 < q ω, i.e.
Subject to this assumption, we may then perform a series expansion in , we then match orders of . At order 0 , As we have no way to solve this, we shall assume the order O( 0 ) term is trivial, . Proceeding with this assumption, the next few orders are: The leading order approximation of the steady state is thus This corresponds to the "strongly supercritical regime" discussed by Benilov et al., [3,6]. A similar analysis applies to the case n = 2, for which an approximate steady state solution is This regime was studied in [22], where the expansion was justified by a contraction mapping argument. This motivates the numerical implementation of this approach as a tool for computing such steady states.

Description of the Algorithm.
In what follows, we focus on the case n = 3, although it can be easily adapted to other cases. We seek to rewrite (3) in the form where L is a linear operator, f contains the nonlinear terms and g is a a driving term. This will motivate the iteration scheme To begin, we divide (3) by u 3 , to obtain (∂ 3 θ + ∂ θ )u = sin θ + qu −3 − ωu −2 . As ∂ 3 θ + ∂ θ has a non trivial kernel, the lefthand side cannot be inverted without careful projection. Pukhnachov resolved this problem by adding and subtracting appropriate terms to the equation. First, we introduce the variable v, v ≡ u − q ω .

Then
Lv The right hand side of (28) contains a driving term, sin θ, and an essentially nonlinear term of order O(v 2 ). Rearranging (28), In Fourier space, this is Plots of solutions computed using this algorithm appear in Figures 3 and 4.

3.2.
Performance of the Spectral Iterative Method. When our algorithm successfully converges, it does so quite rapidly, as shown in Figure 5, where we plot the error between iterates v (j) and the final iterate. We see that the error is proportional to e −αN , although it is clear that the rate, α, varies with q and ω. Indeed, as shown in Figure 6, there are choices of these parameters for which our algoirthm diverges; this is consistent with the threshold (4). While this will  be further discussed in Section 4, we belive this is closely related to the transition between the "small amplitude regime", and an essentially nonlinear regime.
However, when we do observe convergence, there is a very high degree of precision, even using a modest number of grid points. Figure 7, shows how quickly the maximum and minimum of the numerical solution converge to one with 10 5 grid points as the number of grid points increases. Again, solutions with values of q and ω closer to (4) converge less rapidly; the solution in black is much closer to this threshold than the other three.   Benilov, Benilov, & Kopteva computed the solutions directly by a finite difference discretization with a Newton algorithm, [6].
An advantage of our approach is that it converges quickly and is extremely easy to program. An implementation of the algorithm in Matlab is available at http://hdl.handle.net/1807/30015.
Our approach bears some resemblance to Petviashvilli's algorithm and the related spectral renormalization procedure, which have become quite popular within the nonlinear dispersive wave community for computing solitons, [1,17,18]. For such equations, including Korteweg -de Vries and nonlinear Schrödinger/Gross-Pitaevskii, solitary wave solutions solve semilinear elliptic equations of the form For power nonlinearities, f (Q) = Q p , Petviashvilli's algorithm finds a solution by seeking a fixed point of an iteratively renormalized nonlinear equation in the Fourier domain, is the j-th approximation, and Λ j is a renormalization constant at this iterate. Our scheme, (31) lacked any such renormalization constant, Λ j . This raises the question of whether or not such a constant could improve our approach. Our results show this is not the case. Introducing the constant we can alter our iteration scheme to be where γ is a constant which we may alter to stabilize or accelerate the scheme. An outcome of our experiments is that not only does our method converge with γ = 0, there does not appear to be appreciable benefit with any other value; see Table 1. For this reason, we do not employ Λ j , and only deem our algorithm an iterative spectral method. Benilov, Benilov & Kopteva, [6] tried fixing ω, and numerically explored the relationship between the flux and the mass as they varied. They found non-uniqueness of the solutions, in the sense that for a certain range of M , there were multiple solutions with different values of q. Alternatively, for a certain range of q, there were multiple solutions with different values of M . See Figure 14 of their work for a phase diagram, and see our Figure 10 for a similar diagram. This non-uniqueness affects our algorithm. For a given (ω, q), the solution that the method converges to may not be the solution we desire, as the following example shows. Using, AUTO, [14], if we seek a solution with M = 1 and ω = .09, we will obtain the solution pictured in Figure 8 (a). This solution has flux parameter q = 0.0114039.
Suppose we now fix ω = .09, and continue the solution by lowering the mass. As seen Figure 8 (b), there will be a second solution, with the same q, but slightly smaller mass. For q = 0.01140392, the second solution is at M = 0.912782 Thus, for a given (ω, q), there may be multiple solutions, each with a different mass-which solution does the iterative spectral algorithm converge to? As Figure 8 (a) shows, it converges to the solution with the smaller mass. The agreement between the AUTO solution and our method's solution is quite good; after the AUTO solution is splined onto the regular grid of the spectral solution, the pointwise error is O (10 −7 ).
The inability of our algorithm to recover the M = 1 solution can be traced to the development of a linear instability in the iteration scheme. Let v (j) = v + p (j) , where v is a solution, and p (j) is a perturbation at iterate j. If we linearize (30)  about v, we have If the spectrum ofL v lies strictly inside the unit circle, the perturbation will vanish, and the solution v will be stable with respect to the iteration. However, if there are points in the spectrum of this operator which lie outside the unit circle, the solution will be unstable to the iteration, and we would expect our algorithm to diverge. Though we shall not attempt to prove properties of the spectrum ofL, we can discretize our problem and compute eigenvalues of the induced matrix. Computing the eigenvalues of a discrete approximation ofL v at the two AUTO solutions, we find that the matrix corresponding to the M = 1 solution does, in fact, have an eigenvalue outside the unit circle, λ unstable = 1.017059; see Figure 9. The discretized operator is defined as with the differentiation matrix D, a dense Toeplitz matrix for band limited interpolants, as in [25]. The value of the unstable eigenvalue converges quite quickly as we increase the number of grid points, as shown in Table 2.
The (M, q) phase space of solutions with fixed ω can be more tortuous than that shown in Figure 8 (b). Indeed, as we send ω → 0, a singular limit, there may not just be multiple solutions of different mass and the same flux, but also multiple solutions of different flux and the same mass. Such an example appears in Figure  10. These cases are more extensively explored in [6] and in Badali et al, [4].    where the curves of constant mass cross, there will be two or more solutions with the same ω and q values. All curves of solutions are constrained by (4), as expected. Also note that at fixed q, for sufficiently large ω, the relationship is nearly linear. Computed using AUTO.

Non-Existence Results
As mentioned in the introduction, in [12], the authors prove that solutions to (3) do not exist for all (ω, q) pairs. Indeed, they show there are no solutions when (4) holds. Interestingly, this relation is stated entirely in terms of q and ω, and makes no consideration of the mass, M .
We now explore, numerically, (4). For fixed mass, we compute the associated solution at a variety of ω values, determining q. The results are appear in Figure  11. As (ω, q) cannot uniquely determine u, we rely on AUTO and its continuation algorithms, rather than our own iterative scheme.
Examining the figure, we remark on the two features found in each fixed mass curve. Far to the right of the non-existence line, the ω-q relationship appears to be linear. Indeed, for moderate mass values, the relationship appears to be of the form (21). It is in this "small amplitude" regime that our iterative spectral algorithm succeeds.
As ω continues to decrease, the scaling changes to accomodate (4), and appears to follow a curve like q ∝ ω β , for a value of β > 3/2. For all masses, the deviation from (4) is less than an order of magnitude. Thus, there is not an obvious small parameter from which to formulate a series expansion. We also remark that there are many intersections between the curves of constant mass; the lack of uniqueness of a solution for a given (ω, q) is quite common. As ω and q approach this essentially nonlinear regime, our algorithm becomes limited. As discussed in the preceding, section unstable eigenvalues can appear, precluding convergence to solutions with a desired mass, and the rate of convergence slows.
Turning to Figure 12, we see the tendency of the solutions to form singularities, vertical asymptotes, as ω tends towards zero for fixed q. These were computed using our spectral iterative algorithm; the asymptotes occurred where the algorithm failed to converge in a reasonable number of iterations. As these figures show, it is possible to get quite close to the theoretical threshold of (4). This suggests that the constraint may be sharp, and this warrants further investigation. Even if the constant is not correct, the scaling appears to be sharp, and the inequality plays an important role in separating two physical regimes.  10 as ω is varied. The dashed lines denote the observed threshold at which our spectral iterative algorithm fails and the theoretical threshold, (4) beyond which no solutions should exist. Computed using the spectral iterative method.

Discussion
In this work, we have given a simple proof of the positivity of steady state solutions of a thin film equation with rotation, subject to modest assumptions, provided n ≥ 2. It may be possible to employ additional regularity properties of the steady state solution to lower this threshold; we only relied on Sobolev embeddings in our result.
We have also formulated a very easy to implement algorithm for computing the steady state solutions for a given (ω, q) pair. However, in cases where multiple solutions are known to exist, an instability can develop about one of the solutions, precluding convergence to this solution. We conjecture that the preference is related to the proximity to the non-existence curve, where the parameters have left the small amplitude regime. Finding a way to stabilize the spectral algorithm in this regime remains an open problem.
Finally, our computations suggest further exploration of the non-existence curve is warranted. Figures 11 and 12 suggest that the nonexistence relationship may be sharp. Furthermore, the phase diagram indicates that this inequality separates two physical regimes. In the first, where ω and q are nearly linearly related, the solution is nearly constant, and can be approximated by (26). In order to accomodate the nonexistence constraint, the phase diagram bends, and we enter a regime where the parameters are nonlinearly related, and the steady state solutions cease to be well described by unimodal bumps.