Steady states and dynamics of a thin film-type equation with non-conserved mass

TL;DR

This paper investigates the steady states and dynamics of a one-dimensional thin film equation with non-conserved mass, revealing novel behaviors such as multiple steady states, bifurcations, rupture singularities, and limit cycles.

Contribution

It introduces a modified thin film model breaking gradient flow structure, uncovering new solution classes and complex dynamics not seen in conserved cases.

Findings

Existence of two classes of coexisting steady states.

Identification of bifurcation structures and stability regimes.

Observation of finite-time rupture and limit cycle phenomena.

Abstract

We study the steady states and dynamics of a thin film-type equation with non-conserved mass in one dimension. The evolution equation is a nonlinear fourth-order degenerate parabolic PDE motivated by a model of volatile viscous fluid films allowing for condensation or evaporation. We show that by changing the sign of the non-conserved flux and breaking from a gradient flow structure, the problem can exhibit novel behaviors including having two distinct classes of coexisting steady state solutions. Detailed analysis of the bifurcation structure for these steady states and their stability reveal several possibilities for the dynamics. For some parameter regimes, solutions can lead to finite-time rupture singularities. Interestingly, we also show that a finite amplitude limit cycle can occur as a singular perturbation in the nearly-conserved limit.

Figures

Equations

$$\mathcal{E}[\varphi ]={\iint }_{D}U(\varphi )+\frac{1}{2}|\nabla \varphi {|}^{2}dA.$$

$$\mu [\varphi ]\equiv \frac{\delta \mathcal{E}}{\delta \varphi }=U^{\prime }(\varphi )-{\nabla }^2\varphi .$$

$$\frac{\partial \varphi }{\partial t}=\nabla \cdot \left(M_{\mathrm{c}}(\varphi )\nabla \mu [\varphi ]\right)-M_{\mathrm{nc}}(\varphi )\mu [\varphi ],$$

$$\mathcal{M}={\iint }_D\varphi dA\frac{d\mathcal{M}}{dt}=-{\iint }_D{M}_{\mathrm{nc}}\mu dA.$$

$$\frac{d\mathcal{E}}{dt}=-{\iint }_D{M}_{\mathrm{c}}|\nabla \mu {|}^2+M_{\mathrm{nc}}{\mu }^{2}dA\le 0,$$

$$\frac{\partial \varphi }{\partial t}=D\frac{{\partial }^2\mu }{\partial x^2}-\mu ,\mu =U^{\prime }(\varphi )-\frac{{\partial }^2\varphi }{\partial x^2}.$$

$$\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(h^3\frac{\partial p}{\partial x}\right)-\frac{\beta }{h+K}p,p=\Pi (h)-\frac{{\partial }^{2}h}{\partial x^2},$$

$$\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(h^3\frac{\partial }{\partial x}\left[\Pi (h)-\frac{{\partial }^{2}h}{\partial x^2}\right]\right),$$

$$\frac{\partial h}{\partial t}$$

$$p$$

$$\Pi (h)$$

$$\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(h^3\frac{\partial }{\partial x}\left[\Pi (h)+\frac{\gamma }{2h^2}-\frac{{\partial }^{2}h}{\partial x^2}\right]\right)+\gamma \Pi (h),$$

$${\partial\phi\over\partial t}=\nabla\cdot(M_{\mathrm{c}}{}(\phi)\nabla\hat{\mu})+\hat{S}(\phi)\quad\mbox{with $\hat{\mu}=\hat{U}^{\prime}-\nabla^{2}\phi$},$$

$$\mathcal{E}={\int }_0^L\frac{1}{2}{\left(\frac{\partial h}{\partial x}\right)}^2+U(h)dx,\frac{d\mathcal{E}}{dt}=-{\int }_0^L{h}^3{\left(\frac{\partial p}{\partial x}\right)}^{2}dx+\gamma {\int }_0^L{p}^{2}dx,$$

$$\mathcal{M}={\int }_0^{L}hdx,\frac{d\mathcal{M}}{dt}=\gamma {\int }_0^L\Pi (h)dx.$$

$$0=\frac{d}{dx}\left(h^3\frac{dp}{dx}\right)+\gamma p,p=\Pi (h)-\frac{d^{2}h}{d{x}^2},$$

$$\frac{d^{2}h}{d{x}^2}-\Pi (h)=0.$$

$$\Pi (h)=0$$

$$U(h)=\int \Pi (h)dh=-\frac{1}{2h^2}+\frac{1}{3h^3}-\bar{\mathcal{P}}h.$$

$$\Pi (h_{\text{peak}})={\mathcal{P}}_{\max }-\bar{\mathcal{P}},{\mathcal{P}}_{\max }=\frac{27}{256}>0.$$

$${\bar{H}}_{-}=1+\bar{\mathcal{P}}+4{\bar{\mathcal{P}}}^2+O({\bar{\mathcal{P}}}^3).$$

$${\bar{H}}_{+}=$$

$$\lambda =\left[-M_c(\bar{H}){\left(\frac{2k\pi }{L}\right)}^2+\gamma \right]\left[{\Pi }^{\prime }(\bar{H})+{\left(\frac{2k\pi }{L}\right)}^2\right],$$

$${\bar{\ell}_{+}}^{\,\gamma}=\frac{2\pi\bar{H}_{+}^{3/2}}{\gamma^{1/2}}\qquad\mbox{for $\bar{H}=\bar{H}_{+}$,}$$

$${\bar{\ell}_{-}}^{\,\gamma}=\frac{2\pi\bar{H}_{-}^{3/2}}{\gamma^{1/2}}\qquad\mbox{for $\bar{H}=\bar{H}_{-}$.}$$

$${{\bar{\ell }}_{+}}^p=\frac{2\pi }{\sqrt{-{\Pi }^{\prime }({\bar{H}}_{+})}}.$$

$$\frac{dh}{dx}=s\frac{ds}{dx}=\Pi (h).$$

$$H_{\max}=(6\bar{\mathcal{P}})^{-1}+1+O(\bar{\mathcal{P}})\qquad\mbox{for $\bar{\mathcal{P}}\to 0$}$$

$$\ell (h_{\min })={\int }_0^{\ell (h_{\min })}dx=2{\int }_{h_{\min }}^{h_{\max }}\frac{1}{\sqrt{2U(h)-2U(h_{\min })}}dh.$$

$$⟨h⟩=\frac{1}{L}{\int }_0^{L}h(x)dx,$$

$$H_A(x)\ge {\tilde{H}}_A(x)=\max \left(h_{\max }+\frac{\Pi (h_{\max })}{2}{\left(x-\frac{1}{2}L\right)}^2,{\bar{H}}_{-}\right).$$

Full text

1 Introduction

Phase field models are widely used in many branches of continuum mechanics to describe pattern formation and evolving interfaces with respect to an order parameter $\phi({\bf x},t)$ that could represent a volume fraction of one phase in a binary mixture [24, 48, 18], or the density of healthy tissue in a biological system [21], for example.

Such models can be formulated in terms of a free energy, which in simple cases will have contributions from a potential energy of homogeneous phases, $U(\phi)$, and an interfacial energy from the formation of gradients,

[TABLE]

A chemical potential can then be defined from the variational derivative of the energy,

[TABLE]

As described by Thiele in [60, 63], a general dissipative evolution equation for $\phi$ can then be written as

[TABLE]

where $M_{\mathrm{c}}(\phi)$ and $M_{\mathrm{nc}}(\phi)$ are non-negative mobility functions for the conserved and non-conserved parts of the dynamics, respectively. These coefficients define a spatial flux following Fick’s law, ${\bf J}=-M_{\mathrm{c}}\nabla\mu$, and a mass loss rate,

[TABLE]

This model monotonically dissipates the energy as a Lyapnunov functional,

[TABLE]

which shows that all equilibrium states must correspond to constant chemical potentials, and further if mass is not conserved ($M_{\mathrm{nc}}\neq 0$), they must have $\mu\equiv 0$. We note that the $U(\phi)$ term in the energy may include both an intrinsic potential energy $\tilde{U}(\phi)$ and a linear term, $U(\phi)=\tilde{U}(\phi)-\mu_{0}\phi$. This leads to $\mu[\phi]=\tilde{\mu}[\phi]-\mu_{0}$, where $\tilde{\mu}$ represents an intrinsic chemical potential, so that for equilibrium states in the non-conserved system the condition $\mu=0$ corresponds to a constant intrinsic potential with $\tilde{\mu}=\mu_{0}$.

One model fitting into this framework is the mixed Allen-Cahn/Cahn-Hilliard equation (AC/CH) [40, 27, 26, 33] that describes microscopic pattern formation mechanisms like surface diffusion, adsorption and desorption [31]. The basic form of this model has constants for both mobilities $M_{c}$ and $M_{nc}$ ($M_{\mathrm{c}}\equiv D>0$ and $M_{\mathrm{nc}}\equiv 1$) and the potential energy functional is a symmetric double well, $U(\phi)=(1-\phi^{2})^{2}$, written in one spatial dimension as

[TABLE]

The classic Cahn-Hilliard equation [48, 47] for phase separation of binary mixtures is recovered by eliminating the non-conserved flux, $M_{\mathrm{nc}}\equiv 0$. Many variants of the mixed AC/CH model have been used to study more complex interfacial dynamics, one such example is the coupled system describing a model of tumor growth investigated in [21]. Studies have also addressed mathematical properties of solutions when the mobility coefficient is a degenerate function of $\phi$ in the Cahn-Hilliard equation [13] and for the AC/CH model [66]. There is a vast body of literature on existence, regularity and attracting states for the Cahn-Hilliard equation [40]. Some results on steady states in convective Cahn-Hilliard models subject to external fields have also been obtained [34]. Fewer results have been obtained specifically for the AC/CH equation [27, 32, 33, 66], but many extensions with systems of equations having phase-field models coupled to other physical effects have been used in applications [9]. The Cahn-Hilliard equation with other forms of non-conserved terms (not proportional to the chemical potential) have been used to describe linear adsorption/desorption (the Cahn-Hilliard-Oono equation) and other applications [40], but will not have an energy dissipation result like (2).

While the derivation of lubrication models for coating flows of thin viscous films on solid substrates follows from an asymptotic reduction of the Navier-Stokes equations for low Reynolds numbers [53, 11, 46, 50], the resulting thin film equations also fit the form (1) [41, 60, 61, 63, 62], with the order parameter representing the film height, $\phi\to h(x,t)\geq 0$. Likewise, the role of the chemical potential is taken on by the hydrodynamic pressure $\mu\to p$, which will have contributions from surface tension and a hydrostatic potential function.

One class of models derived for evaporating thin films of volatile fluids on heated substrates [39, 2, 1, 3, 42] is consistent with the form of equation (1), also see references listed in [61]. For Newtonian fluids on solid substrates with no-slip boundary conditions, the mobility functions for the conserved and non-conserved terms are $M_{\mathrm{c}}(h)=h^{3}$ and $M_{\mathrm{nc}}(h)=\beta/(h+K)$, where the constants $\beta>0$ and $K>0$ are set by material properties of the fluid and a thermodynamic-kinetic condition at the fluid-vapor interface [1]. Alternative forms for the mobility $M_{c}(h)$ have been discussed for slip models in [44]. The $M_{c}(h)$ function makes thin film models comparable to degenerate-mobility Cahn-Hilliard equations [16, 13, 66]. The resulting model is

[TABLE]

where $\Pi(h)=U^{\prime}(h)$ is a generalized disjoining pressure that incorporates both spatial wetting properties of the substrate and the thermodynamics driving evaporation, $\Pi(h)=\tilde{\Pi}(h)-\bar{\mathcal{P}}$ [30]. The contribution from the standard disjoining pressure function $\tilde{\Pi}(h)$ gives the pressure due to the hydrophobic or hydrophilic properties of the substrate, as typically described in relation to contact angles of droplets [22, 14, 58]. For partially-wetting liquids, the form $\tilde{\Pi}(h)=h^{-3}-h^{-4}$ has been frequently used and provides a positive lower bound for the film thickness [5, 23, 30]. The pressure offset $\bar{\mathcal{P}}$ then encodes the difference between the fluid temperature and the surrounding vapor, determining whether evaporation or condensation is favored [30]. For $\beta\geq 0$ since this model matches (1) it will have an energy dissipation equation of the form (2).

In the special case $\beta=0$, the model (3) reduces to the mass-conserving thin film equation where no evaporation or condensation occurs,

[TABLE]

which has been studied extensively in mathematical fluid dynamics [43, 11, 53, 46] and PDE analysis [37, 38, 49, 56, 6]. For $\beta=0$, [5] gives a proof showing that given appropriate initial conditions, there is global existence of smooth positive solutions for all times. This proof can be extended to (3) with $\beta>0$ [30, 28]. However these proofs do not hold for $\beta<0$ and hence global existence is not guaranteed for that case.

There have also been studies using other physically-based models of evaporating fluid films [8, 52, 51, 7, 64, 61] that have used forms of the evaporative flux that are not consistent with (1). For instance, the linear evaporative flux used in [7] can yield a pattern-forming instability. In [64] a base model for an evaporating film of the form (1) was used, but the addition of an inhomogeneous term representing an influx of fluid from a porous substrate to maintain steady droplet solutions, breaking the variational form. Some PDE analysis has also addressed unstable thin film equations with linear or nonlinear absorption (or “proliferation”) terms [19, 36, 57, 20]. For these equations the evolution of the energy does not follow (2) and different approaches to the analysis of solutions must be used for each model [17].

In [29] it was shown that if the non-conserved term in (3) is modified to have the opposite sign, then the PDE yields finite-time rupture singularities driven by the non-conserved loss term. Here we will more fully explore the other consequences that occur in a model that breaks from the gradient flow framework of (1).

In particular, we study the nonlinear partial differential equation in one dimension on a periodic domain $0\leq x\leq L$,

[TABLE]

with $\bar{\mathcal{P}}$ being a pressure-offset constant and $\gamma>0$ scaling the strength of the non-conserved effects.

This system retains the mobility $M_{\mathrm{c}}=h^{3}$ and the generalized disjoining pressure $\Pi(h)$ from the volatile thin film equation (3), but uses a negative constant for the mobility $M_{\mathrm{nc}}=-\gamma$; apart from the sign-change, this is analogous to (3) with $K\gg h$. For a thorough discussion of mass-conserving and mass non-conserving models in different contexts with variational or non-variational contributions, we refer readers to [17]; our model with $\gamma>0$ belongs to their mass non-conserving non-variational case.

Note that (5) can also be written as

[TABLE]

namely as a Cahn-Hilliard equation with a modified potential and an added source term. It can be shown that to recast (1) in the form

[TABLE]

the mobility $M_{\mathrm{nc}}{}$ must be a constant. To define a gradient flow with energy dissipation (2), the source and modified potential must be related by $\hat{U}^{\prime\prime}=-\hat{S}^{\prime}/M_{\mathrm{nc}}{}+M_{\mathrm{nc}}{}/M_{\mathrm{c}}{}$ for some $M_{\mathrm{nc}}{}>0$. This is a less convenient form than (1), but again serves to show that added lower order terms can change significantly the properties of the mass-conserving model. Rather than attempting to describe the influence of broad classes of source terms, we formulated (5) to take advantage of the conserved and non-conserved fluxes both being expressed simply in terms of a single pressure functional, but showing the consequences of not being a gradient flow.

Starting from positive, finite-mass initial data $h_{0}(x)>0$, we will show that the interactions between the conserved and non-conserved terms in (5) yield novel behaviors that are not possible in (3). In particular we will see that there exist two distinct classes of nontrivial steady states, and that finite-time rupture and limit cycles can occur.

The structure of this paper is as follows. In Sec. 2 general properties and the steady state structure of the PDE (5) are outlined. In Sec. 3 spatially uniform states and their instabilities with respect to spatial perturbations are investigated. Properties of nonuniform second-order and fourth-order states and the corresponding bifurcations will be discussed in Sec. 4 and Sec. 5. Stability of these states will be analyzed in Sec. 6, followed by discussions of limit cycle dynamics in Sec. 7 and the formation of finite time singularities in Sec. 8. Numerical simulations showing dynamic transitions among different steady states are presented in Sec. 9. Concluding notes and a discussion of the remaining open questions are presented in Sec. 10.

2 Properties of PDE (5)

The energy functional and its rate of change for equation (5) are given by

[TABLE]

where $p=\delta\mathcal{E}/\delta h$ is given by (5b) and the potential $U=\int\Pi(h)dh$ with $\Pi(h)$ defined in (5c). For $\gamma>0$, (7)2 shows that the contribution from the second integral can overcome the dissipation from the first integral, and hence the energy can evolve non-monotonically.

The total mass of the solution and its rate of change are

[TABLE]

This indicates that if $\Pi(h)\leq 0$ (or $\Pi(h)\geq 0$) everywhere in the domain, then $\gamma>0$ yields a monotone decreasing (or increasing) mass over time. When $\Pi(h)$ changes sign, the evolution of the mass depends on the form of the solution $h$ and the system parameter $\bar{\mathcal{P}}$, similar to the volatile thin film model (3) studied in [30].

Now we focus on the positive steady states $H(x)>0$ of the model (5). Unlike most mass-conserving thin film equations, the interplay between the conserved and non-conserved terms in the equation (5) leads to more interesting sets of equilibria. By setting the time derivative in (5) to be zero, we get the ODE system for all steady state solutions of the model as

[TABLE]

which we study subject to periodic boundary conditions. This system is equivalent to a fourth-order nonlinear ODE for $h$. If we set the pressure $p\equiv 0$ then (9) reduces to the second-order differential equation

[TABLE]

By setting the derivative term in (10) to be zero, we can further reduce the second-order equation to an algebraic equation

[TABLE]

for spatially-uniform steady states. In particular, we will refer to these steady states as

• •

Uniform steady states, $\bar{H}$, satisfying (11).

• •

Second-order nonuniform steady states, $H_{A}(x)$, satisfying (10) but not (11).

• •

Fourth-order nonuniform steady states, $H_{B}(x)$, satisfying (9) but not (10).

There have been systematic studies of second-order steady states of both volatile and nonvolatile thin film equations [55, 5, 60]. Investigations on fourth-order steady states are given in [64, 17]. However, the present study is the first, to our knowledge, to classify equilibrium solutions of a thin film-type equation into both second- and fourth-order states. In our framework, the second-order states $H_{A}$ and the uniform states $\bar{H}$ are associated with a constant zero pressure $p$, and the fourth-order states $H_{B}$ have nontrivial pressures.

The system parameter $\gamma$ is important for the existence and multiplicity of steady states of the model. For $\gamma\leq 0$, only the uniform and second-order nonuniform steady states can exist [30]. That is, since the energy (7) is a Lyapunov functional in this case, the steady states are given by extrema of the energy. Since the contributions from each integral in (7)2 are non-negative, they must equal to zero independently when the extrema of the energy is attained. Namely by setting ${d\mathcal{E}}/{dt}=0$ in equation (7)2, from the first integral one obtains ${\partial p}/{\partial x}=0,$ which indicates that $p$ is a constant over the domain. Equivalently, for $\gamma=0$ the governing equation (5) reduces to a conserved equation, and the equilibria are given by solutions with a constant pressure, $p(x,t)\equiv P$. For $\gamma<0$ the second integral in (7) then leads to $p=0.$ Therefore based on the form of the dynamic pressure $p$ in (5b), the equilibria satisfy the second-order nonlinear ODE (10).

For $\gamma>0$, we will show that the fourth-order equilibria satisfying the fourth-order ODE (9) co-exist with the second-order steady states of (5) in sufficiently large domains. This is a key difference between the steady state structure of our problem and that of classic conserved thin film models. The second-order steady states are solutions of the mass-conserving model and are independent of $\gamma$ while the fourth-order steady states depend on balances of the conserved and non-conserved fluxes.

We start with a brief discussion of the spatially uniform and second-order steady states. Since these solutions satisfy the second-order ODE (10), their properties are similar to those of the steady states of the model (5a) for $\gamma<0$ that have been studied in [30].

3 Spatially uniform steady states and their stability

The spatially uniform steady states of equation (5) are determined by the critical points of the potential $U(h)$ used in (7). For clarity, we explicitly write out the potential $U(h)$ as the integral of $\Pi(h)$ from (5c),

[TABLE]

Fig. 1 (left) shows a plot of $U(h)$ with $\bar{\mathcal{P}}=0.05$ which has a local minimum at $\bar{H}_{-}$ and a local maximum at $\bar{H}_{+}$, both of which satisfy (11).

This potential $U(h)$ has a unique inflection point at $h_{\text{peak}}=4/3$ (where $\Pi(h)$ has its global maximum), and on the semi-infinite range $h>h_{\mathrm{peak}}$, $U^{\prime\prime}(h)<0$ (see Fig. 1 (left)). This can be regarded as a degenerate case of the bi-stable double well potential used in Cahn-Hilliard models, which have a finite range of $h$ satisfying $U^{\prime\prime}(h)<0$, called the spinodal range.

The number of spatially uniform steady states of equation (5) depends on the value of $\bar{\mathcal{P}}$. Fig. 1 (right) depicts a plot of $\tilde{\Pi}(h)=h^{-3}-h^{-4}$ which has its unique maximum at $h_{\mathrm{peak}}$, yielding

[TABLE]

If $\bar{\mathcal{P}}<0$, there is one spatially uniform steady state $\bar{H}_{-}<h_{\text{peak}}$. For $0<\bar{\mathcal{P}}<\mathcal{P}_{\max}$, there are two uniform states $\bar{H}_{-}<h_{\text{peak}}<\bar{H}_{+}$; we will call this the critical range of the pressures $\bar{\mathcal{P}}$. As $\bar{\mathcal{P}}$ approaches $\mathcal{P}_{\max}$, the two states coincide at the double root $h=h_{\text{peak}}$. For $\bar{\mathcal{P}}>\mathcal{P}_{\max}$, there are no uniform steady states.

These ranges of $\bar{\mathcal{P}}$ are important for the analysis of the steady states and the dynamics of the model (5). For $\bar{\mathcal{P}}\leq\mathcal{P}_{\max}$, the uniform steady state $\bar{H}_{-}$ for $\bar{\mathcal{P}}\to 0$ is given by

[TABLE]

For $0<\bar{\mathcal{P}}\leq\mathcal{P}_{\max}$ the other uniform steady state $\bar{H}_{+}$ can be written as

[TABLE]

While $\bar{H}_{+}$ has a leading order dependence on $\bar{\mathcal{P}}$, $\bar{H}_{-}$ has a weaker dependence on $\bar{\mathcal{P}}$ with a saddle-node bifurcation occurring at $\bar{\mathcal{P}}=\mathcal{P}_{\max}$.

For $0<\bar{\mathcal{P}}\leq\mathcal{P}_{\max}$, the coexisting uniform steady states $\bar{H}_{+}$ and $\bar{H}_{-}$ give rise to novel and interesting solution structures and dynamics. Therefore, for the rest of this paper we focus on the critical range.

Now we consider the linear stability of a uniform steady state $\bar{H}$ with respect to an infinitesimal Fourier mode disturbance, $h(x,t)=\bar{H}+\delta e^{i2k\pi x/L}e^{\lambda t}$, where $k=0,1,2,3,\cdots$ and $\lambda$ describes the growth rate for perturbations starting from the initial amplitude $\delta\ll 1$. Expanding (5) about the uniform steady state $h=\bar{H}$ then gives the $O(\delta)$ equation

[TABLE]

where the first factor corresponds to an operator that includes the mobility functions from both conserved and non-conserved parts of (5a), and the second factor comes from the linearized pressure operator set by (5b) [30]. The dispersion relation (16) shows that the stability of the uniform steady states depends on both the parameter $\gamma$ and the domain size.

Fig. 2 gives a plot of the dispersion relation (16). For $\bar{H}=\bar{H}_{-}$, we have $\Pi^{\prime}(\bar{H}_{-})>0$ and the state $\bar{H}_{-}$ is long-wave unstable with respect to perturbations with $0\leq 2k\pi/L<[\gamma/M_{c}(\bar{H}_{-})]^{1/2}$. For $\bar{H}=\bar{H}_{+}$, since $\Pi^{\prime}(\bar{H}_{+})<0$, the state $\bar{H}_{+}$ is stable in the long-wave limit $L\to\infty$. The form of the curve depends on the parameter $\gamma$, and for small $\gamma$ (as in Fig. 2) the state $\bar{H}_{+}$ is unstable for $[\gamma/M_{c}(\bar{H}_{+})]^{1/2}<2k\pi/L<[-\Pi^{\prime}(\bar{H}_{+})]^{1/2}$. For $k=0$, (16) reduces to $\lambda=\gamma\Pi^{\prime}(\bar{H})$ corresponding to the spatially-uniform rate of change of mass per unit length in (8). We note that while the dispersion relation $\lambda(k)$ can be long-wave unstable, the influence of the fourth-order term always yields strong damping for large $k$. For small $k$ the form of the curves are similar to the type-II and type-III curves in [12]. These behaviors are comparable with dispersion relations for other volatile thin film models with gradient and non-gradient dynamics forms [60, 7].

The growth rate $\lambda$ changes sign at roots of the two factors in (16). In terms of the domain size $L$, these roots correspond to critical lengths at which the uniform state $\bar{H}$ changes stability. In particular, for $M_{c}(H)=H^{3}$ the first factor gives roots satisfying $\bar{H}^{3}(2k\pi/L)^{2}=\gamma$. We express these roots as $L=k\bar{\ell}^{\gamma}$ for $k=1,2,3,\cdots$ in terms of primary critical lengths, given by

[TABLE]

and

[TABLE]

Both ${\bar{\ell}_{-}}^{\,\gamma}$ and ${\bar{\ell}_{+}}^{\,\gamma}$ strongly depend on the parameter $\gamma$ and involve the interaction of both the conserved and non-conserved fluxes.

The roots of the second factor in (16) satisfy $(2k\pi/L)^{2}=-\Pi^{\prime}(\bar{H})$. Since $\Pi^{\prime}(\bar{H}_{-})>0$, the $\bar{H}_{-}$ state yields no real roots. But for $\bar{H}=\bar{H}_{+}$, $\Pi^{\prime}(\bar{H}_{+})<0$ and we obtain another family of critical lengths $L=k{\bar{\ell}_{+}}^{\,p}$ based on the linearized pressure, where the primary critical length is

[TABLE]

These critical lengths are also useful in characterizing the dependence of the stability of the uniform states $\bar{H}_{-}$ and $\bar{H}_{+}$ on the domain size $L$. In the limit $L\to 0$, both states are linearly stable to spatial perturbations. The state $\bar{H}_{-}$ is stable for $0<L<{\bar{\ell}_{-}}^{\,\gamma}$ and unstable for $L>{\bar{\ell}_{-}}^{\,\gamma}$, while the stability of the state $\bar{H}_{+}$ depends on the values of ${\bar{\ell}_{+}}^{\,p}$ and ${\bar{\ell}_{+}}^{\,\gamma}$.

At these critical lengths, $k{\bar{\ell}_{-}}^{\,\gamma}$, $k{\bar{\ell}_{+}}^{\,\gamma}$, and k${\bar{\ell}_{+}}^{\,p}$, nonuniform steady state solutions bifurcate from the uniform states. These will be the branches of second-order and fourth-order steady states. Moreover, these critical lengths give the only possible bifurcation points from the uniform states. Other apparent crossings between branches of nonuniform and uniform states at points not in this set are due to the projection of the full bifurcation structure onto the two-dimensional diagrams that we will be using, showing the average of the solution vs. the domain length. For an example of this in the next section, see Fig. 3 (left).

4 Second-order nonuniform steady states

Here we briefly describe the nonuniform second-order steady states. They are like steady states of broad families of conserved thin film equations that have been extensively studied in [37, 38] and elsewhere. The second-order ODE (10) can be written as a phase plane system

[TABLE]

Equilibrium points of this system correspond to the spatially uniform states $\bar{H}_{-}$, $\bar{H}_{+}$ satisfying (11). Linearizing the system around these equilibrium points using $h=\bar{H}+\delta e^{\sigma x}$ yields $\sigma=\pm\sqrt{\Pi^{\prime}(\bar{H}_{-})}$ for $\bar{H}=\bar{H}_{-}$, and $\sigma=\pm\mathrm{i}\sqrt{-\Pi^{\prime}(\bar{H}_{+})}$ for $\bar{H}=\bar{H}_{+}$. Hence $h=\bar{H}_{-}$ is a hyperbolic saddle point, and $h=\bar{H}_{+}$ is a center point. There is a homoclinic orbit through $\bar{H}_{-}$ with $\bar{H}_{-}\leq h\leq H_{\max}$ where $U(H_{\max})=U(\bar{H}_{-})$ and

[TABLE]

that represents the maximal size of a single droplet state on the whole real axis. In the phase plane this orbit encloses the continuous family of periodic solutions centered around $\bar{H}_{+}$.

For any fixed value of $\bar{\mathcal{P}}$ in the critical range $0<\bar{\mathcal{P}}\leq\mathcal{P}_{\max}$, selecting a minimum height $h_{\min}$ in the range $\bar{H}_{-}\leq h_{\min}\leq\bar{H}_{+}$ determines a periodic steady state solution $H_{A}(x)$ which has a maximum height $h_{\max}$ that satisfies $U(h_{\max})=U(h_{\min})$. Since we can write the first integral of (20) as $s^{2}=2(U(h)-U(h_{\min}))$, the period $\ell(h)$ of the steady state is given by

[TABLE]

In the limit case when $h_{\min}=\bar{H}_{-}$, we have the solitary droplet solution with the length $\ell\to\infty$.

Small amplitude second-order steady states bifurcate from the spatially-uniform steady state $\bar{H}_{+}$. Note that the second-order states are only related to the pressure (5b) and do not depend on the mobility, and their minimum period of oscillations is given by $L={\bar{\ell}_{+}}^{\,p}$ in (19). For $L\geq{\bar{\ell}_{+}}^{\,p}$, we define $h=H_{A}(x)$ as the principal second-order steady state with a single maximum on the domain $0\leq x\leq L$. With the average of the solution defined by

[TABLE]

Fig. 3 (left) shows the average of the solution, $\langle H_{A}\rangle$, for $\bar{\mathcal{P}}=0.05$ parametrized by $L$ bifurcating supercritically from $\bar{H}_{+}$ at $L={\bar{\ell}_{+}}^{\,p}$. Corresponding to the dots marked in Fig. 3 (left), typical profiles of steady states $H_{A}(x)$ with periods $L=30$ and $L=50$ are presented in Fig. 3 (right) using rescaled spatial variables. A bifurcation diagram similar to Fig. 3 (left) for steady states of a thin volatile film model can be found in [60].

To get an estimate of $\langle H_{A}\rangle$ for the limit $L\to\infty$, we use the result that if $h_{\max}\gg\bar{H}_{-}$ the core structure of $H_{A}(x)$ can be approximated by a truncated parabolic profile [5, 23] centered at $x=L/2$,

[TABLE]

Here $h_{\max}$ is the maximum of $H_{A}$ attained at $x=L/2$ and satisfies $\Pi(h_{\max})=d^{2}{H}_{A}/dx^{2}|_{x=L/2}$ from (10). For $L\to\infty$, $H_{A}$ approaches the homoclinic solution, and we can use $h_{\max}\approx H_{\max}$. As a consequence, a lower bound for the $\langle H_{A}\rangle$ branch as $L\to\infty$ is given by

[TABLE]

For $\bar{\mathcal{P}}=0.05$, this corresponds to the approximation $\langle H_{A}\rangle\approx 1.06+44.4/L$, which is plotted in Fig. 3 (left) in a dashed curve beneath the $\langle H_{A}\rangle$ branch. A plot of the truncated parabolic profile $\tilde{H}_{A}$ is also included in Fig. 3 (right) which shows a good agreement with the core structure of $H_{A}(x)$ for $L=50$.

For the domain size $L\geq 2{\bar{\ell}_{+}}^{\,p}$, multiple periodic nonuniform second-order steady states can coexist in the system, as families of second-order solutions $H_{A,k}(x)$ with $k$ periods bifurcate from $\bar{H}_{+}$ at $L=k{\bar{\ell}_{+}}^{\,p}$ for $k=2,3,\cdots$ We will use $\bar{\mathcal{P}}=0.05$ as a typical value for the parameter $\bar{\mathcal{P}}$ for the rest of the paper.

We will focus on the primary bifurcation point $L={\bar{\ell}_{+}}^{\,p}$. Using the rescaling $x=LX$, $H_{A}(x)=H(X)$ and perturbing the domain size $L$ around the bifurcation point ${\bar{\ell}_{+}}^{\,p}$ by a small parameter $\epsilon\ll 1$, $L={\bar{\ell}_{+}}^{\,p}+\epsilon$, we rewrite the ODE (10) as

[TABLE]

Following the local bifurcation analysis in [65, 5], we expand the solution in the neighborhood of ${\bar{\ell}_{+}}^{\,p}$ as

[TABLE]

where $\delta\ll 1$ is a small perturbation parameter, and we need to determine the relation between $\delta$ and $\epsilon$. Substituting (27) in (26) yields a harmonic oscillator equation at $O(\delta)$,

[TABLE]

which for ${\bar{\ell}_{+}}^{\,p}={2\pi}/\sqrt{-\Pi^{\prime}(\bar{H}_{+})}$ from (19) has the solution $H_{1}(X)=\mathcal{A}\cos(2\pi X),$ where $\mathcal{A}$ is the amplitude. The $O(\delta^{2})$ equation

[TABLE]

gives the solution

[TABLE]

In obtaining this form, we have assumed that $\epsilon\ll O(\delta)$ in order to exclude a resonant term from (29) that would have made it impossible to find a periodic solution.

Assuming that $\epsilon=O(\delta^{2})$, at $O(\delta^{3})$ the equation for $H_{3}(X)$ is $\hat{\mathcal{L}}H_{3}=\mathcal{R}_{3}(\bar{H}_{+},H_{1},H_{2})$, where $\mathcal{R}_{3}$ is the right-hand side inhomogeneous term. From the solvability condition for the $O(\delta^{3})$ equation, $\int_{0}^{1}\mathcal{R}_{3}H_{1}(x)~{}dx=0$, and the condition that the amplitude $\mathcal{A}$ of the $O(\delta)$ perturbation is real, we determine that the parameter $\epsilon$ is positive, $\epsilon>0$.

This indicates that a supercritical bifurcation occurs at $L={\bar{\ell}_{+}}^{\,p}$, and the amplitude $\mathcal{A}$ of the leading-order perturbation is

[TABLE]

The local structure of $\langle H_{A}\rangle$ near the bifurcation point $L={\bar{\ell}_{+}}^{\,p}$ is then given by

[TABLE]

Equivalently, the approximation $H_{A}(X)\sim\bar{H}_{+}+\mathcal{A}\cos(2\pi X)(L-{\bar{\ell}_{+}}^{\,p})^{1/2}$ holds in the limit $L\searrow{\bar{\ell}_{+}}^{\,p}$; plotting the extrema of the solution vs. $L$ would then show the local structure as a supercritical pitchfork bifurcation. Large-amplitude solutions of the strongly nonlinear equation (26) can also be constructed by matched asymptotic expansions [5].

Without loss of generality, the linear approximation (31) can be extended to

[TABLE]

for the branches of second-order states $H_{A,k}$ with multiple periods bifurcating from $\bar{H}_{+}$. These approximations in the neighborhood of bifurcation points $L={\bar{\ell}_{+}}^{\,p}$ and $L=2{\bar{\ell}_{+}}^{\,p}$ are plotted against numerical solutions in Fig. 3 (left). This bifurcation diagram will be revisited later for the regime where fourth-order nonuniform steady states also exist in the system.

The slopes of these linear approximations are positive since $\Pi^{\prime\prime}(\bar{H}_{+})>0$ and $\Pi^{\prime}(\bar{H}_{+})<0$ for $\bar{\mathcal{P}}=0.05$. The disjoining pressure has an inflection point at $h_{\mathrm{I}}=5/3$ where $\tilde{\Pi}^{\prime\prime}(h_{\mathrm{I}})=0$, with $\mathcal{P}_{\mathrm{I}}=54/625=0.0864$. For $\bar{\mathcal{P}}$ in the range $\mathcal{P}_{\mathrm{I}}<\bar{\mathcal{P}}<\mathcal{P}_{\max}$, the uniform state $\bar{H}_{+}$ satisfies $\Pi^{\prime\prime}(\bar{H}_{+})<0$ (see Fig. 1 (right)). This choice of $\bar{\mathcal{P}}$ will lead to the approximation (31) with a negative slope and the average $\langle H_{A,k}\rangle\leq\bar{H}_{+}$ for all second-order nonuniform solutions.

The results about the local structure of the second-order steady states are valid for any $\gamma\neq 0$. However, the stability of these states depend on the value of $\gamma$, which will be studied in Sec. 6. As a reminder, the appearance of each $H_{A,k}$ branch crossing $\bar{H}_{+}$ at some point to the right of $k{\bar{\ell}_{+}}^{\,p}$ for $\bar{\mathcal{P}}$ in the range $0<\bar{\mathcal{P}}<\mathcal{P}_{\mathrm{I}}$ is a consequence of the fact that Fig. 3 (left) is a two-dimensional projection of the full structure, as described in Sec. 3. These branches must pass by each other in parameter space since there are no bifurcation points of $\bar{H}_{+}$ at these values of $L$.

5 Fourth-order nonuniform steady states

Next we investigate the family of fourth-order nonuniform steady states that satisfy the fourth-order ODE (9) but not the second-order ODE (10). Equation (9) can be written as a fourth-order autonomous system

[TABLE]

with equilibrium points of the system corresponding to spatially uniform states, with $\Pi(h)=0$ (and $s=p=q=0$). As in the second-order case (20) these uniform states are $h=\bar{H}_{-}$ and $h=\bar{H}_{+}$. But now, linearization of the system around these equilibrium points with $h=\bar{H}+\delta e^{\sigma x}$ leads to $\sigma=\pm\sqrt{\Pi^{\prime}(\bar{H}_{-})}$, $\pm\mathrm{i}\sqrt{\gamma/\bar{H}_{-}^{3}}$ for $\bar{H}=\bar{H}_{-}$, and $\sigma=\pm\mathrm{i}\sqrt{-\Pi^{\prime}(\bar{H}_{+})}$, $\pm\mathrm{i}\sqrt{\gamma/\bar{H}_{+}^{3}}$ for $\bar{H}=\bar{H}_{+}$. Therefore $h=\bar{H}_{-}$ is a saddle-focus and $h=\bar{H}_{+}$ is a focus-focus type equilibrium point [35].

We note that following our choice of terminology from Section 2, for the fourth-order steady states, we will consider only solutions with nontrivial pressures, thus excluding the second-order and uniform steady states.

5.1 Primary bifurcations from uniform states

Small amplitude fourth-order nonuniform steady states bifurcate from both spatially-uniform states $\bar{H}_{-}$ and $\bar{H}_{+}$, and their minimum periods of oscillations are given by $L={\bar{\ell}_{-}}^{\,\gamma}$ in (18) and $L={\bar{\ell}_{+}}^{\,\gamma}$ in (17), respectively. Similar to the second-order case, we denote $H_{B}(x)$ as the principal branch of fourth-order states with period $L$, and denote the family of fourth-order states with $k$ periods as $H_{B,k}(x)$, $k=2,3,\cdots$,

For $\gamma=0.4$, Fig. 4 (left) depicts the first two branches of fourth-order states $H_{B}(x)$ and $H_{B,2}(x)$ in a bifurcation diagram parametrized by $\left\langle H\right\rangle$. These solutions start from the spatially-uniform state $\bar{H}_{-}$ at bifurcation points $L=k{\bar{\ell}_{-}}^{\,\gamma}$, and end at the uniform state $\bar{H}_{+}$ at $L=k{\bar{\ell}_{+}}^{\,\gamma}$. Typical plots of these states on the principal branch $H_{B}$ with periods $L=11.5$, $L=20$, and $L=32.8$ are shown in Fig. 4 (right) which correspond to the marked dots in Fig. 4 (left). Note that the $H_{A}$ nonuniform states are not included in this figure.

To determine the local bifurcation structure near a critical period $L=\ell_{c}$, we apply the approach used in Sec. 4 again. Applying the rescaling $x\to LX$ to (9) with the perturbation $L=\ell_{c}+\epsilon$ for $\epsilon\ll 1$ leads to the rescaled fourth-order ODE on $0\leq X\leq 1$

[TABLE]

In the neighborhood of $\ell_{c}$ we expand the steady state $H_{B}(X)$ as in (27), $H_{B}(X)=\bar{H}+\delta H_{1}(X)+\delta^{2}H_{2}(X)+\delta^{3}H_{3}(X)+O(\delta^{4})$ where $\delta\ll 1$. The $O(\delta)$ equation of the expansion of (34) leads to the critical period $\ell_{c}=2\pi\bar{H}^{3/2}/\sqrt{\gamma}$ as in (17) and (18).

Similar to the construction of the second-order solutions $H_{A}(x)$, the form of the fourth-order solutions for $\delta\to 0$ is given by

[TABLE]

where $\mathcal{A}$ is the amplitude of the leading-order perturbation to $\bar{H}$. We omit the details of calculating the $H_{2}(X)$ solution and move directly to the Fredholm solvability condition for the $O(\delta^{3})$ equation which determines that

[TABLE]

where

[TABLE]

where

[TABLE]

[TABLE]

Since $\mathcal{A}^{2}\geq 0$, the sign of $\mathcal{W}$ determines the relation between the small parameters $\epsilon$ and $\delta$. In particular, for $M_{c}(\bar{H})=\bar{H}^{3}$ and $\bar{H}=\bar{H}_{-}$, we have $\mathcal{W}>0$ for all $\gamma>0$, which yields $\epsilon=\delta^{2}$. Therefore, the average of the solution $\langle H_{B}\rangle$ near the supercritical bifurcation point $\ell_{c}={\bar{\ell}_{-}}^{\,\gamma}$ is approximated by

[TABLE]

For $\bar{H}=\bar{H}_{+}$, there are two cases depending on the value of the parameter $\gamma$ relative to functions of $\bar{H}_{+}$:

• •

For values of $\gamma$ yielding $\mathcal{W}<0$ then we take $\epsilon=-\delta\leq 0$ to indicate a subcritical bifurcation with

[TABLE]

• •

For values of $\gamma$ yielding $\mathcal{W}>0$ then $\epsilon=\delta\geq 0$ for a supercritical bifurcation with

[TABLE]

For $\gamma=0.4$ with $\mathcal{W}<0$, a comparison between $\langle H_{B}\rangle$ and their linear approximations near these bifurcation points are shown in Fig. 4 (left). We also note that the critical value $\gamma=-\Pi^{\prime}(\bar{H}_{+})\bar{H}_{+}^{3}/4$ is a root of the equation ${\bar{\ell}_{+}}^{\,\gamma}=2{\bar{\ell}_{+}}^{\,p}$ (and a zero of $\mathcal{W}(\gamma)$) where ${\bar{\ell}_{+}}^{\,p}$ and ${\bar{\ell}_{+}}^{\,\gamma}$ are defined in (19) and (17). Therefore we expect the local structure of the $H_{B}$ branch to change qualitatively as the critical period ${\bar{\ell}_{+}}^{\,\gamma}$ exceeds $2{\bar{\ell}_{+}}^{\,p}$.

5.2 Secondary bifurcations from nonuniform second-order states

The family of fourth-order states $H_{B,k}(x)$ may also undergo secondary bifurcations from the nonuniform second-order states $H_{A,k}(x)$. As an example, Fig. 5 shows the profiles and the $(H,H_{x})$ phase portraits of the second-order and fourth-order steady states, $H_{A,2}(x)$ and $H_{B}(x)$, that coexist at $L=61.7$. The profiles of the two steady states are close to each other, with the main peak of $H_{B}(x)$ being higher than that of $H_{A,2}(x)$, and the secondary peak of $H_{B}(x)$ is lower than that of $H_{A,2}(x)$, as in a period-doubling bifurcation. In this section, we use Floquet theory [25, 59, 10] to identify these secondary bifurcations.

We consider a fixed domain size $L$ near a critical period $\ell_{c}$ at which a fourth-order state $H_{B}$ bifurcates from a branch of second-order states $H_{A}$. From (9), the fourth-order state $H_{B}$ corresponds to a non-trivial pressure function $P_{B}(x)=\Pi(H_{B})-d^{2}H_{B}/dx^{2}$. In contrast, a second-order state $H_{A}$ corresponds to a trivial pressure $p\equiv 0$. We expand the fourth-order state and its corresponding pressure $(H_{B},P_{B})$ around the principal second-order state $H_{A}(x)$ and its corresponding zero pressure using

[TABLE]

where the small parameter $\delta\ll 1$ is set by the distance $|L-\ell_{c}|$. Substituting these expansions into (9) and using (10) we obtain the $O(\delta)$ linearized problem

[TABLE]

A secondary bifurcation of $H_{B}$ from $H_{A}$ will be identified if (40) has a non-trivial $L$-periodic solution $(\tilde{H},\tilde{P})$, $\tilde{H}\not\equiv 0$ and $\tilde{P}\not\equiv 0$.

For a given $L$-periodic $H_{A}$ state, the problem (40) is a fourth-order differential equation with periodic coefficients set by $H_{A}$. While a periodic solution to (40) needs to satisfy the full system, we can get necessary (but not sufficient) conditions for the existence of such periodic solutions by first looking at the $2\times 2$ system for $\mathbf{Y}(x)=[\tilde{P}(x),\tilde{Q}(x)]^{T}$ from the first equation in (40),

[TABLE]

We numerically solve this system as an initial value problem, and obtain the principal fundamental matrix solution $\mathbf{\Phi}(x)$ which satisfies

[TABLE]

Here we use a predetermined $H_{A}$ with its peak located at $x=L/2$ (see Fig. 3 (right)), and $\mathbf{I}_{2}$ denotes the $2\times 2$ identity matrix. Evaluating the matrix solution $\mathbf{\Phi}(x)$ at $x=L$, we obtain the monodromy matrix $\mathbf{B}_{2}=\mathbf{\Phi}(L)$ whose eigenvalues are characteristic multipliers of equation (41) [59, Section 3.6]. Based on Floquet theory [25, Chapter II section 1], if the matrix $\mathbf{B}_{2}$ has an eigenvalue $\rho$ that satisfies $\rho^{k}=1$ for a positive integer $k$, then there exists a $kL$-periodic solution $\mathbf{Y}$ of ODE (41). In particular, when $\mathbf{B}$ has an eigenvalue $\rho=1$, it corresponds to an $L$-periodic nonuniform solution $\mathbf{Y}$ of (41).

Fig. 6 (right) shows the dependence of characteristic multipliers $\rho$ on the domain size $L$ for $\gamma=0.4$. The “egg-shaped” portions of the plot correspond to pairs of real multipliers, and the remaining portions represent complex conjugates pairs of multipliers $\rho,\bar{\rho}$ satisfying $|\rho|=1$. As $L$ increases the structure of the set of multiplier values appears to converge to a repeating pattern. Our numerical study also shows that at the critical domain sizes corresponding to a characteristic multiplier $\rho=1$, the monodromy matrix $\mathbf{B}_{2}$ takes the form

[TABLE]

where the constants $a,b\neq 0$, indicating that the multiplier $\rho=1$ is a double eigenvalue of geometric multiplicity one.

While the Floquet analysis for the second-order system (41) identifies multiple critical domain sizes $L$ with $\rho=1$, we will show that only the first critical $L$ corresponds to a secondary bifurcation of $H_{A}$. To exclude the other spurious critical $L$, we apply Floquet analysis to the full system in (40), and consider the fourth-order differential equation system for $\mathbf{Z}(x)=[\tilde{H}(x),\tilde{W}(x),\tilde{P}(x),\tilde{Q}(x)]^{T}$,

[TABLE]

Again the principal fundamental solution matrix $\mathbf{\Psi}$ of the system (44) satisfies

[TABLE]

where the solution matrix $\mathbf{\Psi}=\left[\psi_{1},\psi_{2},\psi_{3},\psi_{4}\right]$, with the $i$-th column vector $\psi_{i}(x)=[\tilde{H}_{i}(x),\tilde{W}_{i}(x),\tilde{P}_{i}(x),\tilde{Q}_{i}(x)]^{T}$, and $\mathbf{I}_{4}$ denotes the $4\times 4$ identity matrix. Since the differential equations for $\tilde{P}$ and $\tilde{Q}$ do not depend on $\tilde{H}$ or $\tilde{W}$, the monodromy matrix $\mathbf{B}_{4}=\mathbf{\Psi}(L)$ of equation (44) can be written as a block matrix of four $2\times 2$ matrices,

[TABLE]

where

[TABLE]

where the constant $c\neq 0$. At the critical domain sizes identified in Fig. 6 (right), the matrix $\mathbf{B}_{2}$ is given by (43). We numerically calculate the eigenvalues of the monodromy matrix $\mathbf{B}_{4}$ and plot the corresponding characteristic multipliers $\rho$ in Fig. 6 (left). In addition to the characteristic multipliers that have appeared in Fig. 6 (right) for the second-order system (41), there exists another double multiplier $\rho=1$ of geometric multiplicity $1$ for all $L$. This additional $\rho=1$ corresponds to the translational mode $\tilde{H}=dH_{A}/dx$ that lies on the branch of second-order states $H_{A}$ and does not yield any non-trivial $\tilde{P}$ solutions.

Since ${\psi}_{1}(x)$ and ${\psi}_{2}(x)$ only contain trivial $\tilde{P}$ solutions, in order to construct a periodic solution to (44) with a non-trivial $\tilde{P}$, one needs to use a linear combination of ${\psi}_{3}(x)$ and ${\psi}_{4}(x)$ which satisfy ${\psi}_{3}(0)=[0,0,1,0]^{T}$ and ${\psi}_{4}(0)=[0,0,0,1]^{T}$, $\mathbf{Z}(x)=c_{3}{\psi}_{3}(x)+c_{4}{\psi}_{4}(x)$. The periodicity of $\mathbf{Z}$ requires $\mathbf{Z}(0)=\mathbf{Z}(L)$, which leads to

[TABLE]

Note that at the critical periods, $\mathbf{B}_{2}$ takes the form (43), therefore we have $\tilde{P}_{3}(L)=\tilde{Q}_{4}(L)=1$, and (48) leads to $c_{4}\tilde{P}_{4}(L)=c_{3}\tilde{Q}_{3}(L)=0$. Moreover, since $(\tilde{P}_{4}(L),\tilde{Q}_{3}(L))=(0,a)$ or $(b,0)$ with $a,b\neq 0$, the coefficients $c_{3}$ and $c_{4}$ satisfy $c_{3}=0$ or $c_{4}=0$. This indicates that the periodic solution $\tilde{H}$ is given by a multiple of $\tilde{H}_{3}(x)$ or $\tilde{H_{4}}(x)$.

Therefore, it suffices to check whether $\tilde{H_{3}}(L)=0$ or $\tilde{H_{4}}(L)=0$ to determine if a candidate critical period of the system (41) indeed corresponds to a periodic solution $\tilde{H}$ of (44) with a non-trivial $\tilde{P}$. The plot of $\tilde{H}_{3}(L)\tilde{H}_{4}(L)$ as a function of $L$ in Fig. 6 (right) shows that only the smallest critical period $L=\ell_{AB}^{\,(1)}$ associated with $\rho=1$ corresponds to such a periodic solution. This domain size yields a symmetry-preserving perturbation solution $\tilde{H}$ with respect to the reflectional symmetry of the second-order state $H_{A}(x)$ about $x=L/2$, and corresponds to a secondary bifurcation from the principal second-order state $H_{A}$ to the principal fourth-order state $H_{B}$.

For other characteristic multipliers satisfying $\rho^{k}=1$ with higher values of $k$, more perturbation solutions $(\tilde{H},\tilde{P})$ with period $kL$ exist in the system. They can lead to more interesting symmetry-preserving and symmetry-breaking secondary bifurcations, similar to the bifurcation of $H_{B}$ from $H_{A,2}$ observed in Fig. 5. We will not attempt an exhaustive investigation of these bifurcations in this paper.

5.3 Overall structure of the solution branches

So far we have worked out some local properties of the $H_{A,k}$ and $H_{B,k}$ branches. We will now consider some overall aspects of how these branches interact with each other.

Recall that the $H_{A,k}$ branches always start from a supercritical bifurcation at $L=k{\bar{\ell}_{+}}^{\,p}$ from the uniform state $\bar{H}_{+}$ and have $\langle H_{A}\rangle$ asymptotically approaching $\bar{H}_{-}$ as $L\to\infty$ (see Fig. 3 (left)). These bifurcation points $L=k{\bar{\ell}_{+}}^{\,p}$ are set by $\bar{\mathcal{P}}$ and are independent of $\gamma$. In contrast, the bifurcation points $L=k\bar{\ell}^{\,\gamma}_{\pm}$ at which the $H_{B,k}$ branches bifurcate from $\bar{H}_{\pm}$ depend strongly on $\gamma$. Therefore, different configurations of the $H_{A,k}$ and $H_{B,k}$ branches can be expected for different values of $\gamma$.

For simplicity we focus on the principal branch of fourth-order states $H_{B}(x)$ (we will briefly mention the generalization for $H_{B,k}$ with $k=2,3,\cdots$ later) and discuss the following two basic cases:

• •

Case I: Primary-Primary branches – branches starting and ending from primary bifurcations with constant states $\bar{H}_{\pm}$ at $\bar{\ell}^{\,\gamma}_{\pm}$.

• •

Case II: Primary-Secondary branches – branches starting at a primary bifurcation from $\bar{\ell}^{\,\gamma}$, and ending at a secondary bifurcation with a $H_{A,k}$ branch.

5.3.1 Case I: Primary-Primary branches

We start with the simplest case of an $H_{B}$ branch that starts and ends from primary bifurcations, and does not intersect with any $H_{A}$ branches. To ensure that there are no intersections, $\gamma$ must satisfy ${\bar{\ell}_{+}}^{\,\gamma}<{\bar{\ell}_{+}}^{\,p}$. Having $\gamma$ in this range was the first case discussed in Sec. 5.1. An example of such an $H_{B}$ branch is shown in Fig. 7 (left) with $\gamma=0.6$, where the $H_{B}$ branch starts from a supercritical bifurcation at $L={\bar{\ell}_{-}}^{\,\gamma}$, and ends at a subcritical bifurcation at $L={\bar{\ell}_{+}}^{\,\gamma}$. The $H_{B}(x)$ solutions on this branch all have a single local maxima.

A primary-primary $H_{B}$ branch may also undergo a secondary bifurcation as a transcritical crossing with the principal branch $H_{A}$. This regime corresponds to the value of $\gamma$ satisfying ${\bar{\ell}_{+}}^{\,p}<{\bar{\ell}_{+}}^{\,\gamma}<2{\bar{\ell}_{+}}^{\,p}$, so that $H_{B}$ intersects the principal second-order branch $H_{A}$, but not other second-order branches. Using the local properties of the $H_{B}$ branch discussed in section 5.1, having $\gamma$ satisfy the upper bound ${\bar{\ell}_{+}}^{\,\gamma}=2{\bar{\ell}_{+}}^{\,p}$ corresponds to the primary bifurcation at $L={\bar{\ell}_{+}}^{\,\gamma}$ changing from a subcritical bifurcation to a supercritical when as ${\bar{\ell}_{+}}^{\,\gamma}$ exceeds $2{\bar{\ell}_{+}}^{\,p}$. Fig. 7 (right) shows an example of this case with $\gamma=0.4$, where the $H_{B}$ branch starts from $\bar{H}_{-}$, ends at $\bar{H}_{+}$, and crosses the $H_{A}$ branch at a transcritical secondary bifurcation point $L=\ell_{AB}^{\,(1)}$. This secondary bifurcation at $L=\ell_{AB}^{\,(1)}$ has been identified in section 5.2. The fourth-order states on this $H_{B}$ branch also have a unique maximum (see Fig. 4 (right)) and coincide with the second-order state at $L=\ell_{AB}^{\,(1)}$.

5.3.2 Case II: Primary-Secondary Branches

More interesting primary-secondary $H_{B}$ branches emerge for small $\gamma$ values satisfying ${\bar{\ell}_{+}}^{\,\gamma}>2{\bar{\ell}_{+}}^{\,p}$. One such example for $\gamma=0.05$ is shown in Fig. 8 (top). In this case, the overall $H_{B}$ branch takes the form of two disconnected primary-secondary branches, where the first branch starts from the uniform state $\bar{H}_{-}$ and ends at the $H_{A,2}$ branch, and the second branch connects the uniform state $\bar{H}_{+}$ and the $H_{A,3}$ branch.

In particular, the first branch starts at the primary bifurcation at $L={\bar{\ell}_{-}}^{\,\gamma}$, ends at the secondary bifurcation at $L=\ell_{AB}^{\,(2)}$, and crosses the principal second-order branch $H_{A}$ at a transcritical bifurcation point $L=\ell_{AB}^{\,(1)}$ similar to the primary-primary branch shown in Fig. 7 (right). Note that the fourth-order states $H_{B}$ on this branch connect the $H_{A}$ states with a single peak and the $H_{A,2}$ states with two peaks. That is, the properties of $H_{B}$ change qualitatively as the branch is followed from $L={\bar{\ell}_{-}}^{\,\gamma}$ to $L=\ell_{AB}^{\,(2)}$. Fig. 9 (left) depicts typical profiles of $H_{B}(x)$ on a closeup view of the bifurcation diagram from Fig. 8 (top). It shows that a secondary peak develops in the $H_{B}$ profiles moving along this branch towards $L=\ell_{AB}^{\,(2)}$, and the two peaks become identical when $L=\ell_{AB}^{\,(2)}$ is reached (also see Fig. 5). This observation indicates that the $H_{B}$ solutions from the pitchfork bifurcation at $L=\ell_{AB}^{\,(2)}$ break the discrete $L/2$ translation symmetry of the $H_{A,2}$ solutions. The point where the secondary peak first develops is above $\ell_{AB}^{\,(1)}$ and is distinct from that point; it has not been labeled in the figures. The fold point $L=\ell_{B}^{\mathrm{fold},1}$ on this branch does not change the qualitative structure of the profiles, but does influence the stability of the solutions as a saddle-node bifurcation.

The second $H_{B}$ branch connects the supercritical primary bifurcation at $L={\bar{\ell}_{+}}^{\,\gamma}$ and the secondary bifurcation at $L=\ell_{AB}^{\,(3)}$. No transcritical secondary bifurcations are observed on this branch. Typical profiles of $H_{B}(x)$ on this branch are plotted in Fig. 9 (right). At the pitchfork bifurcation $L=\ell_{AB}^{\,(3)}$, where $H_{B}$ bifurcates from $H_{A,3}$, the discrete $L/3$ translation symmetry of $H_{A,3}$ is broken. This bifurcation also breaks most of the reflection symmetries about the extrema of $H_{A,3}$. The $H_{B}$ solutions have one large primary peak and two smaller secondary peaks that shrink as $L$ approaches ${\bar{\ell}_{+}}^{\,\gamma}$, and transition to single-peak solutions as the branch connects to $\bar{H}_{+}$ at the pitchfork bifurcation $L={\bar{\ell}_{+}}^{\,\gamma}$. Similar to the first $H_{B}$ branch, the fold point $L=\ell_{B}^{\mathrm{fold},2}>{\bar{\ell}_{+}}^{\,\gamma}$ on this branch impacts the stability of the solutions, but does not change their symmetry. The transition point for solutions changing from three peaks to one peak has not been labeled but was found to occur after the fold point.

To better characterize the coexisting states in this regime, we also plot the extremal values of the pressures $p$ in Fig. 8 (middle) corresponding to the steady states from Fig. 8 (top). The $p\equiv 0$ trivial branch corresponds to the spatially uniform steady states $\bar{H}_{-}$, $\bar{H}_{+}$, and the second-order steady states $H_{A,k}(x)$. The branches with $p\not\equiv 0$ exist for ${\bar{\ell}_{-}}^{\,\gamma}<L<\ell_{B}^{\mathrm{fold},1}$ and $\ell_{AB}^{\,(3)}<L<\ell_{B}^{\mathrm{fold},2}$, and represent the two branches of fourth-order steady states $H_{B}$ that undergo primary bifurcations from $\bar{\ell}^{\,\gamma}_{\pm}$ and secondary bifurcations from $\ell_{AB}^{\,(k)}$ where $k=1,2,3$. This is an alternative way of detecting transcritical bifurcation points involving $H_{B}$ that does not require doing the full linear stability analysis.

Using the approach developed in Sec. 5.2, we also identify these secondary bifurcations $\ell_{AB}^{\,(k)}$. Fig. 8 (bottom) shows the real part of the characteristic multipliers $\rho$ for the second-order system (41) for $\gamma=0.05$. The three sections correspond to a partial structure of the multipliers obtained for the states $H_{A}$, $H_{A,2}$, and $H_{A,3}$, respectively. The full structures for these states are similar to the one presented in Fig. 6 (right) for $\gamma=0.4$. Then for $k=1,2,3$, we solve the fourth-order system (44) using $H_{A,k}$ as the base state, and obtain the corresponding $\tilde{H}_{3}(L)\tilde{H}_{4}(L)$ as functions of $L$. Recall that an admissible critical domain size $L=\ell_{AB}^{\,(k)}$ for a secondary bifurcation from $H_{A,k}$ of period $\ell_{AB}^{\,(k)}/k$ needs to satisfy the following two conditions for the base state $H_{A,k}$:

• •

The fourth-order system (44) gives $\tilde{H}_{3}(L)\tilde{H}_{4}(L)=0$;

• •

There exists a characteristic multiplier $\rho=1$ for the second-order system (41).

The rescaled $\tilde{H}_{3}(L)\tilde{H}_{4}(L)$ curves for $k=1,2,3$ in Fig. 8 (bottom) show that their values reach zero, along with the corresponding characteristic multiplier $\rho=1$, at the bifurcation points $L=\ell_{AB}^{\,(k)}$ identified in the top and middle panels. It is also worth mentioning that the bifurcation points $\ell_{AB}^{\,(k)}$ with $k>1$ can also be folded back into the characteristic multiplier plots of the principal state $H_{A}$ instead of $H_{A,k}$, since a multiplier satisfying $\rho^{k}=1$ at $L=\ell_{AB}^{\,(k)}/k$ implies the existence of a periodic solution for $L=\ell_{AB}^{\,(k)}$.

There are very likely more complicated cases involving higher order bifurcations, yielding period doubling of $H_{B}$ solutions, and an extensive investigation of all the interesting structures is beyond the scope of this study. The rest of the paper will focus on the stability and dynamical behavior of the principal branches of the solutions.

6 Stability analysis of nonuniform steady states

To examine the stability of the steady states, we consider an $L$-periodic positive steady state $H(x)$ over the domain $0\leq x\leq L$ and perturb it by setting $h(x,t)\sim H(x)+\delta\Psi(x)e^{\lambda t}$, where $\delta\ll 1$ and $\Psi(x)$ is also $L$-periodic. We linearize equation (5) around the base state $H(x)$, a solution of (9), and obtain the $O(\delta)$ equation

[TABLE]

where the linear operator $\mathscr{L}$ is given by

[TABLE]

If any eigenvalue $\lambda$ of problem (49) has a positive real part, then the steady state $H(x)$ is unstable. Since the second-order steady states satisfy $p(H_{A})\equiv 0$, the linear operator $\mathscr{L}$ in (50) reduces to a simpler operator $\mathscr{L}_{A}$ for $H=H_{A}(x)$

[TABLE]

For the case $\gamma\leq 0$ the instability of the second-order nonuniform steady states has been studied in [30]. Similar stability analysis for second-order steady states in related equations has been performed in [36, 17]. Here we will extend this study to the case $\gamma>0$ which reveals critical domain sizes that give rise to limit cycles and other interesting dynamics. For simplicity, for the rest of the paper we only focus on the case where only a single-period solution fits in the domain, but these results can be easily extended to states with multiple periods. For instance, symmetry breaking, droplet merging, and leveling dynamics arising from a second-order state with two periods have been numerically investigated in [30].

First we numerically solve the eigenproblem (49) for both nonuniform steady states $H_{A}$ and $H_{B}$ for a range of domain sizes $L$. For all $L$, a translational mode $\lambda^{T}=0$ with $\Psi=H^{\prime}(x)$ is present on the periodic domain. The dependence of the dominant eigenvalues on $L$ is plotted in Fig. 10, showing that each state has a dominant real eigenvalue crossing zero at the bifurcation point $\ell_{AB}^{\,(1)}$. This indicates that $\ell_{AB}^{(1)}$ labelled in Fig. 8 is a transcritical bifurcation point. In addition, both plots show pairs of real eigenvalues transitioning to become complex conjugate pairs that yield secondary Hopf bifurcations that have not been labelled in Fig. 8. Further, the fold points shown in Fig. 8 indicate saddle-node bifurcations, where a real eigenvalue crosses zero. The dependence of the dominant eigenvalues on the system parameter $\gamma$ is illustrated in Fig. 11 for a second-order state $H_{A}$ with a fixed domain size $L$.

These numerical results indicate that the stability of the steady states depends strongly on both $L$ and $\gamma$. Therefore, stability in plots of the solution branches parametrized by the domain size $L$ can change significantly when $\gamma$ is changed. For example, Fig. 11 shows that the stability of the $H_{A}$ solution with a fixed domain size undergoes qualitative changes as $\gamma$ varies, despite the fact that $H_{A}$ itself is independent of $\gamma$.

Next we study the stability of $H_{A}(x)$ in the limit of weak non-conserved flux for $\gamma\to 0$. Since the stability of these states is closely related to the stability of equilibria of the mass-conserving equation (5) with $\gamma=0$, we present below a brief overview of the $\gamma=0$ case.

For $\gamma=0$, the linear operator $\mathscr{L}_{A}$ in (51) is self-adjoint with respect to a weighted $H^{-1}$ norm [65, 4]. Hence from spectral theory for $\gamma=0$ the state $H_{A}$ on a compact domain is associated with a spectrum that is real and discrete. Moreover, it has been shown in [30] that for $\gamma=0$ the state $H_{A}$ changes stability at a critical $L=\ell^{*}$, where $\ell^{*}$ depends on the parameter $\bar{\mathcal{P}}$. That is, for ${\bar{\ell}_{+}}^{\,p}<L<\ell^{*}$ the state $H_{A}$ is unstable, and for $L>\ell^{*}$ it is stable.

In the non-conserved case with small $\gamma>0$, the stability of the state $H_{A}$ depends on both $L$ and $\gamma$. For a fixed period $L=29.2$ in the range ${\bar{\ell}_{+}}^{\,p}<L<\ell^{*}$, we numerically solve the eigenproblem (49) with $\mathscr{L}=\mathscr{L}_{A}$ for the state $H_{A}(x)$. The dominant eigenvalues of $H_{A}$ parametrized by the parameter $\gamma$ are plotted in Fig. 11. The figure shows that in addition to the translational eigenvalue $\lambda^{T}=0$, for small positive $\gamma>0$, $H_{A}(x)$ has two unstable eigenvalues $\lambda^{P},\lambda^{V}>0$. In the limit $\gamma\to 0^{+}$, $\lambda^{P}$ approaches zero, and $\lambda^{V}$ approaches a positive value given by the unstable eigenvalue of $H_{A}$ in the $\gamma=0$ case. At $\gamma_{c}^{-}\approx 0.0059$, the two eigenvalues merge together and become a pair of complex conjugate eigenvalues $\lambda^{C},\bar{\lambda}^{C}$ with positive real parts. At the critical value $\gamma_{c}^{+}\approx 0.041$, the complex conjugate pair cross the imaginary axis and a Hopf bifurcation occurs. The real part of the complex conjugate pair is positive for $\gamma_{c}^{-}\leq\gamma\leq\gamma_{c}^{+}$, and is negative for $\gamma_{c}^{+}\leq\gamma\leq\gamma_{c}^{*}$, where $\gamma_{c}^{*}\approx 0.2$. Therefore this suggests a limit cycle bifurcating from the steady state $H_{A}(x)$ at the Hopf bifurcation point $\gamma=\gamma_{c}^{+}$.

In summary, on the range ${\bar{\ell}_{+}}^{\,p}<L<\ell^{*}$ Hopf bifurcations occur at $\gamma=\gamma_{c}^{+}(L)$ with $\Re(\lambda^{C}(\gamma))=0$. To express this in terms of critical domain sizes, we define two critical periods $\ell_{AC}^{\,(1)}$ and $\ell_{AC}^{\,(2)}$ as functions of $\gamma$ such that $\ell_{AC}^{\,(1,2)}(\gamma_{c}^{+})=L$. The two curves $\ell_{AC}^{\,(1)}(\gamma)$ and $\ell_{AC}^{\,(2)}(\gamma)$ merge at $\gamma^{*}\approx 0.041$. The region bounded by $\ell_{AC}^{\,(1,2)}(\gamma)$ on $0<\gamma<\gamma^{*}$ is shown in the bifurcation diagram parametrized by $(\gamma,L)$ in Fig. 12.

7 Limit cycle dynamics

We will show that stable limit cycles occur in the finite region marked as $H_{C}$ above in Fig. 12. The dynamics of the observed stable limit cycles can be described by plotting some global properties, such as the average of the solution $\langle h\rangle$ or the extreme values of the solution profile, $h_{\min}=\min_{x}h$ and $h_{\max}=\max_{x}h$, as functions of time. To understand the significance of the stability of the second-order steady states in this regime, we perform PDE simulations starting from initial conditions of the form

[TABLE]

where $\Psi(x)$ is an unstable eigenmode of $H=H_{A}$ associated with an eigenvalue $\lambda$ with $\Re(\lambda)>0$, and $\|\Psi\|_{2}=1$, $\Psi(0)\geq 0$ and $|\delta|\leq 1$. At the Hopf bifurcation $\gamma=\gamma_{c}^{+}$ a stable limit cycle bifurcates from $H_{A}$ in the form

[TABLE]

where the period of the limit cycle oscillations is determined by the parameter $\gamma$.

Fig. 13 (left) shows the evolution of two such limit cycle simulations parametrized by $\langle h\rangle$ with $\gamma=0.02$ and $\gamma=0.006$, where both $\gamma$ values satisfy $\gamma<\gamma_{c}^{+}$. It is observed that the period of the limit cycle with $\gamma=0.02$ is significantly smaller than the one with $\gamma=0.006$. For $\gamma^{+}_{c}<\gamma<\gamma_{c}^{*}$, the real part of the complex conjugate pair $\lambda^{C},\bar{\lambda}^{C}$ becomes negative (see Fig. 11), and the corresponding simulation with $\gamma=0.05$ in Fig. 13 (right) shows that the initial spatial perturbations exhibit oscillatory decay over time. For $\gamma>\gamma_{c}^{*}$, the dominant eigenvalues of the steady state $H_{A}(x)$ are real and negative. The corresponding PDE simulation in Fig. 13 (right) with $\gamma=0.21>\gamma_{c}^{*}$ shows that the dynamic solution converges monotonically to the stable steady state $H_{A}(x)$.

Using numerical continuation with respect to the parameter $\gamma$, we also identify a continuous family of stable limit cycles for $0<\gamma<\gamma_{c}^{+}$. Fig. 14 shows the dependence on the parameter $\gamma$ of two properties of the limit cycles, the profile magnitude $h_{\max}-h_{\min}$ and the average of the solution $\langle h\rangle$. As $\gamma\to 0$, the magnitude of the limit cycle oscillations approaches a positive limit magnitude with an increasing period. For $\gamma>\gamma_{c}^{+}$ the stable limit cycle vanishes, and the dynamic solution converges to the stable steady state $H_{A}(x)$ (see Fig. 13 (right)).

Typical relaxation ocillator-type limit cycles in fast/slow dynamical systems [45] consist of four stages (shown in Fig. 15), where the slow-time solution evolves on a slow manifold in stages $(1)$ and $(3)$, and alternates with fast dynamics occurring in stages $(2)$ and $(4)$. We begin by determining the slow-time solution. By rescaling the original PDE (5) using

[TABLE]

we get the full PDE system for $(H,P)$

[TABLE]

Using regular perturbation expansions for both $H$ and $P$,

[TABLE]

[TABLE]

and collecting the $O(1)$ and $O(\gamma)$ terms after substituting the expansions in (54), we obtain the leading order slow system for $\gamma\to 0$

[TABLE]

and the order $O(\gamma)$ system

[TABLE]

The leading order system (55) gives a quasi-steady solution, and imposing periodic boundary conditions leads to a spatially-uniform leading order pressure term $P_{0}\equiv P_{0}(T)$. With different values of $P_{0}$ there are coexisting solutions to the second-order ODE,

[TABLE]

which have been studied extensively in [30]. In addition to the constant solution $\bar{H}$ that satisfies $P_{0}=\Pi(\bar{H})$, a family of nonuniform solutions $H_{0}(x)$ also exists for a range of $P_{0}$. Fig. 16 (left) shows a bifurcation diagram of these solutions parametrized by their average $\langle H\rangle$ with the domain size $L=29.2$. The branch of nonuniform solutions $H_{0}(x)$ bifurcates from the constant states $\bar{H}$ at a critical $P_{0}=P_{c}$, which corresponds to a critical solution average $P_{c}=\Pi(\langle H_{c}\rangle)$. Moreover, there exists a critical $P^{*}>0$ such that for $P_{0}>P^{*}$ the average of the nonuniform solution $\langle H_{0}\rangle$ monotonically increases as $P_{0}$ decreases. While for $P_{c}<P_{0}<P^{*}$ the nonuniform solution has a decreasing average $\langle H_{0}\rangle$ as $P_{0}$ decreases. As a result, for any average in the range $\langle H_{c}\rangle<\langle H\rangle<\langle H^{*}\rangle$, there exist two nonuniform solutions $H_{0}$ satisfying (57) with different values of $P_{0}$. With the solution average satisfying $\langle H\rangle<\langle H_{c}\rangle$, only one nonuniform solution $H_{0}$ coexists with the constant solution $\bar{H}$.

Integrating both sides of the $O(\gamma)$ equation (56)1 over the periodic domain, one obtains

[TABLE]

This allows us to qualitatively determine the dynamics of the average $\langle H_{0}\rangle$ of the solution on the slow manifold, as $\langle H_{0}(T)\rangle$ is increasing in time for $P_{0}>0$, and is decreasing for $P_{0}<0$. Restricting the solution to stay on the slow manifold, then starting from any quasi-steady solution $H_{0}(x)$, with $P_{0}<0$ the solution will be driven to the critical point $(\bar{H}_{+},0)$ in Fig. 16 (left), and for $P_{0}>0$ the solution evolves towards the other critical point $(\langle H^{*}\rangle,P^{*})$. Fig. 16 (right) plots the limit cycle dynamics shown in Fig. 15 parametrized by the average pressure $\langle p\rangle=L^{-1}\int_{0}^{L}p(h)~{}dx$ and the average $\langle h\rangle$ of the dynamic solution. It shows that in the slow stages $(1)$ and $(3)$ the PDE solution evolves on the slow manifold in the direction based on the sign of $\langle p\rangle$ as described above.

When the two critical points $(\bar{H}_{+},0)$ and $(\langle H^{*}\rangle,P^{*})$ are approached, the dynamic solution jumps off the slow manifold and is governed by the fast-time problem, which can be described using the original PDE (5a) and (5b). For $\gamma\to 0$ we apply the regular perturbation expansions for both $h$ and $p$,

[TABLE]

[TABLE]

and obtain the leading order fast problem

[TABLE]

which is identical to the mass-conserving thin film equation (4). This indicates that for the fast dynamics in the limit cycle, the leading order average of the solution $\langle h_{0}\rangle$ is constant in time when $p_{0}$ is evolving, which explains the motions in the fast stages (2 and 4) in Fig. 16 (right).

8 Finite-time singularity formation

Next, we illustrate how the linear instability of a fourth-order steady state can lead to another fundamental mode of the dynamics of (5) – the formation of a finite-time rupture singularity [15]. Fig. 17 (left) depicts a typical numerical simulation of (5) for $\gamma=0.05$ and $L=50$ with initial data $h_{0}(x)=H_{B}(x)+0.03\Psi_{B}(x)$, where $\Psi_{B}(x)$ is an unstable eigenmode associated with the complex conjugate pair eigenvalues $\lambda_{B}\approx 0.0017\pm\mathrm{i}0.0013$. After initial transients, the minimum decreases and approaches $h\to 0$ at an isolated point $x_{c}$ at a finite critical time, $t_{c}$. Since $\Pi(h)\to-\infty$ as $h\to 0$, this is a singularity of the dynamics and the solution cannot be continued beyond the critical time. The simulation also suggests that the singularity formation is localized, that is, away from $x_{c}$, the solution remains smoothly and slowly evolving even as properties become singular in a small neighborhood of $x_{c}$ as $t\to t_{c}$.

In [29] it was shown that first-kind self-similar dynamics lead to this singularity, and we briefly summarize key points here. Let $\tau$ be the time remaining until rupture, $\tau=t_{c}-t$, we can seek a focusing self-similar solution for rupture at $x_{c}$ of the form

[TABLE]

where the positive scaling exponents $\alpha,\beta$ must be determined. Substituting this ansatz into (5) will not balance all terms in the PDE, but by considering the limit $\tau\to 0$ with $H=O(1)$, we can seek an asymptotically self-similar solution that satisfies a leading-order dominant balance of terms. Since the singularity cannot form without the non-conserved term [30], we know the balance must involve part of the $\gamma p$ term balancing the rate of change and other terms as $\tau\to 0$; at leading order this determines the scaling exponents as $\alpha=1/5$ and $\beta=3/10$ and reduces to the ODE for $H(\eta)$,

[TABLE]

The localized nature of the self-similar dynamics implies slow evolution away from $x_{c}$, namely $h_{t}=O(1)$ for $|x-x_{c}|=O(1)$ as $\tau\to 0$ [65]. This yields a far-field boundary condition on the similarity solution,

[TABLE]

Using a one-parameter shooting method with $H(0)=H_{0}$ and $H^{\prime}(0)=0$ we numerically solve the symmetric similarity equation (61) subject to (62). For the unique value of $H(0)=H_{0}^{*}\approx 0.711$, we obtained a self-similar solution $H(\eta)$ of this second order boundary value problem, corresponding to the effective PDE

[TABLE]

This is in contrast with problems leading to a fourth-order problem for similarity solutions, where multiple solutions can co-exist [65].

A comparison between the PDE profiles scaled using $h_{\min}(t)=\min_{x}h(x,t)$ and the similarity solution in Fig. 17 (right) shows that the scaled PDE solution converges to $H(\eta)$ as $t\to t_{c}$. As rupture is approached, the energy (7) will be dominated by the local contribution near the rupture point, $\mathcal{E}\sim{1\over 3}\int h^{-3}\,dx=O(\tau^{-3/10})\to\infty$. The fact that the energy is increasing (and diverging) as $t\to t_{c}$ again shows that influence of the non-conserved term in breaking the gradient flow structure of (1) can be dominant in the dynamics. In [29], a fixed value for $\bar{\mathcal{P}}$ was used in the range $\bar{\mathcal{P}}>\mathcal{P}_{\max}$ so no steady states would exist and rupture would occur from all generic initial conditions. In this paper, we have $\bar{\mathcal{P}}$ in the critical range, and in the next section we will discuss how different dynamics and transitions occur depending on the initial data.

9 Numerical simulations of dynamic transitions

We have shown that the existence and stability of the coexisting steady states crucially depend on the domain size and the system parameters. The complex bifurcation structures also gives the potential for interesting dynamic transitions among these states. For instance, Fig. 8 (top) depicts a sequence of bifurcation points on the domain size $L$, ${\bar{\ell}_{+}}^{\,p}<{\bar{\ell}_{-}}^{\,\gamma}<\ell_{AB}^{\,(1)}<2{\bar{\ell}_{+}}^{\,p}<\cdots$, dividing the system into different regimes. With the domain size in each $L$-intervals, equation (5) has distinctive coexisting steady states and dynamic behaviors.

Here we focus on the parameters used in Fig. 8 and illustrate typical dynamics of the equation (5) for two choices of $L$, $L=40$ and $L=50$. In both cases, the principal second-order and fourth-order steady states $H_{A}(x)$ and $H_{B}(x)$ coexist with the two constant equilibria $\bar{H}_{-}$ and $\bar{H}_{+}$. However, we will show that the stability of these states change, leading to qualitatively different dynamics.

To study the effects of the dominant eigenmodes of a steady state $H(x)$, where $H(x)=H_{A}(x)$ or $H(x)=H_{B}(x)$, for $|\delta|\ll 1$ we perform PDE simulations of (5) starting from perturbed initial conditions (52), where $\Psi(x)$ is a normalized eigenmode associated with an unstable eigenvalue $\lambda$ and satisfies $\|\Psi\|_{2}=1$ and $\Psi(0)>0$. Consequently, the linearized dynamics takes the form $h(x,t)=H(x)+\delta\Psi(x)e^{\lambda t}$ when close to the steady state $H(x)$. The growth or decay of the perturbation in time can be quantified by the $L^{2}$-distance of the dynamic solution $h(x,t)$ and the corresponding steady state $H(x)$,

[TABLE]

For the first simulation, we use the domain size $L=40$ where the second-order steady state $H_{A}(x)$ is linearly stable up to translations, while the fourth-order steady state $H_{B}(x)$ is linearly unstable. As shown in Fig. 10, the largest nonzero eigenvalue of $H_{A}(x)$ is $\lambda_{A}\approx-0.0018$, and the dominant unstable eigenmode of $H_{B}(x)$, $\Psi_{B}(x)$, is associated with $\lambda_{B}\approx 0.017$. Fig. 18 (left) shows the PDE simulation of (5) starting from the initial data (52) with $\Psi=\Psi_{B}(x)$ and $\delta=-0.01$. In this simulation, the dynamic solution evolves away from $H_{B}(x)$ and approaches the stable steady state $H_{A}(x)$. Fig. 18 (right) gives the $L^{2}$-distances between the PDE solution and the two steady states. In the early stage, the distance $\|h(x,t)-H_{B}(x)\|_{2}$ grows exponentially in time following (63) at the growth rate $\lambda=\lambda_{B}$. As the solution approaches the steady state $H_{A}(x)$, the distance $\|h(x,t)-H_{A}(x)\|_{2}$ decays exponentially at the rate $\lambda=\lambda_{A}$. If we perturb the unstable steady state $H_{B}(x)$ in the opposite direction, and specify the initial data using (52) with $\Psi=\Psi_{B}(x)$ and $\delta=0.01$, then finite-time rupture phenomenon as shown in Fig. 17 will occur.

Figure 19 (left) gives a transition-state diagram that summarizes the dynamics among the coexisting states in this system with $L=40$. In addition to the transitions from $H_{B}(x)$ to $H_{A}(x)$ and from $H_{B}(x)$ to finite-time rupture induced by the most unstable eigenmodes $\Psi_{B}$, it also includes the dynamical transitions starting from spatially-uniform solutions. Spatially-uniform solutions $\bar{h}(t)$ of the PDE (5) are governed by the ODE

[TABLE]

which is obtained by dropping the spatial derivative terms from (5). Since $\gamma>0$, the evolution of such solutions is determined by the sign of $\Pi(\bar{h})$. Therefore, for $\bar{\mathcal{P}}$ in the critical range, spatially uniform perturbations around the uniform steady state $\bar{H}_{-}$ can lead the solution to approach either the uniform state $\bar{H}_{+}$ or the quenched state $h\equiv 0$. We also note that unbounded growth, $\bar{h}\to\infty$, is not possible in (64) in the critical range. PDE simulations starting from spatially-perturbed initial data near the $\bar{H}_{+}$ state lead to convergence to the stable steady state $H_{A}(x)$. Finite-time rupture is also numerically observed when the initial data is given by $\bar{H}_{-}$ with spatially nonuniform perturbations.

In contrast, for the domain size $L=50$ the four coexisting steady states $H_{A}(x)$, $H_{B}(x)$, $\bar{H}_{-}$ and $\bar{H}_{+}$ are all unstable. Fig. 19 (right) presents the corresponding transition-state diagram for $L=50$, showing that all of these steady states with perturbations can lead to finite-time rupture. Again we include arrows from $\bar{H}_{-}$ to $\bar{H}_{+}$ and $h\equiv 0$ that correspond to transitions induced by spatially uniform perturbations governed by (64).

Moreover, we include an arrow from $\bar{H}_{+}$ to $H_{B}(x)$ to represent transient dynamics between the two states. An example of such dynamics is depicted in Fig. 20 (top) where the initial perturbations around $\bar{H}_{+}$ lead the PDE solution to evolve towards the fourth-order state $H_{B}(x)$, whose instability then yields finite-time singularity at the edges of the domain. The $L^{2}$-distance between the PDE solution and the uniform steady state $\bar{H}_{+}$ in Fig. 20 (bottom left) shows that $\|h(x,t)-\bar{H}_{+}\|_{2}$ grows monotonically as the $h(x,t)$ evolves away from $\bar{H}_{+}$ with a growth rate $\lambda_{+}\approx 0.00412$. This value agrees with the prediction (16) from the stability analysis for the state $\bar{H}_{+}$. Meanwhile, the distance $\|h(x,t)-H_{B}(x)\|_{2}$ decreases for $0<t<2515$, indicating that the PDE solution evolves towards $H_{B}(x)$ in the early stage, but eventually develops a finite time singularity at a critical time $t_{c}\approx 3317.1$. The self-similar nature of the rupture solution as $t\to t_{c}$ is captured in Fig. 20 (bottom right) which follows $h_{\min}(t)=O\left((t_{c}-t)^{1/5}\right)$ from (60) with $\alpha=1/5$.

The eigenvalue plots in Fig. 10 also show that the system can have an unstable $H_{A}$ state and a coexisting stable $H_{B}$ state for a range of $L$ (say for $L=57$), which will yield a transition from $H_{A}$ to $H_{B}$ similar to the one shown in Fig. 18. For larger domain sizes, more coexisting steady states are involved in the system and other complicated dynamics like coarsening and symmetry breaking are also expected (see [30]); further work is needed on these issues.

10 Conclusions

We have studied interesting equilibrium structures and dynamics in the non-conserved thin film-type equation (5) with small positive $\gamma>0$. Unlike the conserved model (4) (equivalent to (5) with $\gamma=0$) and the volatile thin film model (equivalent to (5) with $\gamma<0$) where only second-order nonuniform steady states exist in the system, our equation also has fourth-order nonuniform steady state solutions. This is a consequence of not having a gradient flow structure.

In a critical range of the parameter $\bar{\mathcal{P}}$ we have investigated various local bifurcation structures that arise from the coexisting uniform and nonuniform steady states. For instance, the interactions between second-order and fourth-order branches yield different primary-primary and primary-secondary bifurcations. Further work is needed to more fully analyze the structure of the fourth-order steady states [54]. Stability analysis of the steady state solutions also reveals distinctive dynamic behaviors driven by the non-conserved effects, including the appearance of limit cycles, finite-time rupture, and dynamic transitions among steady states.

The weakly non-conserved limit $\gamma\to 0$ is of interest to be studied as a global bifurcation from the mass-conserving equation for $\gamma=0$. While we have shown that some properties persist in this limit, other features of the system change dramatically. For instance, second-order states exist for any value for $p$ in the case $\gamma=0$, but they are restricted to having only $p=0$ for $\gamma>0$. The appearance of limit cycles, the fourth-order steady states, and finite-time rupture for $\gamma>0$ also motivates further studies from the perspectives of PDE analysis and dynamical systems theory.

A number of interesting questions regarding the model (5) remain to be solved. First, while this paper focuses on the weak non-conserved effect limit with $\gamma\to 0$, we are also interested in the transient pattern formation under strong non-conserved influences (with large $\gamma$ values).

Moreover, we have used the Floquet theory to numerically identify the secondary bifurcations from the second-order states. A more careful asymptotic analysis with respect to spatial perturbations is needed to generalize these results to higher-order secondary bifurcations. For instance, we observe that the local structure for $\tilde{H}_{3}(L)\tilde{H}_{4}(L)$ for a base state $H_{A,k}$ takes the form of $\tilde{H}_{3}(L)\tilde{H}_{4}(L)=O([L-\ell_{AB}^{\,(k)}]^{k})$ (see Fig. 8 (bottom)).

The current analysis is limited to steady states and their stability on relatively small domains. For a larger-scale domain where more non-trivial steady states coexist, more complex bifurcation structures and rich pattern formation are expected and need careful investigation.