2 Two-peak interfacial dynamical system

Figure 1. Schematic of the system, consisting of two liquid layers sandwiched between a hot plate at temperature $T_{H}$ and a passive atmosphere at a far-field temperature $T_{\infty }$ . The dashed lines represent the interface positions in the unperturbed base state. $\unicode[STIX]{x1D6FE}_{1}(T)$ is the (temperature-dependent) interfacial tension of the inter-liquid interface, while $\unicode[STIX]{x1D6FE}_{2}(T)$ is the surface tension of the free interface of liquid 2. $\unicode[STIX]{x1D70C}_{i}$ , $\unicode[STIX]{x1D707}_{i}$ and $D_{i}$ denote densities, viscosities and thermal diffusivities, respectively.

Consider first two immiscible liquids, layered between a hot plate at temperature $T_{H}$ and a passive gas at a far-field temperature of $T_{\infty }$ , as shown in figure 1. Liquid 1, adjacent to the plate, has a mean thickness $H_{1}$ , while liquid 2 has a mean thickness $H_{2}-H_{1}$ . The interfacial tensions of both interfaces are temperature dependent: $\unicode[STIX]{x1D6FE}_{i}=\unicode[STIX]{x1D6FE}_{i\infty }-\unicode[STIX]{x1D6FD}_{i}(T-T_{\infty })$ with $(i=1,2)$ . The densities, viscosities and thermal diffusivities of the liquids are denoted by $\unicode[STIX]{x1D70C}_{i}$ , $\unicode[STIX]{x1D707}_{i}$ and $D_{i}$ . Then focus on the case of thin viscous films, for which the unsteady and convective contributions to the momentum and energy balances can be neglected. Under these conditions, the two interfaces are unstable solely to long-wavelength disturbances, and $\unicode[STIX]{x1D706}$ , the length scale of longitudinal variations along $x$ , is much larger than the local film thickness, which is of order $H_{2}$ . Taking advantage of this separation of length scales, we reduce the Navier–Stokes and thermal energy balance equations to a pair of long-wave evolution equations for the interface positions ( $h_{1}$ and $h_{2}$ , cf. figure 1). These equations, non-dimensionalized with length scale $l_{c}=H_{2}$ , velocity scale $u_{c}=\unicode[STIX]{x1D6FE}_{2\infty }/\unicode[STIX]{x1D707}_{2}$ and time scale $t_{c}=\unicode[STIX]{x1D706}/u_{c}$ , are presented below:

Here, $Q_{1}$ and $Q_{2}$ denote the respective phase flow rate. $Ma=(T_{H}-T_{\infty })\unicode[STIX]{x1D6FD}_{2}/\unicode[STIX]{x1D6FE}_{2\infty }$ is the Marangoni number, and $Bo=H_{2}^{2}\unicode[STIX]{x1D70C}_{2}g/\unicode[STIX]{x1D6FE}_{2\infty }$ is the Bond number. $Bi=UH_{2}/k_{2}$ is the Biot number, wherein $U$ is the heat transfer coefficient between liquid 2 and the passive gas and $k_{2}$ is the thermal conductivity of liquid 2. $\unicode[STIX]{x1D707}_{r}$ , $\unicode[STIX]{x1D70C}_{r}$ , $\unicode[STIX]{x1D6FD}_{r}$ and $k_{r}$ are ratios of the respective properties of liquid 1 to liquid 2, and $\unicode[STIX]{x1D6FE}_{r}=\unicode[STIX]{x1D6FE}_{1\infty }/\unicode[STIX]{x1D6FE}_{2\infty }$ .

These reduced equations are derived using the weighted residual integral boundary layer method (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000; Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2013), assuming a small Reynolds number, $Re\sim O(H_{2}/\unicode[STIX]{x1D706})$ , and a moderate Prandtl number, $Pr\sim O(1)$ (or the Péclet number $Pe=RePr\sim O(H_{2}/\unicode[STIX]{x1D706})$ ). For simplicity, we have dropped all terms of order $(H_{2}/\unicode[STIX]{x1D706})^{2}$ and smaller, i.e. momentum and thermal inertia, longitudinal thermal and viscous diffusion, and normal viscous surface stresses are neglected. The transverse temperature profile is a linear conduction profile that is slaved to the local interface positions. Equivalent equations can be derived using the long-wave expansion method (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Nepomnyashchy, Simanovskii & Legros Reference Nepomnyashchy, Simanovskii and Legros2012).

The two fluid system described by these equations is unstable to the Marangoni instability if $Ma>0$ ( $T_{H}>T_{\infty }$ , $\unicode[STIX]{x1D6FD}>0$ ). The Rayleigh–Taylor instability destabilizes the free interface if $Bo>0$ (gravity acts in the positive $z$ direction), while the inter-liquid interface is destabilized if $Bo>0$ and $\unicode[STIX]{x1D70C}_{r}>1$ , or $Bo<0$ and $\unicode[STIX]{x1D70C}_{r}<1$ . For the case of equal mean thicknesses of the liquids, the unstable quiescent base state corresponds to $(h_{1},h_{2})=(0.5,1)$ . Equations (2.1) are linearized about this state, considering infinitesimal perturbations of the form $({\hat{h}}_{1},{\hat{h}}_{2})\exp (\text{i}\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D70E}t)$ . On solving the resulting eigenvalue problem for $\unicode[STIX]{x1D70E}$ , one obtains the linear growth rate ( $\unicode[STIX]{x1D70E}_{r}=\text{Re}[\unicode[STIX]{x1D70E}]$ ) as a function of the wavenumber $\unicode[STIX]{x1D6FC}$ (the corresponding plot is called the dispersion curve).

Figure 2. Two-peak linear growth curves (dispersion curves), obtained from a normal-mode linear stability analysis of (2.1), considering a quiescent base state with flat liquid layers of equal depth. The solid (blue) and dashed (red) lines correspond to the two eigenmodes, with in-phase and out-of-phase interface deflection respectively. Results are shown for two sets of parameter values, which constitute the two principal cases studied in this work. (a) TP1: $H_{1}/H_{2}=0.5$ , $Bo=0.1$ , $Ma=0.5$ , $\unicode[STIX]{x1D6FE}_{r}=1.5$ , $\unicode[STIX]{x1D6FD}_{r}=0.85$ , $\unicode[STIX]{x1D707}_{r}=3$ , $\unicode[STIX]{x1D70C}_{r}=27.697$ , $Bi=1.4$ , $k_{r}=1.25$ ; (b) TP2: $H_{1}/H_{2}=0.5$ , $Bo=0$ , $Ma=0.2$ , $\unicode[STIX]{x1D6FE}_{r}=0.243$ , $\unicode[STIX]{x1D6FD}_{r}=2$ , $\unicode[STIX]{x1D707}_{r}=1.2$ , $\unicode[STIX]{x1D70C}_{r}=2$ , $Bi=1$ , $k_{r}=2$ .

A special feature of this system is that, for appropriate combinations of fluid properties, the dispersion curve can have two peaks with comparable growth rates. Dispersion curves for two such cases, henceforth referred to as TP1 and TP2, are presented in panels 2(a,b). Here, for each wavenumber $\unicode[STIX]{x1D6FC}$ , there are two values of $\unicode[STIX]{x1D70E}$ corresponding to two distinct modes. One mode (solid blue line) corresponds to in-phase interface deflections, wherein a crest of the top interface coincides with a crest of the bottom interface ( ${\hat{h}}_{1}$ and ${\hat{h}}_{2}$ have the same sign). The other mode (red dashed line) corresponds to out-of-phase interface deflections, wherein a crest of the top interface coincides with a trough of the bottom interface ( ${\hat{h}}_{1}$ and ${\hat{h}}_{2}$ have the opposite sign). In both cases, TP1 and TP2, there exists a range of wavenumbers (where the solid and dashed curves meet) for which these two modes are oscillatory and their eigenvalues form a complex-conjugate pair. However, the peaks of the dispersion curves are located in regions where the modes are non-oscillatory ( $\text{Im}[\unicode[STIX]{x1D70E}]=0$ ). Here, the in-phase deflection mode is dominant, as its growth rate is significantly larger than that of the out-of-phase deflection mode. Therefore, in what follows we focus solely on the in-phase deflection mode (solid blue line).

Numerical solutions of (2.1) are obtained using a second-order finite difference scheme for spatial discretization. Time integration is carried out with the NDSolve routine of Mathematica^®, which implements adaptive time stepping. The simulations are terminated once the thickness of either fluid reaches 0.01 times its initial value. By this stage, interface deformation ought to have progressed far enough for the nonlinear pattern to fully emerge. Moreover, the interfaces would be close to either pinching off or rupturing at the wall, and following the evolution further would require the introduction of intermolecular forces, such as van der Waals forces (Israelachvili Reference Israelachvili2011). These forces are important only for ultra-thin films and do not affect the pattern that is selected as the peak-modes of figure 2 transit from linear to nonlinear growth.

3 Nonlinear reinforcement of the fastest-growing linear mode

In this section, we use a simple example to show how the fastest-growing linear mode dominates interfacial patterns in the case of a single-peak dispersion curve. For this, we consider the Marangoni instability of a single fluid film, in the absence of gravity, obtained by setting $h_{1}=0$ , $\unicode[STIX]{x1D6FE}_{r}=0$ , $\unicode[STIX]{x1D6FD}_{r}=0$ , $\unicode[STIX]{x1D707}_{r}=1$ , $\unicode[STIX]{x1D70C}_{r}=1$ , $k_{r}=1$ and $Bo=0$ in (2.1). The corresponding dispersion curve has a single peak, as depicted in panel 3(a). The originally flat and static interface is subjected to two separate perturbations: one corresponding to the peak-mode ( $\exp (\text{i}\unicode[STIX]{x1D6FC}_{m}x)$ ) and another corresponding to a slower mode of twice the wavelength ( $\exp (\text{i}\unicode[STIX]{x1D6FC}_{m}x/2)$ ). Both calculations are carried out on a periodic domain of length $4\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{m}$ . Although the initial form of the perturbed interface is quite different in the two cases (compare the two panels in figure 3 b), the pattern that emerges at later times is indistinguishable (compare figure 3 c,d).

Figure 3. Dominance of the fastest-growing linear mode when the dispersion curve has a single peak. (a) Linear dispersion curve, corresponding to a single fluid on a hot plate ( $Ma=2$ , $Bi=1$ , $Bo=0$ ), with a peak at $\unicode[STIX]{x1D6FC}_{m}=0.613$ . (b) Spatial form of linear modes corresponding to wavenumbers $\unicode[STIX]{x1D6FC}_{m}$ (top panel) and $\unicode[STIX]{x1D6FC}_{m}/2$ (bottom panel). The corresponding nonlinear evolution is shown in (c,d). The non-dimensional time is mentioned above each panel.

How is it that the peak-mode is able to dominate the pattern even when it is absent from the initial perturbation? This occurs due to nonlinearities that generate the peak-mode $\unicode[STIX]{x1D6FC}_{m}$ as a higher harmonic of the initial longer-wave disturbance $\unicode[STIX]{x1D6FC}_{m}/2$ . Once generated, the $\unicode[STIX]{x1D6FC}_{m}$ mode quickly out paces the $\unicode[STIX]{x1D6FC}_{m}/2$ mode and dominates the emerging pattern. In this problem, the nonlinearity originates from (i) the Marangoni surface stress that depends on interfacial temperature gradients, which in turn vary in response to the evolving interface ( $Ma\unicode[STIX]{x2202}_{x}[T_{i}(h_{i})]$ in (2.1)) and (ii) viscous stresses that retard fluid flow, to an extent that increases as the film approaches the no-slip wall (terms of the form $h_{i}^{2}/2$ and $h_{i}^{3}/3$ in (2.1)). On imposing the initial perturbation of wavenumber $\unicode[STIX]{x1D6FC}_{m}/2$ , Marangoni forces begin to drain fluid symmetrically from under the single trough (cf. bottom panel of figure 3 b). However, viscous stresses, which retard this flow, are strongest under the trough. Consequently, the largest outflow occurs not directly under the trough, but at two points located symmetrically on either side. The film height decreases faster at these locations, and the original trough eventually transforms into a crest, sandwiched between two new secondary troughs (cf. top panel of figure 3 d). This process is aided by the fact that the Marangoni flow near the original trough changes direction along with the slope of the interface, so that fluid drains out of the secondary troughs into the newly formed crest. In this manner, shorter-scale structures are generated spontaneously, a short time after imposing the initial perturbation. These include the $\unicode[STIX]{x1D6FC}_{m}$ mode, which then grows much faster and quickly dominates the evolution of the interface.

The physical origin of nonlinearities capable of generating shorter-wavelength modes is different in each problem. In the capillary break-up of jets, where there are no wall effects, capillary pressure due to transverse curvature, inertia and the presence of an outer viscous fluid are all sources of structure-generating nonlinearity (Eggers Reference Eggers1997; Pozrikidis Reference Pozrikidis1999). Whatever their origin, these nonlinearities act to reinforce the fastest linear mode so that it ultimately dominates the pattern, regardless of the amplitude of the peak-mode in the initial disturbance. This explains why it is possible to anticipate the dominant length scale in nonlinear patterns from the peak of a linear dispersion curve and why the purely linear theory of interfacial pattern formation put forth by Rayleigh (Reference Rayleigh1879a ) found reasonable agreement with experiment.

This section has shown how nonlinear effects enhance the growth of the mode with the fastest linear growth rate. In the next section, we eliminate this linear growth bias by considering situations wherein the dispersion curve has two peaks with equal growth rates. The basic question we ask is: will both peak-modes be present in the final pattern or will nonlinearity select one of them?

4 Nonlinear selection in the context of two-peak dispersion curves

This section presents results of numerical simulations of (2.1), for cases that have dispersion curves with two peaks of comparable height. To focus attention on the peak-modes, the simulations are carried out in a finite periodic domain of length equal to the wavelength of the longer-wave peak-mode, given by $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{L}$ where $\unicode[STIX]{x1D6FC}_{L}$ is the wavenumber of the longer-wave peak. In the simulations, we perturb the $\unicode[STIX]{x1D6FC}_{L}$ peak, as well as a commensurate mode, $\unicode[STIX]{x1D6FC}_{S}=n\unicode[STIX]{x1D6FC}_{L}$ , which is closest to the shorter-wavelength peak. The peak-modes in figure 2 are very close to being commensurate, with $n=2$ for case TP1 and $n=4$ for case TP2. (This feature was a major reason for selecting these particular cases for study.) We consider small initial perturbations of the form:

where the eigenvectors $({\hat{h}}_{1L},{\hat{h}}_{2L})$ and $({\hat{h}}_{1S},{\hat{h}}_{2S})$ are of unit norm. We compare the interfacial patterns that emerge after exciting (i) only the longer-wave peak ( $A_{0}=0.001,B_{0}=0$ ), (ii) only the shorter-wave peak ( $A_{0}=0,B_{0}=0.001$ ) and (iii) both peaks in equal proportion ( $A_{0}=B_{0}=0.001$ ). The two peak-mode disturbances are allowed to have a phase difference $\unicode[STIX]{x1D719}$ . In what follows, we present results for cases TP1 and TP2 (panels (a,b) of figure 2), which serve to illustrate all the key findings of this work. The effect of extending the domain length and admitting modes longer than $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{L}$ is examined later, in § 6.

Figure 4 presents results for case TP1 (cf. panel 2 a), wherein the two peak-modes have equal linear growth rates. Column (a) shows the evolution of the interfacial pattern that results from perturbing the longer-wave peak-mode, $\unicode[STIX]{x1D6FC}_{L}=0.477$ . Column (b) shows the result of perturbing the shorter-wave peak-mode, $\unicode[STIX]{x1D6FC}_{S}=2\unicode[STIX]{x1D6FC}_{L}=0.954$ . The corresponding eigenvectors are $({\hat{h}}_{1L},{\hat{h}}_{2L})=(0.253,0.968)$ and $({\hat{h}}_{1S},{\hat{h}}_{2S})=(0.820,0.573)$ . Apart from the difference in wavelengths, the qualitative nature of the two patterns is very different. The top interface is more active in the case of the $\unicode[STIX]{x1D6FC}_{L}$ perturbation, and the two interfaces meet each other and pinch off, before the bottom interface meets the wall (figure 4 a). In contrast, for the $\unicode[STIX]{x1D6FC}_{S}$ disturbance, the bottom interface deforms more than the top interface, leading to rupture of the bottom interface at the wall (figure 4 b). When both the peak-modes are perturbed simultaneously, these two contrasting patterns will grow together and possibly compete. The principle of the fastest-growing linear mode suggests that both modes should be present in the final pattern, because the two peak-modes have identical linear growth rates. However, the result shown in column (c) of figure 4 is in complete contradiction with this line of reasoning. The shorter-wavelength $\unicode[STIX]{x1D6FC}_{S}$ mode is seen to completely dominate the interfacial pattern – all evidence of the longer-wavelength peak-mode is completely eliminated, even at relatively early times.

Figure 4. Nonlinear selection among the two peak-modes, corresponding to case TP1 (cf. figure 2 a). Columns (a)–(d) depict the evolution of different initial perturbations: (a) only the longer-wavelength peak-mode ( $\unicode[STIX]{x1D6FC}_{L}=0.477$ ), (b) only the shorter-wavelength peak-mode ( $\unicode[STIX]{x1D6FC}_{S}=2\unicode[STIX]{x1D6FC}_{L}=0.954$ ), (c) equal in-phase superposition of peak-modes ( $A_{0}=0.001,B_{0}=0.001,\unicode[STIX]{x1D719}=0$ in (4.1)), (d) equal, out-of-phase superposition of peak-modes ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ). The eigenvectors corresponding to the two peak-modes are $({\hat{h}}_{1L},{\hat{h}}_{2L})=(0.253,0.968)$ and $({\hat{h}}_{1S},{\hat{h}}_{2S})=(0.820,0.573)$ . The non-dimensional time corresponding to each interfacial pattern is mentioned above the respective panel.

Figure 5. The influence of phase difference and a counter-example to Rayleigh’s principle in the context of case TP2 (cf. figure 2 b). Each panel shows the final pattern obtained just before terminating the simulation, along with the corresponding non-dimensional time (mentioned above each panel). In (a–d), the two excited peak-modes, $\unicode[STIX]{x1D6FC}_{L}=0.196$ and $\unicode[STIX]{x1D6FC}_{S}=4\unicode[STIX]{x1D6FC}_{L}=0.784$ , have equal linear growth rates of 0.00043. The eigenvectors corresponding to the two peak-modes are $({\hat{h}}_{1L},{\hat{h}}_{2L})=(0.165,0.986)$ and $({\hat{h}}_{1S},{\hat{h}}_{2S})=(0.993,0.120)$ . The perturbations in each panel are: (a) only $\unicode[STIX]{x1D6FC}_{L}$ , (b) only $\unicode[STIX]{x1D6FC}_{S}$ , (c) in-phase superposition ( $A_{0}=0.001,B_{0}=0.001,\unicode[STIX]{x1D719}=0$ in (4.1)), (d) out-of-phase superposition ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ). Panels (e–h) show how the final pattern for equal in-phase perturbations ( $A_{0}=0.001,B_{0}=0.001,\unicode[STIX]{x1D719}=0$ ) changes as the growth rate of the shorter-wave peak-mode ( $\unicode[STIX]{x1D6FC}_{S}$ ) is reduced in four stages (by increasing the interfacial tension ratio $\unicode[STIX]{x1D6FE}_{r}$ ): (e) 0.00042, (f) 0.00040, (g) 0.00038, (h) 0.00033.

The situation becomes even more intriguing when the effect of a phase difference between the two peak-mode perturbations is considered. Figure 4(d) shows that the pattern changes completely when the perturbations are out of phase ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ in (4.1)), compared to the case when they are in phase ( $\unicode[STIX]{x1D719}=0$ , cf. figure 4 c). The nature of pinch off in figure 4(d) is now governed by the $\unicode[STIX]{x1D6FC}_{L}$ mode, although contributions from both peak-modes are present (with $\unicode[STIX]{x1D6FC}_{L}$ primarily influencing the top interface and $\unicode[STIX]{x1D6FC}_{S}$ impacting the bottom interface). Linear stability calculations are unaffected by phase differences, and consequently, the linear dispersion curve is powerless to explain such phase-dependent pattern evolution.

In hindsight, the influence of a phase difference can be understood by comparing the early evolution of the two independent perturbations to the peak-modes. The $\unicode[STIX]{x1D6FC}_{L}$ mode causes the interfaces to move downward at the centre of the domain (figure 4 a), whereas the $\unicode[STIX]{x1D6FC}_{S}$ mode causes the interfaces to move upward at the same location (figure 4 b). Therefore, the two modes compete when they are excited in phase, with the result that $\unicode[STIX]{x1D6FC}_{S}$ dominates (figure 4 c). On the other hand, out-of-phase disturbances can grow together and coexist in the final pattern, as seen in figure 4(d). Thus, the influence of phase differences is tied to the interaction between crests and troughs of the two peak-mode patterns. By this logic, phase differences should cease to matter when the two peak-modes have very different wavelengths. In this case, a trough of the longer-wave pattern will overlap with several crests and troughs of the shorter-wave pattern, and the net interaction of the two patterns should not be significantly influenced by the phase difference. This holds true for case TP2 (wherein $\unicode[STIX]{x1D6FC}_{S}=4\unicode[STIX]{x1D6FC}_{L}$ , cf. panel 2 b), as demonstrated by panels 5(a–d). Panels 5(a,b) show the final pattern resulting from independent perturbations of the $\unicode[STIX]{x1D6FC}_{L}$ and $\unicode[STIX]{x1D6FC}_{S}$ modes respectively. Panels 5(c,d) show the results of in-phase and out-of-phase superposed peak-mode perturbations. As expected, the phase difference does not significantly affect the qualitative nature of the pattern, which is completely dominated by the shorter-wave $\unicode[STIX]{x1D6FC}_{S}$ peak-mode.

The erstwhile dominant peak-mode (shorter-wavelength peak) continues to dominate even when its growth rate is slightly reduced. This result is demonstrated by panel 5(e), which corresponds to a slightly modified case of TP2 (cf. panel 2 b), in which the $\unicode[STIX]{x1D6FC}_{S}$ mode has a slightly smaller growth rate than the $\unicode[STIX]{x1D6FC}_{L}$ mode (cf. the caption of figure 5). This is a definite counter-example to Rayleigh’s paradigm of the fastest-growing linear mode, because in this case the mode with a smaller linear growth rate dominates the pattern! Of course, as the growth rate of the $\unicode[STIX]{x1D6FC}_{S}$ mode is decreased further (panels 5 f–h), the system eventually behaves like one with a single-peak dispersion curve and the $\unicode[STIX]{x1D6FC}_{L}$ mode sets the pattern (compare panels 5 a,h). A similar result is found for in-phase perturbations of the peak-modes of case TP1. The shorter $\unicode[STIX]{x1D6FC}_{S}$ mode (cf. figure 4 b) continues to dominate even when its growth rate is slightly reduced. Further reduction, however, eventually leads to the longer $\unicode[STIX]{x1D6FC}_{L}$ peak-mode (cf. figure 4 a) dominating.

5 Low-order ODE model: a truncated Galerkin approach

The previous section has demonstrated that interfacial pattern evolution cannot be anticipated from linear stability theory alone, in cases where the dispersion curve has two peaks of comparable height. Numerical investigations tell us that either peak-mode can dominate the pattern and that the outcome may also depend on the phase difference between initial perturbations, but not always. Anticipating the final pattern requires one to account for nonlinear interactions between the peak-modes. At the same time, full numerical simulations are computationally expensive and may be unfeasible in many situations, especially when different types of perturbations must be considered. An alternative is to develop a low-order ordinary differential equation model that accounts for only the key interactions between peak-modes, in order to predict which of the two dominates. In systems close to the onset of instability, at codimension-two supercritical points, such amplitude equations have been derived using asymptotic techniques like multiple scale expansions (Cross & Hohenberg Reference Cross and Hohenberg1993) or the centre manifold reduction (Carini, Auteri & Giannetti Reference Carini, Auteri and Giannetti2015; Roberts Reference Roberts2015). In this work, we are primarily concerned with systems that are far above the instability threshold and undergoing subcritical bifurcations. Moreover, we wish to describe the evolution of the most unstable peak-modes and not those close to the cutoff neutral wavenumber. Therefore, rather than applying a strict asymptotic procedure, we derive a low-order model by carrying out a Galerkin projection of the governing equations (2.1) onto a set of judiciously selected modes (Balmforth, Provenzale & Whitehead Reference Balmforth, Provenzale, Whitehead, Balmforth and Provenzale2001).

By comparing various truncated models with full simulations of the partial differential equation (PDE) model (2.1), we have found that the lowest-order ODE model capable of adequately describing the system is one based on three modes – the two peak-modes and a third mode that facilitates their interaction. The interface heights are represented as a combination of these modes, as follows:

Here $\unicode[STIX]{x1D6FC}_{S}=n\unicode[STIX]{x1D6FC}_{L}$ , where $n$ is an integer, as we are considering commensurate peak-modes. The first two modes, with amplitudes $A(t)$ and $B(t)$ , correspond to the longer-wave and shorter-wave peak-modes respectively. The third mode, with amplitude $C(t)$ , has the eigenvector of the longer-wave peak-mode ( ${\hat{h}}_{1L},{\hat{h}}_{2L}$ ) but the spatial variation of the shorter-wave peak-mode ( $\text{e}^{(\text{i}\unicode[STIX]{x1D6FC}_{S})x}$ ). This mode is especially important when the two eigenvectors ( ${\hat{h}}_{1L},{\hat{h}}_{2L}$ ) and ( ${\hat{h}}_{1S},{\hat{h}}_{2S}$ ) are orthogonal. In this case, a nonlinear term (quadratic if $n=2$ , as in case TP1) acting on the $\unicode[STIX]{x1D6FC}_{L}$ mode will not make any contribution to the equation for the $\unicode[STIX]{x1D6FC}_{S}$ mode, after Galerkin projection. The third mode, on the other hand, will always be affected nonlinearly by the $\unicode[STIX]{x1D6FC}_{L}$ mode while simultaneously contributing to the equation for the $\unicode[STIX]{x1D6FC}_{S}$ mode, thereby facilitating communication between the two peak-modes. The importance of considering this third interaction mode is illustrated in figure 8, to be introduced later.

Expression (5.1) is substituted into (2.1), which is then projected onto the three modes. Projection involves taking the scalar product of the equations with the mode, followed by spatial integration. After equating real and imaginary parts, this results in six nonlinear coupled ODEs for the amplitudes $A(t)$ , $B(t)$ and $C(t)$ and the phases $\unicode[STIX]{x1D719}_{1}(t)$ , $\unicode[STIX]{x1D719}_{2}(t)$ and $\unicode[STIX]{x1D719}_{3}(t)$ . For simplicity, only the lowest-order nonlinear terms that enable the peak-modes to interact are retained. For case TP1, second-order terms are required, because the peak-mode wavenumbers are in a 1 : 2 ratio, whereas, terms up to fourth order must be retained for case TP2 where the peak-modes are in a 1 : 4 ratio (cf. panels 2 a,b). (It was verified that higher-order terms do not modify the prediction of the dominant mode.) The following initial conditions were used, to simulate the evolution of an initial perturbation to both peak-modes: $A=B=0.001$ , $C=0$ , $\unicode[STIX]{x1D719}_{1}=0$ , $\unicode[STIX]{x1D719}_{2}=\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D719}_{3}=0$ . Here, $\unicode[STIX]{x1D719}$ is the phase difference between the initial peak-mode perturbations, as defined in (4.1), and the disturbance amplitude of 0.001 is the same as that used in the simulations of the PDE model (2.1).

Figure 6. Low-dimensional description of the evolution of the peak-modes for case TP1 of panel 2(a). $A(t)$ is the amplitude of the longer-wave peak-mode while $B(t)$ is the amplitude of the shorter-wave peak-mode. The solid lines represent the amplitudes of the peak-modes (cf. (5.1)) extracted from the simulation of the PDE model (2.1), while the dashed lines are the predictions of the three-mode ODE model. Panels (a,b) depict the phase plane trajectory and temporal evolution of the amplitudes, for the case of in-phase peak-mode perturbations ( $\unicode[STIX]{x1D719}=0$ in (4.1)). Panels (c,d) correspond to out-of-phase peak-mode perturbations ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ in (4.1)). The associated spatial evolution is presented in panels 4(c,d), respectively. In (a,c), time progresses along the phase plane trajectories, starting at the origin. The dotted lines in these panels are diagonals.

Figure 7. Low-dimensional description of the evolution of the peak-modes for case TP2 of panel 2(b). $A(t)$ is the amplitude of the longer-wave peak-mode while $B(t)$ is the amplitude of the shorter-wave peak-mode. The solid lines represent the amplitudes of the peak-modes (cf. (5.1)) extracted from the simulation of the PDE model (2.1), while the dashed lines are the predictions of the three-mode ODE model. Panels (a,b) depict the phase plane trajectory and temporal evolution of the amplitudes, for the case of in-phase peak-mode perturbations ( $\unicode[STIX]{x1D719}=0$ in (4.1)). Panels (c,d) correspond to out-of-phase peak-mode perturbations ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ in (4.1)). The final spatial patterns corresponding to these two cases are presented in panels 5(c,d), respectively. The dotted lines in (a,c) are diagonals.

Figure 6 compares the prediction of this three-mode ODE model (dashed lines) with the amplitudes extracted from simulations of the PDE model (solid lines), for case TP1 (cf. panel 2 a). Panels 6(a,c) depict phase plane plots for in-phase ( $\unicode[STIX]{x1D719}=0$ ) and out-of-phase ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ) initial peak-mode perturbations, while panels 6(b,d) show the corresponding temporal evolution. The ODE model is able to capture the qualitative features of the evolution of the peak-modes, including the influence of the phase difference. It successfully predicts the dominance of the shorter-wave peak-mode for in-phase perturbations (panel 6 a), and the co-existence of both peak-modes for out-of-phase perturbations (with the longer-wave peak-mode playing a larger role, cf. panel 6 c). Quantitative accuracy of the model is, however, restricted to small values of the amplitudes (panels 6 b,d). As the amplitudes of the peak-modes increase, the neglected high-order interaction terms (cubic and greater in this case) become increasingly relevant. Moreover, higher harmonics of the peak-modes get excited and may begin to influence the growth of the peak-modes. These effects cannot be captured by the truncated three-mode ODE model. The fact that the ODE model is able to predict the dominant mode, in spite of these limitations, implies that pattern selection is largely determined by direct nonlinear interaction between the peak-modes, at a relatively early stage of their evolution.

Figure 7 presents results for case TP2 (cf. panel 2 b). The ODE model successfully predicts the dominant mode to be the shorter-wave peak-mode, for both in-phase (panels 7 a,b) and out-of-phase perturbations (panels 7 c,d). If the interfacial tension ratio is increased so as to slightly decrease the linear growth rate of the shorter-wave peak-mode, the ODE model continues to predict the dominance of the shorter-wave mode, in agreement with panels 5(c,e).

Similar to the results for case TP1 (cf. figure 6), the ODE model is quantitatively accurate only for small amplitudes. Beyond $t\approx 11\,000$ , the ODE prediction for $B(t)$ increases much more rapidly than the PDE results (cf. panels 7 b,d). In fact the ODE solution displays a singularity at $t\approx 11\,800$ . This non-physical behaviour may be cured by including additional modes in the ODE model. However, this is not necessary if the primary goal is to predict which peak-mode is dominant in the final pattern.

Figure 8. Demonstration of the importance of the interaction mode (of amplitude $C(t)$ in (5.1)) in the reduced-order ODE model. Phase plane orbits (a) and the temporal evolution (b) of the peak-modes are presented, for the system in panel 2(b), considering in-phase perturbations. (a) The solid lines are obtained from a simulation of the PDE model (2.1) and are identical to the solid lines in panels 7(a,b). The dashed lines are predictions of a two-mode ODE model that does not include the interaction mode. These predictions should be compared with the dashed lines in 7(a,b), obtained from the three-mode ODE model which retains the interaction mode. While the three-mode ODE model is qualitatively accurate and predicts the dominance of the shorter-wave peak-mode, the two-mode ODE model fails entirely and predicts the opposite result.

A particularly encouraging aspect of the ODE model is its ability to capture the influence of the phase difference between perturbations, which significantly impacts the final pattern in case TP1 (compare panel 4 c with 4 d and panel 6 a with 6 c). As the model contains only three modes, it may be possible to gain insight into inter-peak interaction and the role of phase differences, by a direct examination of the ODEs. To obtain the equations for case TP1, the parameter values corresponding to case TP1 (cf. caption of panel 2 a) are substituted into the full PDE model (2.1), which is then projected onto the three modes (5.1). The calculations are carried out with the aid of symbolic computing in Mathematica^®. The equations take on a relatively simple form if we restrict ourselves to in-phase and out-of-phase perturbations ( $\unicode[STIX]{x1D719}=0,\unicode[STIX]{x03C0}$ in (4.1)):

Here, $\unicode[STIX]{x1D70E}$ and $\unicode[STIX]{x1D714}$ are the individual linear growth rates of each mode and $m_{i}$ are positive constants. The values of these coefficients are specific to the parameter values of case TP1: $\unicode[STIX]{x1D70E}=0.008$ , $\unicode[STIX]{x1D714}=-0.063$ , $m_{1}=0.137$ , $m_{2}=0.053$ , $m_{3}=0.0003$ , $m_{4}=0.002$ , $m_{5}=0.012$ . The coefficient $\unicode[STIX]{x1D70E}$ is the same for the two peak-mode amplitudes $A$ and $B$ because we have used parameter values corresponding to case TP1, for which the peak-modes have equal linear growth rates.

The other terms in (5.2) account for inter-mode interaction, and include nonlinearities up to second order (the lowest order required for nonlinear interaction between the peak-modes, which are in a 1 : 2 wavelength ratio). Particularly important are the third terms on the right-hand side of (5.2a ) and (5.2b ). For in-phase perturbations ( $\unicode[STIX]{x1D719}=0$ ), these terms show that the two peak-modes compete, with the growth of $B$ suppressing that of $A$ (via $-m_{2}AB$ ) and vice versa (via $-m_{4}A^{2}$ ). Because $m_{2}>m_{4}$ , in this case, the shorter-wave mode dominates the pattern. For out-of-phase perturbations ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ), on the other hand, the sign of these terms reverse, indicating that the two modes cooperate and allow each other to coexist in the final pattern. This competition/cooperation stems from the manner in which crests and troughs of the two peak-modes overlap, as was discussed pictorially in § 4.

We have tested the three-mode ODE model for other cases as well, and have found similar positive results, with regard to the qualitative nature of the peak-mode interaction and the prediction of the dominant mode. For the sake of brevity, we only present one additional example in appendix A, for which the peak-mode wavelengths are in a 1 : 3 ratio. In this particular case, both modes play nearly equal roles in establishing the final pattern and neither is able to dominate.

We conclude this section by demonstrating the importance of including the interaction mode, of amplitude $C(t)$ , in the ODE model. This mode is particularly important for case TP2 (cf. panel 2 b), because the eigenvectors corresponding to the peak-modes are close to being orthogonal $({\hat{h}}_{1L},{\hat{h}}_{2L})=(0.165,0.986)$ and $({\hat{h}}_{1S},{\hat{h}}_{2S})=(0.993,0.120)$ . In this case, although fourth-order nonlinearities would generate a harmonic of the longer-wave mode (of wavenumber $\unicode[STIX]{x1D6FC}_{L}$ ) with a spatial variation matching the shorter-wave mode (of wavenumber $4\unicode[STIX]{x1D6FC}_{L}$ ), its projection onto the latter would be very small. Thus a two-mode ODE model, with just the peak-modes, will not be able to capture their interaction. This is shown by the results in figure 8, which corresponds to in-phase perturbations of the peak-modes of case TP2. The two-mode model completely fails, as it predicts the dominance of the longer-wave peak-mode, when in fact the shorter-wave mode dominates.

This deficiency of the two-mode model has been overcome by including just one extra mode – the interaction mode, which combines the spatial variation of the shorter-wave peak-mode ( $\unicode[STIX]{x1D6FC}_{S}$ ) with the eigenvector of the longer-wave peak-mode ( ${\hat{h}}_{1L},{\hat{h}}_{2L}$ ) (cf. (5.1)). The longer-wave mode can strongly excite the interaction mode, which in turn influences the evolution of the shorter-wave peak-mode via linear coupling. For case TP1, this interaction is represented by the second terms on the right-hand side of (5.2c ) and (5.2b ). These terms do not play a major role in case TP1 because the eigenvectors are far from being orthogonal, $({\hat{h}}_{1L},{\hat{h}}_{2L})=(0.253,0.968)$ and $({\hat{h}}_{1S},{\hat{h}}_{2S})=(0.820,0.573)$ , which enables a strong direct interaction between the peak-modes. In contrast, the interaction mode is important in case TP2 where direct interaction is weak, and it is only by including the interaction mode that the ODE model is able to successfully identify the dominant peak-mode (compare panel 7 a with 8 a).

6 Coexisting patterns on long domains

In the numerical simulations of § 4, the non-dimensional domain length ( $L$ ) was chosen to be just large enough to accommodate the longer-wave peak-mode, i.e. $L=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{L}$ . Increasing the domain length further admits modes of larger wavelengths. Due to their smaller growth rates, one may expect them to be dominated by the peak-modes. However, these long-wavelength modes may be able to affect the interfacial pattern by modulating the local phase difference between the peak-modes. In situations where the pattern depends on the phase difference, such as case TP1 (cf. panels 4 c,d), this phase modulation could result in the emergence of coexisting patterns on long domains.

Figure 9. Emergence of coexisting patterns on a long domain, $L=4\times 2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{L}$ . (a) The spatially varying phase difference $\unicode[STIX]{x1D719}(x)$ used in the initial perturbation (4.1) to the peak-modes. (b) The coexisting pattern that results from perturbing the peak-modes of case TP1 (panel 2 a) with the spatially varying phase difference given in (a). (c) The result of perturbing all admissible unstable modes of case TP1 (with wavenumbers $\unicode[STIX]{x1D6FC}=2m\unicode[STIX]{x03C0}/L$ , $m=1,2,3...$ ), without any phase differences. Here the relative phase between the peak-modes is modulated naturally, due to the excitation of the longest wavelengths in the domain. (d) Result of a similar perturbation as that in (c), but for case TP2 (cf. panel 2 b) wherein the phase difference between peak-mode perturbations does not affect the pattern appreciably.

To demonstrate this, we consider an extended domain of $L=4\times 2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{L}$ for case TP1 (cf. panel 2 a). This domain is four times that of figure 4. Both peak-modes are perturbed with equal amplitudes ( $A_{0}=B_{0}=0.001$ in (4.1)), but with a spatially varying phase difference $\unicode[STIX]{x1D719}(x)$ such that the peak-modes are in phase on the left half of the domain and out of phase on the right half. The function $\unicode[STIX]{x1D719}(x)$ is plotted in panel 9(a) and its exact form is given below:

The interfacial pattern that results from this perturbation is shown in panel 9(b). The left half of the domain resembles the pattern in panel 4(c), wherein the bottom interface approaches rupture at the wall. The right half of the domain resembles panel 4(d), wherein the two interfaces approach pinch off. The patterns transition from one to the other at the centre of the domain. The shorter-wave peak dominates on the left half (cf. panel 4 a) while the signature of the longer-wave peak is apparent on the right half (cf. panel 4 b), in accordance with the influence of the phase difference on inter-peak competition (cf. § 5).

That such modulation of the local phase difference between peak-modes can arise naturally, due to the influence of very-long-wavelength modes, is demonstrated by panel 9(c). Here all admissible wavelengths ( $\unicode[STIX]{x1D6FC}=2m\unicode[STIX]{x03C0}/L=m\unicode[STIX]{x1D6FC}_{L}/4$ , with $m=1,2,3\ldots$ ), with positive linear growth rates, are excited equally and in phase. The signature of the longer-wave peak-mode can be seen in the centre of the domain while the shorter-wave peak-mode dominates near the two ends of the periodic domain.

Panel 9(d) shows the results of a similar calculation for case TP2 (cf. panel 2 b) where the phase difference between peak-mode perturbations does not appreciably affect the pattern (cf. panel 5 c,d). Here, although the interface is less deformed in certain regions of the domain, the shorter-wave peak-mode (panel 5 b) dominates throughout and there is no signature of the longer-wave peak-mode (panel 5 a). These results show that coexisting patterns are likely to arise for cases wherein the phase difference between peak-modes significantly influences their interaction.

It should be noted that systems undergoing supercritical bifurcations also exhibit coexisting patterns, which are called domain structures (Raitt & Riecke Reference Raitt and Riecke1995, Reference Raitt and Riecke1997). This was demonstrated by Raitt & Riecke (Reference Raitt and Riecke1995) for a fourth-order (in space) Ginzburg–Landau equation that has a two-peak linear dispersion curve. Here the peak-modes saturate to nonlinear states with distinct wavelengths, which exist alongside one another in different regions of the domain. In the absence of a rupture-like catastrophe, these systems can evolve indefinitely. Therefore, for a coexisting pattern to be observed it is necessary for the state with coexisting patterns to be stable to perturbations. Otherwise, one pattern would invade the other and lead to a uniformly patterned state. This issue does not arise in the interfacial system under study, because rupture/pinch off occurs before the patterns arising in different regions can influence each other significantly.