1. Introduction
Falling liquid films are thin viscous liquid flows driven down an inclined plate by gravity. For several decades, they have been under extensive investigation due to their important role in various mechanical and chemical applications, such as complex flow coating (Weinstein & Ruschak Reference Weinstein and Ruschak2004) and cooling of mechanical and electronic systems (Hu et al. Reference Hu, Yang, Yu and DAI2014). For instance, the wavy nature of falling films organised into localised structures or solitary waves (Miyara Reference Miyara2000) were found to increase the heat/mass transfer rate by 
$10\unicode{x2013}100\,\%$ compared with flat films (Goren & Mani Reference Goren and Mani1968). Therefore, understanding the mechanisms leading to the complex dynamics of falling liquid films are crucial to the development of various industrial applications. Studying the linear stability is considered the first step in understanding the disordered spatiotemporal state that characterises falling liquid films. Investigating the onset of the instability and the domains in which the flow is stable or unstable is crucial to understand this kind of flows.
Mostly, the stability of falling liquid films with respect to infinitesimal perturbations is studied using the modal approach, in which the eigenmodes of an operator describing the evolution of the perturbations are examined separately. This approach relies on decomposing the flow quantities 
${{\boldsymbol q}}$ into a steady part 
$\bar {\boldsymbol q}$ (usually called base state) and an unsteady part 
$\tilde {{\boldsymbol q}}$ (perturbations). Traditionally, the base state of a falling film is found using the parallel flow approximation, where the wall-normal velocity component and the spanwise derivatives are assumed to be negligible. With this approximation, the base state forms a semiparabolic velocity profile varying in the wall-normal direction, while it remains constant in the other directions. This base state was first obtained by Nusselt (Reference Nusselt1916), and is referred to as the Nusselt film solution.
In the case of laterally confined falling films, the base state does not only change in the direction normal to the bottom wall, but also in the spanwise direction. Scholle & Aksel (Reference Scholle and Aksel2001) followed by Haas, Pollak & Aksel (Reference Haas, Pollak and Aksel2011) used analytical and experimental methods to investigate the effect of side walls on the base state of falling films. They found two competing effects in the vicinity of the wall: (i) the reduction of the base flow velocity due to the no-slip boundary condition; and (ii) an increase in the surface velocity near the side walls due to the increased film thickness caused by the capillary elevation. However, the latter effect can be ignored if the film thickness is sufficiently large.
The linear stability analysis based on the Nusselt base state is a local stability analysis, since the flow stability at a single point in the downstream is sought. The stability analysis in which the base flow varies in two spatial directions, as in laterally confined falling films, is named biglobal stability analysis. Finally, when the base state is not homogeneous in any of the considered spatial directions, the stability analysis is referred to as global stability analysis; see, for example, Albert, Tezuka & Bothe (Reference Albert, Tezuka and Bothe2014). A quantitative comparison between different stability modal approaches can be found in the work of Juniper, Hanifiand & Theofilis (Reference Juniper, Hanifiand and Theofilis2014).
The pioneering works of Kapitza & Kapitza (Reference Kapitza and Kapitza1948, Reference Kapitza and Kapitza1949) were the first to describe the development of ‘long’ wavelength deformations on the surface of the film. They identified these deformations as hydrodynamic instability, and therefore called it H-mode instability. The Orr–Sommerfeld approach has been followed in much of the work devoted to studying the local stability of falling films with respect to the hydrodynamic instability. Benjamin (Reference Benjamin1957) and Yih (Reference Yih1963) were first to follow this approach as they performed a temporal stability analysis in which the wave velocity and amplification factor for a perturbation with a given wavelength is sought. They identified the onset of the H-mode for small wavenumbers in terms of the critical Reynolds number and inclination angle 
$\beta$ as 
$Re_c = (5/6) \cot (\beta )$, where this condition is only valid when the Reynolds number is based on the average liquid velocity. On the other hand, Krantz & Owens (Reference Krantz and Owens1973) performed a spatial stability analysis where the spatial amplification and the wavelength of a disturbance with a prescribed frequency is sought. A more simplified approach was formulated by Benney (Reference Benney1966), which used long-wave expansion to obtain a temporal evolution equation (Benney's equation) that describes the interface development for small wavenumber perturbations. The evolution equation can be utilised to obtain a dispersion relation of the temporal amplification rate in the regions of small wavenumbers and weak inertia. Various physical effects were considered while studying the local stability of falling films such as heating and phase change (Mohamed & Biancofiore Reference Mohamed and Biancofiore2020), shearing gas (Lavalle et al. Reference Lavalle, Li, Mergui, Grenier and Dietze2019) and flexible inclined surface (Alexander, Kirk & Papageorgiou Reference Alexander, Kirk and Papageorgiou2020). For the wide range of mathematical approaches and physical effects considered while studying the interface evolution and local stability of falling films, refer to the reviews by Oron, Davis & Bankoff (Reference Oron, Davis and Bankoff1997), Craster & Matar (Reference Craster and Matar2009) and Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012).
With regards to the linear stability of laterally confined falling films (biglobal stability), Vlachogiannis et al. (Reference Vlachogiannis, Samandas, Leontidis and Bontozoglou2010) and Georgantaki et al. (Reference Georgantaki, Vatteville, Vlachogiannis and Bontozoglou2011) experimentally studied long-wave perturbations in narrow channels. The earlier showed that the critical Reynolds number (
$Re_c$) increases significantly as the channel width decreases, while it converges to the two-dimensional (2-D) theoretical value as channel width increases. The latter, Georgantaki et al. (Reference Georgantaki, Vatteville, Vlachogiannis and Bontozoglou2011), confirmed the aforementioned observation, and additionally showed that surface tension, opposed to classical stability, has a drastic effect on the long-wave instability threshold. They showed that the deviation in the critical Reynolds number due to spanwise confinement is strongly enhanced as the ratio between capillary and viscous forces increases. Pollak, Haas & Aksel (Reference Pollak, Haas and Aksel2011) investigated the influence of the contact angle at the side walls on the onset of the stability. Experiments showed a significant change in the neutral stability curve between two different contact angles, but with no structural modification on the neutral curves. Most recently, Kögel & Aksel (Reference Kögel and Aksel2020) performed experiments for a wide range of forcing frequencies and Reynolds numbers. They showed that the side walls can have a dramatic effect on the shape of the neutral curve by causing a fragmentation due a selective damping of the instability at moderate frequency range. They concluded that the damping effect is independent of the inclination angle and only depends on the channel width.
While the effect of spanwise confinement on the stability of falling films is not properly investigated from the theoretical perspective, there are several examples in the literature for other flow configurations. Tatsumi & Yoshimura (Reference Tatsumi and Yoshimura1990) formulated two Orr–Sommerfeld-like equations to study the stability of laminar flow in a rectangular duct. Their stability analysis showed that the critical Reynolds number increases monotonically with decreasing the aspect ratio of the duct, thus a stabilising effect is observed. Theofilis, Duck & Owen (Reference Theofilis, Duck and Owen2004) showed similar results solving the full linearised Navier–Stokes equations. More recently, Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018) and Ducimetière et al. (Reference Ducimetière, Gallaire, Lefauve and Colm-Cille2021) studied the effect of spanwise confinement on stratified shear instabilities showing that the lateral confinement has a stabilising effect.
We examine the effect of side walls on the linear stability of falling liquid films by combining theoretical modelling and experimental techniques. This approach to the problem introduces the missing theoretical background for the spanwise confinement effect reported in the literature, which is limited to experiments only. In addition, the theoretical stability model offers the opportunity to explore a wide range of flow parameters, which was not possible due to the experimental techniques restrictions. This monograph is organised as follows: in § 2, we discuss the non-dimensional governing equations alongside the dimensionless parameters. The base state and the linear stability analysis are also presented in the same section. The experimental set-up and the linear stability measurement technique are presented in § 3. Our results consisting of validating the numerical stability model with the experiments followed by extensive theoretical investigation of different aspects of the problem are presented in § 4. Finally, we present our concluding remarks and potential direction of investigations in §§ 5 and 6, respectively.
2. Theoretical formulation
The theoretical formulation of the problem is presented in this section. We start with listing the non-dimensional governing equations, and the associated walls and free surface boundary conditions in § 2.1. Afterwards, we present the base state solution (§ 2.2) followed by the methodology to linear stability analysis (§ 2.3). Finally, the numerical approach to solve the problem is presented in § 2.4.
2.1. Non-dimensional governing equations:
We consider a three-dimensional (3-D) liquid film running down a tilted channel due to gravity 
$g$ (figure 1). The channel has a width 
$W$ and forms an angle 
$\beta$ with the horizontal. The density 
$\rho$ and dynamic viscosity 
$\mu$ are constant, while the kinematic viscosity 
$\nu$ is given by 
$\nu = \mu /\rho$. We use the Cartesian coordinates 
$(x,y,z)$, where 
$x$ corresponds to the streamwise direction, 
$y$ is bottom wall normal increasing into the liquid and 
$z$ is in the spanwise direction. The left and right walls of the channel are located at 
$z = 0$ and 
$z = W$, respectively. The following length, time, velocity and pressure scales are introduced in order to non-dimensionalise the system:
\[ \begin{array}{lcccc} & h\rightarrow h_m h^*, & (x,y,z)
\rightarrow h_m (x^*,y^*,z^*), & W\rightarrow h_m
W^*,\\ & t \rightarrow h_m^2/ \nu t^*, & (u,v,w)
\rightarrow \nu /h_m (u^*,v^*,w^*), & p \rightarrow \rho
\nu^2 / h_m^2 p^*,\end{array}
\]
where 
$h_m$ is the mean film thickness. By using these scales and dropping the stars for simplicity, we obtain the dimensionless governing equations, continuity and Navier–Stokes, as follows:
\begin{gather} \partial_x {u} + \partial_y {v} + \partial_z {w} = 0, \end{gather}
\begin{gather}\partial_t u +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {u} ={-} \partial_x p + \nabla^2 {u} + G \sin(\beta), \end{gather}
\begin{gather}\partial_t v +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {v} ={-} \partial_y p + \nabla^2 {v} + G \cos(\beta), \end{gather}
\begin{gather}\partial_t w +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {w} ={-} \partial_z p + \nabla^2 {w}, \end{gather}
where 
$G = g h_m^3/\nu ^2$ is the ratio of inertial to viscous forces. The Reynolds number is a function of the flow rate 
$\dot {q}$ and channel width 
$W$:
\begin{equation} Re = \frac{\dot{q}}{2 \nu W}. \end{equation}
Figure 1. Schematic diagram of liquid film falling down an inclined channel. Here, 
$h(x,z,t)$ is the local film thickness and 
$h_m$ is the mean film thickness.
It is important to clarify that 
$G \sin (\beta )$ does not represent the Reynolds number as in the case of the 2-D problem (
$W \rightarrow \infty$). The quantity 
$G \sin (\beta )$ is always larger than the Reynolds number because the velocity is slowed down in the vicinity of the wall. Here, 
$Re$ approaches 
$G \sin (\beta )$ as the channel width goes to infinity. The dimensionless boundary conditions at the bottom 
$(y=0)$ and side walls 
$(z=0,W)$ read
\begin{gather} \boldsymbol{u}(x,y=0,z) = 0, \end{gather}
\begin{gather}\boldsymbol{u}(x,y,z = (0, W)) = 0. \end{gather}
Moreover, the kinematic and dynamic couplings at the free surface lead to the interface boundary conditions at 
$y=h(x,z,t)$ which are written as
\begin{gather} v = \partial_t h + u \partial_x h + w \partial_z, \end{gather}
\begin{align} p & = \frac{2}{n^2} [ (\partial_x h)^2 \partial_x u + (\partial_z h)^2 \partial_z w + \partial_x h \partial_z h (\partial_z u + \partial_x w) \nonumber\\ & \quad - \partial_x h (\partial_y u + \partial_x v)- \partial_z h (\partial_z v + \partial_y w) + \partial_y v] \nonumber\\ & \quad - \frac{1}{n^3} 3S [ \partial_{xx} h(1+(\partial_z h)^2) + \partial_{zz} \small(1+(\partial_x h)^2) - 2\partial_x h \partial_z h \partial_{xz} h ] , \end{align}
\begin{align} 0 & = \frac{1}{n} [ 2\partial_x h (\partial_y v - \partial_x u) + (1-(\partial_x h)^2)(\partial_y u + \partial_x v) \nonumber\\ & \quad -\partial_z h (\partial_z u + \partial_x w) - \partial_x h \partial_z h(\partial_z v + \partial_y w), \end{align}
\begin{align} 0 & = \frac{1}{n} [ 2\partial_z h (\partial_y v - \partial_z w) + (1-(\partial_z h)^2)(\partial_y w + \partial_z v) \nonumber\\ & \quad - \partial_x h (\partial_z u + \partial_x w) -\partial_x h \partial_z h(\partial_y u + \partial_x v), \end{align}
where 
$n = (1+(\partial _x h)^2 + (\partial _z h)^2)^{1/2}$. The kinematic interface boundary condition (2.4a) governs the relationship between the film thickness and the normal velocity component, while the equilibrium between the pressure and surface tension at the interface is governed by the normal and tangential stress boundary conditions (2.4b)–(2.4d), where the non-dimensional surface tension is written as 
$S = {\sigma h_m}/{3 \rho \nu ^2}$. For the complete derivation of the boundary conditions, see Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012).
Finally, we need to address the wetting effects which are also tied to the contact line behaviour at the side walls where solid, liquid and ambient meet. Formulating an accurate description of the intertwined issues of meniscus and contact line dynamics is a complex task, and has received a considerable amount of attention in the literature (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). For a thin liquid film, wetting effects result in an increase in the local film thickness at the side walls. This increase is defined as 
${\rm \Delta} h = \ell _c \sqrt {2 (1- \sin (\theta ))}$ where 
$\ell _c$ is the general capillary length (
$\ell _c = \sqrt {\sigma /\rho g \cos (\beta )}$). In the scope of this work, the Bond number (
${Bo} = W^2/\ell _c^2$), which measures the magnitude of gravitational forces to surface tension for interface movement, is very high since the capillary length is much smaller than the channel width (
$\ell _c\ll W$). Therefore, wetting effects are not dominant and are limited to a thin region near the side walls and, thus, neglected. Consequently, we implement a simple free-end boundary condition where the contact line can freely slip at the side walls with a contact angle usually chosen as 
${\rm \pi} /2$:
\begin{equation} \partial_z {h}(y=h,z=(0,W)) ={\pm} \cot(\theta) = 0. \end{equation}
For a follow-up discussion on wetting effects and contact line dynamics, the reader is referred to section (4.3). In addition, wetting effects can cause an overshoot in the velocity in the vicinity of the side walls when the increase in the film thickness 
${\rm \Delta} h$ exceeds the mean film thickness 
$h_m$ as shown by Scholle & Aksel (Reference Scholle and Aksel2001). However, this condition is never met in our work and, therefore, this effect is outside the scope our analysis.
2.2. Base flow
The effect of side walls on the base flow with an undisturbed free surface is examined in this section. When the spanwise direction is unbounded (
$W \rightarrow \infty$), the system presented in (2.1) alongside the walls and interface boundary conditions (2.3)–(2.5) has a one-dimensional (1-D) base flow solution as a semiparabolic velocity profile known as the Nusselt film solution (Nusselt Reference Nusselt1916). When the domain is confined in the spanwise direction, the base flow solution is altered because of (i) the no-slip boundary condition on the side walls which causes the velocity to decrease in the vicinity of the side walls and (ii) the wetting effects causing a capillary elevations at the side walls. Regarding the latter effect, we neglect it by setting the contact angle to 
${\rm \pi} /2$ based on the discussion in § 2.1.
Next, by assuming the bottom wall-normal velocity (
$v$) and the streamwise spatial derivatives (
$\partial _x$) to be zero, the base flow of the 3-D problem can be found numerically by solving the following system of equations:
\begin{gather} \partial_{yy} \bar{U} + \partial_{zz}\bar{U} = G \sin(\beta), \end{gather}
\begin{gather}\partial_z \bar{P} = 0. \end{gather}With the walls and interface boundary conditions:
\begin{gather} \bar{U}(y=0,z) = 0, \end{gather}
\begin{gather}\bar{U}(y,z=(0,W)) = 0, \quad \mbox{and} \quad \partial_z \bar{h}(z= (0, W)) = 0, \end{gather}
\begin{gather}\partial_y \bar{U}(y=h,z) = 0, \quad \mbox{and} \quad \bar{P}(y=h,z) = 0. \end{gather} Figure 2(a) shows the base state velocity for a confined channel with 
$W=20$. The colour scale presents the percentage drop in the velocity profile 
$\bar {U}(z,y)$ from the 
$1\textrm{-}{\rm D}$ Nusselt velocity profile 
$\bar {U}_{1\textrm{-}D}(y)$ along the 
$y$-direction. Clearly, the base state is only affected in a small region near the side walls, in which the velocity is slowed down, while it retains the Nusselt film solution in most of the channel width. We approximate this region to be in the order of the mean film thickness 
$h_m$. This is expected since the side walls effect on the velocity is mainly controlled by the ratio between the film thickness and the channel width which is typically small, for example, for a relatively narrow channel of 
$W=20$, this ratio equals 
$0.05$.
Figure 2. (a) Base state velocity for a channel with 
$W=20$ for 
$G \sin (\beta ) = 40$. (b) Normalised local flow rate per unit width 
$\bar {q}$.
Moreover, the local flow rate per unit width (
$\bar {q}$) along the channel width 
$W$ is shown in figure 2(b). The values are normalised with the flow rate per unit width for an infinitely wide channel (
$\bar {q}_{1\textrm{-}D}$), which is, in fact, equal to 
$G \sin (\beta )/3$. The local flow rate does not experience a significant drop even for relatively narrow channels. For example, the flow rate is only decreased by 
$5\,\%$ compared with the 1-D flow rate for a channel width of 
$W=25$. The local flow rate asymptotically reaches 
$\bar {q}_{1\textrm{-}D}$ as the channel width increases. This is again because the effect of side walls on the flow is restricted to a thin region in the vicinity of the walls, while outside this region, the flow obeys the 2-D Nusselt flow.
2.3. Linear stability analysis
The flow field variables 
$\boldsymbol {q}(x,y,z,t) =(u,v,w,p$) in addition to the interface variable 
$h(x,z,t)$ are expanded as a sum of the base flow and an infinitesimal perturbation as
\begin{gather} \boldsymbol{q}(x,y,z,t) = \bar{\boldsymbol{q}}(y,z,t) + \epsilon \tilde{\boldsymbol{q}}(x,y,z,t), \end{gather}
\begin{gather}h(x,z,t) = 1 + \epsilon \tilde{h}(x,z,t). \end{gather} Expansions in (2.8) are then substituted in the governing equations (2.1) and boundary conditions (2.3)–(2.5), which are then linearised for 
$\epsilon \ll 1$ leading to the linearised perturbation equations at 
$O (\epsilon )$:
\begin{gather} \partial_x \tilde{u} + \partial_y \tilde{v} + \partial_z \tilde{w} = 0, \end{gather}
\begin{gather}\partial_t \tilde{u} + \bar{U} \partial_x \tilde{u} + \partial_y \bar{U} \tilde{v} + \partial_z \bar{U} \tilde{w} + \partial_x \tilde{p} - \nabla^2 \tilde{u} = 0, \end{gather}
\begin{gather}\partial_t \tilde{v} + \bar{U} \partial_x \tilde{v} + \partial_y \tilde{p} - \nabla^2 \tilde{v} = 0, \end{gather}
\begin{gather}\partial_t \tilde{w} + \bar{U} \partial_x \tilde{w} + \partial_z \tilde{p} - \nabla^2 \tilde{w} = 0, \end{gather}with the bottom and side wall boundary conditions
\begin{gather} \tilde{\boldsymbol{u}}(x,y=0,z) = 0, \end{gather}
\begin{gather}\tilde{\boldsymbol{u}}(x,y,z=(0, W)) = 0. \end{gather}
For the interface boundary conditions at 
$y= 1+\tilde {h}$, we utilise a Taylor expansion and expand the equations around the undeformed interface 
$y=1$, for example the variable 
$X$ is expanded as follows 
$X|_h = X(1) + \tilde {x}|_1 + DX(1) \tilde {h}$:
\begin{gather} \tilde{v} - \partial_t \tilde{h} - \bar{U} \partial_x \tilde{h} = 0, \end{gather}
\begin{gather}\tilde{p} + \partial_y \bar{P} \tilde{h} - 2 \partial_y \tilde{v} + 3S \nabla^2_{xz} \tilde{h} = 0, \end{gather}
\begin{gather}- \partial_{yy} \bar{U} \tilde{h} - \partial_y \tilde{u} - \partial_x \tilde{v} + \partial_z \tilde{h} \partial_z \bar{U} = 0, \end{gather}
\begin{gather}\partial_y \tilde{w} + \partial_z \tilde{v} - \partial_z \bar{U} \partial_x \tilde{h} =0. \end{gather}Afterwards, the pressure perturbation boundary conditions are needed at the walls in order to close the problem. We use the solution proposed by Theofilis et al. (Reference Theofilis, Duck and Owen2004), which is to introduce the compatibility condition at the walls derived from the Navier–Stokes equations as follows:
\begin{gather} \partial_y \tilde{p} = \nabla^2 \tilde{v} - \bar{U} \partial_x \tilde{v}, \end{gather}
\begin{gather}\partial_z \tilde{p} = \nabla^2 \tilde{w} - \bar{U} \partial_x \tilde{w}. \end{gather} Finally, we introduce Fourier modes to expand the perturbation variables in 
$x$ and 
$t$ as
\begin{gather} \tilde{\boldsymbol{q}}(x,y,z,t) = \boldsymbol{Q}(z,y) \exp({\rm i}kx - {\rm i}\omega t), \end{gather}
\begin{gather}\tilde{h}(x,z,t) = \eta(z) \exp ({\rm i}kx - {\rm i}\omega t), \end{gather}
where 
$k$ is the wavenumber and 
$\omega$ is the angular frequency. For temporal stability analysis, 
$k$ is considered to be real, and we solve for the complex eigenvalue 
$\omega$. The temporal growth rate 
$\omega _i$ is the imaginary part of 
$\omega$. If 
$\omega _i >0$ the perturbation grows in time and the base flow is unstable. Moreover, the phase velocity is denoted as 
$c_r = \omega _r/k$. For spatial stability analysis, 
$\omega$ is assumed real and 
$k$ is complex. Similarly, the flow is unstable if the spatial growth rate 
$\gamma _i$ is larger than zero. Equations (2.13) may describe waveguide modes (Lighthill Reference Lighthill1978), where the wave disturbances are allowed to travel in the longitudinal direction while being confined in the spanwise direction. They can be interpreted physically as plane modes (
$W \rightarrow \infty$) bouncing back and forth between the side walls (Elmore & Heald Reference Elmore and Heald1969). This concept is revisited in the next sections as it leads to interesting results.
2.4. Numerical method
We discretise the domain using 2-D Chebyshev differentiation matrices (Trefethen Reference Trefethen2000). The boundary conditions were implemented by replacing the corresponding rows in the matrices with the discretised boundary conditions. We prefer a temporal stability analysis since it requires less storage and computational power than what is required by the spatial stability problem. We base our results mainly on the temporal growth rate 
$\omega _i$, and utilise Gaster transformation (Gaster Reference Gaster1962) that allows us to relate the temporal and spatial growth rates as
\begin{equation} \gamma = \omega_i \left/ \frac{\partial \omega_r}{\partial k} \right. .\end{equation}Brevdo et al. (Reference Brevdo, Laure, Dias and Bridges1999) showed that this transformation is applicable to falling liquid films since the spatial and temporal growth rates are small. Finally, the eigenvalue problem was solved using the QZ algorithm (Moler & Stewart Reference Moler and Stewart1973).
3. Experimental set-up
We use the very experimental rig as (Pollak et al. Reference Pollak, Haas and Aksel2011; Schörner, Reck & Aksel Reference Schörner, Reck and Aksel2016; Kögel & Aksel Reference Kögel and Aksel2020). Here, we briefly present the experimental apparatus and measurement technique. For more details, please refer to the aforementioned works.
3.1. Experimental apparatus
The experimental set-up we use in this work is sketched in figure 3. The set-up consists of a channel with an aluminium bottom and Plexiglas side walls. The width of the channel is 
$170$ mm, but two movable side walls can be placed within the channel to achieve smaller channel widths. The whole channel is tilted with an inclination angle 
$\beta$ and mounted onto a damped table in order to reduce the effect of vibrations. A collection tank is located at the end of the channel in which the liquid is maintained at a certain temperature 
$T$. An eccentric pump is used to convey the liquid from the collection tank to the inflow tank located at the inlet of the channel. The inflow tank allows the liquid to rest and to smoothly overflow into the channel. Periodic waves are introduced in the flow using a paddle that is permanently dipped in the liquid film. The paddle oscillates at an adjustable frequency 
$f_p$ and amplitude 
$A_p$, causing sinusoidal waves evolving downstream the channel with the same frequency and amplitude. The disturbances generated by the paddle are instantaneously convected downstream the channel, and do not disturb the flow in the inflow tank. The working liquid used in this experiment is Elbesil 140 which is a mixture of silicon oils, Elbesil 50 and Elbesil 200. The physical properties of Elbesil 140 at temperature 
$24\,^{\circ }$C are listed in table 1.
Figure 3. Sketch of the experimental set-up used to measure linear stability.
Table 1. Physical properties of Elbesil 140 at working temperature 
$T = 24\,^{\circ }\mbox {C}$.
3.2. Linear stability measurement
The linear stability of the flow was investigated through measuring the spatial growth rate of the excited waves 
$\gamma$. Two laser beams are reflected at the film's surface at two locations, upstream and downstream. The two laser reflections are projected onto a screen on which they oscillate at the same frequency of the sinusoidal waves generated by the paddle. The oscillations of the laser spots on the screen can can be written as
\begin{gather} A_u(t) = A_u^0 \sin(2{\rm \pi} f t), \end{gather}
\begin{gather}A_d(t) = A_d^0 \sin(2{\rm \pi} f t), \end{gather}
where 
$A^0$ is the amplitude of the signal, and the subscripts 
$u$ and 
$d$ correspond to upstream and downstream locations, respectively. Simply, the oscillation amplitudes 
$A_u^0$ and 
$A_d^0$ are proportional to the oscillation amplitude of the liquid interface slope, which is also proportional to the waves amplitudes. This is true as long as the waves are linear and sinusoidal, for a detailed discussion please see the work of Pollak et al. (Reference Pollak, Haas and Aksel2011). A CCD camera with a 100 frames per second is used to capture the screen for 10 seconds, resulting in 1000 images each measurement. Figure 4(a,b) show the oscillation of the upstream and downstream laser spots, for a paddle frequency 
$f_p$ equals 
$10$ Hz with an amplitude 
$A_p= 0.6$ mm. The oscillations of the laser spots were obtained by applying a Gaussian filter on the images and tracing the centres of the laser spots. The mean film thickness 
$h_m$ was measured to be 
$9$ mm and, thus, the excited waves amplitude is 
$6.6\,\%$ of the film thickness. The excited waves amplitude was chosen large enough to be adequately examined, but small enough to be compatible with the linear stability theory. Subsequently, the amplitudes 
$A_u^0$ and 
$A_d^0$ were found by performing discrete Fourier analysis of the laser spots oscillations traced from the images. Figure 4(c) shows the Fourier spectra for the same paddle frequency and amplitude (
$f_p=10, A_p=0.6$). The spectra show a dominant peak indicating that the signal is monochromatic with a frequency very close to the excitation frequency (maximum difference between the excitation frequency and response frequency is 
$0.1$ Hz).
Figure 4. Normalised upstream and downstream signals and their corresponding discrete Fourier transform for forcing frequency 
$f = 10$ Hz and amplitude 
$A_p = 0.6$ mm. The peaks are located at 
$f=9.91$ Hz.
The downstream amplitude was geometrically corrected since it covers a smaller distance from the reflection point to the screen as follows:
\begin{equation} A^0_{d,{corr}} = A^0_d \left(1 + \frac{{\rm \Delta} x}{x_S}\right), \end{equation}
where 
${\rm \Delta} x = 400$ mm is the distance between the two reflection points, and 
$x_S = 2.966$ m is the distance from the downstream reflection point to the screen. Furthermore, 
$x_1$ is the distance between the channel inlet and upstream laser spot, which we made sure to be sufficiently large to avoid inlet effects on the base flow. Finally, the linear growth rate of the excited waves is obtained as
\begin{equation} \gamma = \frac{\ln(A^0_{d,{corr}} /A_u^0 )}{{\rm \Delta} x}. \end{equation}4. Results
We start with confronting our numerical linear stability with experimental results in § 4.1. Afterwards, in § 4.2, we utilise the numerical model to examine the effect of the spanwise confinement on several aspects of the problem, such as (i) maximum linear growth rate, (ii) types of stability modes and their distinct behaviour, (iii) eigenmodes structure and (iv) eigenvalues spectra.
4.1. Comparison of experiments and theory
In this section, we compare the results obtained by our numerical stability model, our experiments, as well as the experiments conducted by Kögel & Aksel (Reference Kögel and Aksel2020). The comparison is based on the spatial growth rate 
$\gamma$, which is directly found from the experimental measurements as explained in § 3.2. Gaster transformation is utilised to obtain the numerical spatial growth rate from the quantities 
$\omega$ and 
$k$ as shown in (2.14). We also need to obtain the response frequency 
$f$ corresponding to each wavenumber 
$k$, which can be easily found as 
$f=\omega _r(k)/2{\rm \pi}$. The experimental control parameters and their corresponding non-dimensional numbers used in the numerical model are listed in table 2, in which the general capillary length is shown to be much smaller than the channel width (
${Bo} \approx 3850$) and, therefore, the wetting effects can be neglected.
Table 2. Experiment parameters and non-dimensional numbers for the numerical stability model.
Figure 5 shows the spatial growth rate in terms of the response frequency. In general, good agreement is observed between our numerical model and our experiments. Most importantly, the sharp cusp in the growth rate and its corresponding frequency are captured accurately. The minimal shift in the response frequency is due to the accuracy of the Fourier analysis which has an error of 
$\pm 0.1$ Hz. The grey band represents the upper and lower limits found by the numerical model. These limits were obtained based on the tolerances in the fluid properties and experiment control parameters (tables 1 and 2). For example, the upper limit was obtained by choosing the parameters limits that result in the maximum growth rate combining 
$(G_{max},S_{min},\beta _{max},W_{max})$, while the opposite was applied in finding the lower limit.
Figure 5. Spatial growth rate obtained using our theoretical model, our experiments and the experiments of Kögel & Aksel (Reference Kögel and Aksel2020). The colour in background represents the amplitude of 
$\sup _{f_r \in \mathbb{R} } \|\boldsymbol{\mathsf{R}}\|$.
Furthermore, while our theoretical and experimental results match very well for low and high frequencies, an evident discrepancy is observed before and after the cusp around the frequencies 
$5.5$ and 
$7.5$ Hz. We offer a heuristic support for this discrepancy by plotting the supremum of the resolvent norm (
$\sup _{f_r \in \mathbb {R} } \|\boldsymbol{\mathsf{R}}\|$), which indicates the maximum amplification over the excitation frequency 
$f_r$ (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993):
\begin{equation} \sup_{f_r \in \mathbb{R} } \|\boldsymbol{\mathsf{R}}\| = \sup_{f_r \in \mathbb{R} } \| ({-}i f_r \boldsymbol{\mathsf{B}}+ \boldsymbol{\mathsf{A}})^{{-}1}\|, \end{equation}
where 
$\boldsymbol{\mathsf{A}}$ and 
$\boldsymbol{\mathsf{B}}$ are our matrices in the general eigenvalue problem. Since the streamwise spatial derivative is replaced by 
$ik$ based on the modal decomposition in (2.13), this requires sweeping over the frequency 
$f_r$ to find 
$\sup _{f_r \in \mathbb {R} } \|\boldsymbol{\mathsf{R}}\|$ at a certain wavenumber 
$k$. Finally, we find the corresponding response frequency for every wavenumber to be able to compare against the spatial growth rate, as shown in figure 5. It is obvious that the norm of the resolvent appears to be proportional to the discrepancy between the numerical and experimental results. The maximum value of 
$\sup _{f_r \in \mathbb {R} } \|\boldsymbol{\mathsf{R}}\|$ is located before the cusp at which the maximum discrepancy exists ( 
$f \approx 5.5$ Hz), while the lower peak matches the smaller discrepancy region after the cusp ( 
$f \approx 7.5$ Hz). Therefore, the discrepancy between the theoretical and experimental growth rate appears to be maximum where 
$\sup _{f_r \in \mathbb {R} }\|\boldsymbol{\mathsf{R}}\|$ is large. While the experimental growth rate accounts for the maximum amplification of the disturbances, the numerical growth rate is taken directly from the most unstable stability mode.
4.2. Theoretical results
We now study the effect of spanwise confinement on the linear stability from different perspectives. The numerical model allows exploring the parameters space, which was not possible using experimental means. To begin with, it is useful to recall the results of the well-studied 1-D local stability problem (where the spanwise variation is not taken into account) in order to accurately assess the influence of the side walls. The 1-D problem results in only one prominent stability mode, which is responsible for the 1-D hydrodynamic instability. For this reason, is it called the H-mode. (We refer to 1-D H-mode for distinction as 1-D-mode in this paper and preserve the name H-mode for the 2-D equivalent.) We found that our problem results in several stability modes which behave differently depending on the degree of the spanwise confinement. In the following, we first focus on the maximum temporal growth rate regardless of which stability mode is dominant. Afterwards, we quantify the most important stability modes and examine them from different perspectives, in order to find an explanation for the unique effect of the side walls on the stability.
4.2.1. Initial observations: maximum temporal growth rate
We examine the effect of spanwise confinement on the maximum temporal growth rate for a certain range of parameters. Figure 6 shows the temporal growth rate contours as a function of 
$G \sin (\beta )$ and 
$k$ for different 
$W$ and 
$\beta$ values. The white area is stable, while the coloured area is unstable. The solid line is a local minimum in the growth rate along the wavenumber, while the dashed line is the neutral curve for the classical 1-D-mode. In general, we observe one unstable region identical to that of the 1-D problem for 
$W = \infty$. As spanwise confinement becomes stronger, significant modifications are seen on both, the structure and the magnitude of the unstable region, while the local minimum in the growth rate shifts to higher wavenumbers. Moreover, the stabilisation effect shows to be stronger for small inclination angles as seen when comparing the unstable regions of the inclination angles 
$\beta =5$ and 
$\beta =15$.
Figure 6. Temporal growth rate contours in the 
$G\sin (\beta )-k$ space for (a–d) 
$\beta =5$, (e–h) 
$\beta =10$ and (i–l) 
$\beta =15$, with 
$S=1$. The dashed line is the neutral curve of the 1-D problem, while the purple line is a local minimum along the wavenumber.
A local minimum in the growth rate is shown at 
$W=30$ for inclination angle 
$\beta =10$. Decreasing the channel width to film height ratio to 
$W=20$, shifts the local minimum to higher wavenumbers, while the unstable region starts splitting. Further decrease in the channel width creates two unstable regions with a horizontal stable band in between. For a smaller inclination angle (
$\beta =5$), the stabilisation effect is strong enough that the upper unstable region is completely stabilised, while the local minimum still exists. On the other hand, for 
$\beta =15$, the stabilisation effect is less dominant; however, it manages to create the two unstable regions with a much thinner stable band separating them. Note that the mean film thickness (
$h_m$) is kept constant while changing the channel width 
$W$, this is done by changing the flow rate 
$\dot {q}$ accordingly in order to keep the Reynolds number 
$Re$ constant, as (2.2) shows.
A more detailed presentation of the stabilisation behaviour can be shown by plotting the maximum temporal growth rate against the wavenumber for a fixed 
$G\sin (\beta )$. In figure 7(a), the growth rate at intermediate wavenumbers is affected the most, while it approaches the 1-D growth rate at small and large wavenumbers. This agreement decreases as the confinement becomes stronger. A smooth dent in the growth rate curve is seen for moderate confinement (
$W=50$), which slowly transforms into a sharp cusp as the channel becomes narrower. Furthermore, figure 7(b) highlights the effect of the inclination angle on the stabilisation behaviour. It is clear that the flow is stabilised more significantly at smaller inclination angles when comparing the difference between the growth rates for 
$W=30$ (solid) and that of the 1-D problem (dashed). In addition, the band of stabilised wavenumbers and the local minimum are located at higher values for smaller inclination angles. To conclude this section, the fragmentation effect of the unstable region was first observed in the experimental results of Kögel & Aksel (Reference Kögel and Aksel2020); however, no theoretical explanation was offered for this unique stabilisation behaviour. We suggest an explanation for this phenomenon in the following sections using different perspectives of the problem.
Figure 7. (a) Maximum growth rate for different 
$W$ values for 
$G\sin (\beta ) =85$, 
$\beta = 10$ and 
$S=1$. (b) Influence of the inclination angle on the stabilisation effect for 
$W=30$, 
$G =500$ and 
$S=1$. The dashed lines correspond to the 1-D growth rate.
4.2.2. The 2-D stability modes
The initial observations in the previous section were based on the maximum temporal growth rate 
$\omega _i$, regardless of what is the dominant stability mode. In the following, we examine several stability modes, and analyse the effect of confinement on them. The modes are obtained using the arc-continuation technique (Chan & Keller Reference Chan and Keller1982).
Figure 8 shows the growth rate 
$\omega _i$ along the wavenumber 
$k$ for the most dominant stability modes for several degrees of lateral confinement. Starting with very weak confinement (
$W = 500$) as in figure 8(a), all the stability modes except for two, coincide all together matching the 1-D-mode (dashed line). Thus, these modes belong to the H-mode type. The remaining pair of modes matches the H-modes for small wavenumber, but diverges at large values. We call these two modes wall modes (or simply W-modes).
Figure 8. Stability modes along the wavenumber for (a) 
$W=500$, (b) 
$W=40$, (c) 
$W=25$ and (d) 
$W=15$ for the parameters 
$G=500$, 
$\beta =10$ and 
$S=1$. H-modes are blue and W-modes are orange.
For 
$W = 40$, significant changes occur to all the stability modes as seen in figure 8(b). All the modes are stabilised in general, but we observe three distinct stabilisation actions; the most unstable H-mode is mainly stabilised at moderate wavenumbers showing the smooth indentation observed earlier in figure 7(a). More interestingly, the remaining more stable H-modes are strongly stabilised at low wavenumbers, while the opposite is experienced by the W-modes, which show strong stabilisation at large wavenumbers. Moreover, the local minimum of the most unstable H-mode (
$\mathcal {H}_1$) is accurately aligned with the maximum of the second most unstable W-mode (
$\mathcal {W}_2$). This configuration will oddly develop in later stages. From another perspective, the H-modes (except for the most unstable mode) are transformed into waveguide modes. Waveguide modes are similar to the plane 1-D mode, but allowed to propagate in the longitudinal direction only. They are also restricted to a cut-off wavenumber below which they do not propagate, which explains the stabilisation of the H-modes for low wavenumbers.
Figure 8(c) shows similar stabilisation actions experienced by the H-modes and W-modes for 
$W=25$. Intriguingly, an attraction centre involving 
$\mathcal {H}_1$ and 
$\mathcal {W}_2$ modes takes a place leading to the formation of a sharp cusp in 
$\mathcal {H}_1$ and an acute tip in 
$\mathcal {W}_2$. The reason behind this attraction involving those two specific modes among the others will be investigated in subsequent sections. Eventually, further decrease in the channel width (
$W = 15$) shows the final picture in figure 8(d). The concurrence between the 
$\mathcal {H}_1$ cusp and 
$\mathcal {W}_2$ tip sufficiently increases until they meet and switch branches for wavenumbers before the meeting point. At this stage the most unstable H-mode also becomes a waveguide mode. This mode switching behaviour was also observed previously while studying the stability of falling liquid films on flexible substrates when changing the damping ratio (Alexander et al. Reference Alexander, Kirk and Papageorgiou2020).
In summary, the solution of the biglobal stability problem contains several prominent stability modes. For the parameters under considerations, most of these modes belong to the H-mode type, while two of them are different, we name them W-modes. The stabilisation effect at moderate wavenumbers is a result of two factors. First, the spanwise confinement influence the two stability mode types in a different manner, where H-modes are transformed into waveguide modes, and stabilised at low wavenumbers, while W-modes are stabilised at high wavenumbers. Second, an attraction centre involving two different modes, the most unstable H-mode 
$(\mathcal {H}_1)$ and the second most unstable W-mode
$(\mathcal {W}_2)$. This attraction centre strengthens as confinement increases, until a mode switch takes place in the branch with lower wavenumbers, after which the most unstable H-mode becomes a waveguide mode as well.
With the same approach, figure 9 shows the phase velocity (
$c_r = \omega _r/k$) for the stability modes for 
$W=40$ and 
$W=15$. The highlighted curve follows the maximum temporal growth rate along the wavenumber. We observe a similar behaviour to that of the temporal growth rate as a transition in the phase velocity is shown leading to a sharp jump when the channel width is sufficiently small. For 
$W=40$, 
$\mathcal {H}_1$ is still the most dominant stability mode along the wavenumber, and the phase velocity is very close to that of the 1-D mode, while the other H-modes show higher phase velocity. The W-modes phase velocity drop sharply as wavenumber increases since they are stabilised strongly in this range. As channel width decreases to 
$W=15$, the mode switch has occurred and a sudden jump in the phase velocity is observed. The phase velocity of the W-modes does not change significantly, but all the H-modes show a higher phase velocity than that of the 1-D mode, which is a characteristic behaviour of waveguide modes (Elmore & Heald Reference Elmore and Heald1969).
Figure 9. Phase velocity 
$c_r$ for different stability modes along the wavenumber for (a) 
$W=40$ and (b) 
$W=15$ for the parameters 
$G=500$, 
$\beta =10$ and 
$S=1$. The highlighted curve follows the maximum temporal growth rate along the wavenumber.
4.2.3. Perturbation patterns and eigenvectors structure
In this section, we study the perturbation patterns and the eigenvectors structure of different stability modes for several confinement configurations. This will identify the structural nature of the different mode types, and help us better comprehend their particular behaviour as the channel width decreases. First, we compare the interface perturbation field 
$\tilde {h}(x,z)$ of the stability modes at two wavenumbers, i.e. 
$k=0.15$ and 
$k=0.75$. These two values are suitable options since the spanwise confinement affects the modes heavily at the former, but very lightly at the latter wavenumber. Afterwards, we plot the eigenmodes as a function of the wavenumber, along which the stability modes change significantly, in particular for strong confinement.
Figure 10 shows 
$\tilde {h}(x,z)$ for the pair of modes (
$\mathcal {H}_1, \mathcal {W}_2$) for different 
$W$ values. The pair has an even symmetry around the axis 
$z=0$, which offers an explanation for the unique behaviour involving them. For 
$W=500$, the modes show their natural structure at 
$k=0.75$, where 
$\mathcal {H}_1$ is concentrated at the centre of the channel, while 
$\mathcal {W}_2$ is only observed in the vicinity of the side walls, as expected for a mode created by the confinement. At 
$k=0.15$, the two modes are in close proximity, see figure 8(a) and, thus, they affect each other slightly causing 
$\mathcal {H}_1$ to increase its curvature, while 
$\mathcal {W}_2$ spreads more into the domain. These structures are logical since confinement is yet very weak. For 
$W=25$, 
$\mathcal {H}_1$ changes its structure to a straight wave along the spanwise direction at 
$k=0.15$, while it keeps the original structure at 
$k=0.75$ with more curvature. For the same confinement (
$W=25$), 
$\mathcal {W}_1$ is dramatically enhanced at 
$k=0.15$, while it just spreads slightly into the domain at 
$k=0.75$. Finally, 
$W=15$ shows a structure switch between the two modes at 
$k=0.15$. It also shows a strong curvature in 
$\mathcal {H}_1$ at 
$k=0.75$, while stronger spreading for 
$\mathcal {W}_1$.
Figure 10. Interface perturbation field 
$\tilde {h}(x,z)$ for the modes (a) 
$\mathcal {H}_1$ and (b) 
$\mathcal {W}_2$ at different wavenumbers for the parameters 
$G=500$, 
$\beta =10$ and 
$S=1$.
For a complete analysis, we also plot the interface perturbation field 
$\tilde {h}(x,z)$ for the modes pair (
$\mathcal {H}_2, \mathcal {W}_1$), which are the second most unstable H-mode and the most unstable W-mode, respectively (figure 11). These two modes have an odd symmetry, such as 
$\tilde {h}(x,z)=-\tilde {h}(x,-z)$, and show less complex behaviour compared with the other pair (
$\mathcal {H}_1, \mathcal {W}_2$). For 
$W=500$, the spatial structure 
$\mathcal {H}_2$ is one step higher in the harmonic range compared with 
$\mathcal {H}_1$, while 
$\mathcal {W}_1$ is restricted at the walls. As confinement becomes stronger (
$W=25, W=15$), 
$\mathcal {H}_2$ changes structure and shows higher harmonics of 
$\mathcal {W}_1$ at 
$k=0.15$, while it keeps its original structure with a curvature proportional to the confinement degree at 
$k=0.75$. Furthermore, 
$\mathcal {W}_1$ spreads more into the domain for all wavenumbers, but strongly at lower values. For both mode pairs, the symmetry of the perturbation fields are conserved for any 
$W$ value.
Figure 11. Interface perturbation field 
$\tilde {h}(x,z)$ for the modes (a) 
$\mathcal {H}_2$ and (b) 
$\mathcal {W}_1$ at different wavenumbers for the parameters 
$G=500$, 
$\beta =10$ and 
$S=1$.
With regards to the eigenmodes, figure 12 shows the absolute of the normalised 
$\eta (z)$, which is the eigenmode of the interface perturbation field 
$\tilde {h}(x,z)$ (see (2.13b)). The first row belongs to the H-mode 
$\mathcal {H}_1$, while the second row belongs to most unstable W-mode (
$\mathcal {W}_1$). The third and fourth rows show the eigenmodes 
$\mathcal {W}_2$ and 
$\mathcal {H}_2$, respectively. The dashed lines present the location at which the different mode types cross each other, except for 
$\mathcal {H}_1$ at 
$W=40$ and 
$W=25$ where it shows the local minimum. When the spanwise confinement is very weak (
$W=500$), all the eigenmodes are constant along the wavenumber except at small values. The most dominant mode 
$\mathcal {H}_1$ is concentrated at the centre and zero at the walls, while 
$\mathcal {H}_2$ has similar structure but with a higher step in the harmonic range. On the other hand, the W-modes are strictly concentrated at the side walls and zero everywhere else. The stability modes form a ‘discrete spectrum’ consisting of the H-modes only (
$\mathcal {H}_1, \mathcal {H}_2, \mathcal {H}_3$, etc.), while W-modes play no major role. This is in agreement to what is observed earlier in figures 10 and 11.
Figure 12. The real part of normalised 
$\eta (z)$ eigenmode along the wavenumber for different confinement ratios, for the parameters 
$G=500$, 
$\beta =10$ and 
$S=1$.
The unique patterns in the eigenmodes at small wavenumbers for 
$W=500$ spread into higher values as confinement increases to 
$W=40$ and 
$W=25$. The presence of the W-modes become more evident causing significant restructuring of the other eigenmodes along the wavenumber axis. The region in which the eigenmodes are reconstructed extends into higher wavenumbers proportionally to the strength of the confinement. In this limited region (bounded by the dashed line), the stability modes form a different harmonic sequence; 
$\mathcal {W}_1$ is no longer only concentrated at the walls but spreads towards the centre of the channel, while 
$\mathcal {W}_2$ shows similar structure but with the further step in the harmonic range. Moreover, the H-modes are reconstructed to fit the new harmonic sequence; 
$\mathcal {H}_2$ shows a third step in the new harmonic range, while 
$\mathcal {H}_1$ shows a constant maximum along 
$z$, which is the lowest step in the harmonic range. The eigenmodes are slightly affected after the cusp (dashed line), and the original harmonic series is conserved. For 
$W=15$, the same structures spread to higher wavenumbers and become more prominent. In addition, we notice the switch of the eigenmode between 
$\mathcal {H}_1$ and 
$\mathcal {W}_2$ in order to conserve the harmonic series.
In short, when the spanwise confinement is negligible, the H-modes form a harmonic sequence with a fundamental mode that is maximum at the centre and zero at the side walls, while W-modes are insignificant. When confinement becomes stronger, the W-modes presents an opposite (as in maximum at the walls and zero in the centre of the channel) harmonic series at small wavenumbers region. This forces a reconstruction on the H-modes to be compatible in the new harmonic sequence.
4.2.4. Eigenvalue spectra of confined 2-D problem vs 1-D oblique problem
As seen, the W-modes play no significant role for weak spanwise confinement, but experience dramatic change in both growth rate and spatial structure for strong confinement. We also observed that they seemed to developed from H-modes when confinement is weak. In this section, we confront the eigenvalues spectrum of the 2-D problem to that of the oblique unbounded case to show how the two W-modes evolve as confinement becomes stronger. In the oblique simplified problem, the base state is spanwise invariant (
$\bar {U}(y),\bar {P}(y)$), while the perturbation fields have a harmonic structure in the spanwise direction, and expanded as follows:
\begin{gather} \tilde{\boldsymbol{q}}(x,y,z,t) = \boldsymbol{\hat{Q}}(y) \exp ({\rm i} \zeta z) \exp ({\rm i}kx - {\rm i}\omega t), \end{gather}
\begin{gather}\tilde{h}(x,z,t) = \hat{\eta} \exp ({\rm i} \zeta z) \exp ({\rm i}kx - {\rm i}\omega t), \end{gather} where 
$\zeta \in \mathbb {R}$ is the spanwise wavenumber. Equation (4.2) implies that the spanwise direction is unbounded, since any spanwise wavelength is possible for the perturbations (
$\zeta \in \mathbb {R}$). For 
$\zeta =0$, the oblique problem becomes identical to the 1-D problem, while for 
$\zeta \ne 0$, a spanwise variation is allowed for the perturbations without considering the confining effect of the side walls. Similarly to the 1-D problem, only one unstable mode (H-mode) shows in the oblique problem regardless of the value of 
$\zeta$.
Figure 13 compares the eigenvalues spectrum of the 2-D problem with the oblique mode. The two problems are compared at three different wavenumbers for several confinement values. The oblique solutions in all the subfigures are the same since the problem is not a function of the spanwise confinement. For 
$\zeta =0$, the oblique mode (line) matches the 1-D mode. As 
$\zeta$ increases, the oblique modes become more stable since they experience stronger viscous damping. This is a consequence of the increase in the diffusive term weight as the Laplacian term became (
$-k^2+\partial _{yy}-\zeta ^2$) instead of (
$-k^2+\partial _{yy}+ \partial _{zz}$) due to the oblique expansion in (4.2).
Figure 13. The 2-D problem eigenvalues spectrum (hollow markers) vs oblique dispersion relation (lines) as a function of the spanwise wavenumber 
$\zeta$: (a) 
$W=100$, (b) 
$W=40$, (c) 
$W=25$, (d) 
$W=15$ at wavenumber values 
$k=0.15$ (
$\raise 3pt\hbox{$_{_\square}$}$, solid), 
$k=0.35$ (
$\circ$, dashed), 
$k=0.75$ (
$\diamond$, dotted). The red dots mark the 1-D mode, while the blue and orange markers are the 2-D H-modes and W-modes, respectively.
Focusing first on 
$W=100$ (figure 13a), the 2-D modes including both H-modes and W-modes match the oblique dispersion relation at small wavenumber 
$k=0.15$, while being tightly packed around the 1-D mode. For larger wavenumbers, the H-modes follow the oblique mode but the W-modes start departing from the line (
$k=0.35$), until they disappear at 
$k=0.75$. For stronger confinement as in figure 13(b–d), the 2-D modes spread out on the oblique line, while the agreement decreases. The most evident discrepancy between the 2-D modes and the oblique dispersion relation is observed for moderate and high wavenumbers. This discrepancy is a reason of the complex behaviour we observed involving the stability modes, which takes place at moderate wavenumbers and influence higher wavenumbers as well. Interestingly, the 2-D modes including the W-modes are very close with the oblique dispersion relation for strong confinement at small wavenumber 
$k=0.15$. Moreover, the 2-D modes are always more stable than the 1-D mode for all parameters.
4.3. Discussion on surface tension and contact line dynamics
The results presented are limited to cases where surface tension is not the dominant restoring force where 
$S \sim {O} (1)$ in both theoretical analysis and experiments. In such scenarios, the capillary length is very small compared with the channel width, for example Bo is of 
${O}(10^3)$ in experiments. Therefore, neglecting the wetting effects and using the simple free-end contact line boundary condition (2.5) is admissible. This is also consistent with the base flow under consideration in which the interface is assumed to be invariant in the spanwise direction.
Nevertheless, capillary dissipation due to wetting effects and contact line dynamics in confined surface waves is non-negligible in several cases, such as very narrow channels (Monsalve et al. Reference Monsalve, Maurel, Pagneux and Petitjeans2022), micro-gravity conditions (Zhang & Thiessen Reference Zhang and Thiessen2013) or simply when surface tension is the main restoring force, and capillary forces are dominant over viscous forces (Georgantaki et al. Reference Georgantaki, Vatteville, Vlachogiannis and Bontozoglou2011). The last scenario is directly related to this work where it was shown experimentally that surface tension, when dominant, can significantly affect the onset of the long-wave instability in falling films confined by side walls. They suggest that this effect is a result of a 3-D capillary wave attenuation mechanism, which stems from damping the oscillations near the side walls caused by (i) viscous boundary layer and, more relevantly, (ii) the heightened resistance to the depinning of the contact line.
Under the stability model in this work, the free-end boundary condition (2.5) allows the interface to freely slip at the side walls, which eliminates the damping effect due to resistance of contact line depinning. Therefore, the unique capillary damping of oscillations near the side walls, which leads to a delay in the long-wave instability onset, is not captured in our model. In order to be able to analyse scenarios in which surface tension is dominant, a more complex contact line boundary condition is necessary. Hocking (Reference Hocking1987), motivated by experimental evidence suggesting that contact angle is not fixed but depends on contact line motion (Dussan Reference Dussan1979), proposed a contact line boundary condition in which the contact line velocity is proportional to the variation in the contact angle 
$\theta$:
\begin{equation} \frac{\partial h}{\partial t} = \alpha \frac{\partial h}{\partial z}, \end{equation}
where 
$\alpha$ is a slip coefficient that characterises the boundary condition. An important feature of this boundary condition is that it includes two extreme limits of zero and infinite slip coefficients, which correspond to a pinned-end (
$\partial _t h =0$) and a free-end (
$\partial _z h =0$) conditions, respectively. More complex models are present in the literature that incorporate contact angle hysteresis, where 
$\alpha$ is described as function of contact line displacement (Viola, Brun & Gallaire Reference Viola, Brun and Gallaire2018).
Another related issue to be addressed here is the well-known contact line singularity, which results from assuming a moving contact line at the wall alongside the no-slip boundary condition (Dussan & Davis Reference Dussan and Davis1974). This singularity can be avoided by using a pinned-end contact line boundary condition which is physically relevant in specific cases, such as brim-full containers (Benjamin & Scott Reference Benjamin and Scott1979). Apart from that, various remedies were developed to remove the singularity, such as the common approach of using a slip model in which the fluid is assumed to slip in a small neighborhood near the contact line (Dussan Reference Dussan1979). The outer region where the no-slip applies, and the slip region can be perfectly matched as done in Sibley, Nold & Kalliadasis (Reference Sibley, Nold and Kalliadasis2015a). For a comprehensive comparison between approaches to tackle this singularity, please refer to Snoeijer & Andreotti (Reference Snoeijer and Andreotti2013) and Sibley et al. (Reference Sibley, Nold, Savva and Kalliadasis2015b).
In this context, using the range of parameters relevant to this work, we compared the results produced by both limits of Hoking's law (4.3), a pinned contact line vs a free-moving contact line. We found the difference in temporal growth rate to be clearly visible but consistently below 
$5\,\%$. This shows that for the range of parameters we considered, the implementation of a free-moving contact line condition and the associated singularity do not affect the qualitative accuracy of our results. Further investigation of capillary-dominated flows is planned for future work.
5. Conclusions
Using temporal stability analysis and experimental approaches, we study the linear stability of falling liquid films bounded by side walls. The two agree very well for small and large wavenumbers. Noticeable discrepancies were observed before and after the cusp in the growth rate. We showed that the supremum of the resolvent norm is large in these two regions, indicating the presence of stronger receptivity to disturbances. We plan to investigate this further in a future work.
The theoretical stability model showed that the solution of the confined 2-D problem shows several stability modes, unlike the classical 1-D problem where only the H-mode is present. In the 
$k$–
$G \sin (\beta )$ plane, we observe the fragmentation effect on the unstable region which was reported by Kögel & Aksel (Reference Kögel and Aksel2020). The unstable region splits into two smaller regions separated by a stable band at moderate wavenumbers. As confinement becomes stronger, the upper unstable region disappears. The stabilising effect is depending on the inclination angle. As inclination angle decreases, the magnitude of damping increases, and also the band of stabilised wavenumbers shifts towards higher values.
The reason behind this behaviour is the existence of different stability modes. For a certain range of parameters, the 2-D modes except a pair converge to the 1-D H-mode for weak confinement configurations. The unique pair is named W-modes. As confinement strengthens, H-modes and W-modes are stabilised in a different manner along the wavenumber. The former is stabilised strongly at low wavenumbers, while the latter experience stabilisation at high wavenumbers. An attraction centre involving the most unstable H-mode and a more stable W-mode also contributes to the formation of the stabilised wave band. The spatial structure of the eigenmodes shows an alternative even–odd symmetry in both H-modes and W-modes. The spatial structure of the most unstable H-mode shows a maximum in the centre of the channel, while that of the most unstable W-mode shows a maximum at the walls. The less-stable modes are an alternating symmetry with a further degree in the harmonics. The eigenvalues spectra shows that the 2-D eigenvalues follow the 1-D oblique problem for small wavenumbers regardless of the confinement degree, and it deviates at larger wavenumbers as confinement becomes stronger.
6. Future work
The results presented in this work can stimulate future research. The biglobal stability problem results in a multiple stability modes with comparable growth rates. This could lead to very rich nonlinear dynamics, especially that the stability modes are structural harmonics of each other (Ducimetière et al. Reference Ducimetière, Gallaire, Lefauve and Colm-Cille2021). One way to examine this can be through 3-D direct numerical simulations which was previously carried out for falling films but with periodic boundary conditions in the spanwise direction (Dietze et al. Reference Dietze, Rohlfs, Nährich and Scheid2014). Moreover, the excitation introduced by the paddle can be included in the theoretical model as a forcing term. The forced problem can also be solved and compared with the existing autonomous problem as in the work of Fabre et al. (Reference Fabre, Longobardi, Citro and Luchini2020). Finally, there exist a plethora of applications that involve volatile liquid films falling down a heated surface. The instabilities induced by thermocapillary forces (S-mode) or vapour recoil (E-mode) were only examined in the 2-D configurations in the classical Orr–Sommerfeld problem (Mohamed & Biancofiore Reference Mohamed and Biancofiore2020) or in the long-wave expansion regime (Mohamed, Dallaston & Biancofiore Reference Mohamed, Dallaston and Biancofiore2021). The analysis of the effect of spanwise confinement on this kind of instabilities is of great benefit due to the intrinsic spanwise nature of the thermal instabilities.
Acknowledgements
The authors would like to acknowledge M. Märkl for helping in carrying out the experiments. The authors would like to thank the reviewers for their insightful comments and suggestions, which helped improve the quality of this work.
Declaration of interests
The authors report no conflict of interest.
Open access
