# Stability of an elongated thickness fluctuation in a horizontal soap film

Isabelle Cantat, Corentin Tregouet

► To cite this version:

Isabelle Cantat, Corentin Tregouet. Stability of an elongated thickness fluctuation in a horizontal soap film. Journal of Fluid Mechanics, 2025, 1003, pp.A22. 10.1017/jfm.2024.1204. hal-04923106

HAL Id: hal-04923106

https://hal.science/hal-04923106v1

Submitted on 31 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stability of an elongated thickness fluctuation in a horizontal soap film

Isabelle Cantat<fnref n="1" />,† and Coretin Trégouët<fnref n="1" />,<fnref n="2" />

(Received 29 May 2024; revised 25 November 2024; accepted 10 December 2024)

Even though liquid foams are ubiquitous in everyday life and industrial processes, their ageing and eventual destruction remain a puzzling problem. Soap films are known to drain through marginal regeneration, which depends upon periodic patterns of film thickness along the rim of the film. The origin of these patterns in horizontal films (i.e. neglecting gravity) still resists theoretical modelling. In this work, we theoretically address the case of a flat horizontal film with a thickness perturbation, either positive (a bump) or negative (a groove), which is initially invariant under translation along one direction. This pattern relaxes towards a flat film by capillarity. By performing a linear stability analysis on this evolving pattern, we demonstrate that the invariance is spontaneously broken, causing the elongated thickness perturbation pattern to destabilise into a necklace of circular spots. The unstable and stable modes are derived analytically in well-defined limits, and the full evolution of the thickness profile is characterised. The original destabilisation process we identify may be relevant to explain the appearance of the marginal regeneration patterns near a meniscus and thus shed new light on soap-film drainage.

**Key words:** foams, thin films, lubrication theory

## 1. Introduction

In dry liquid foams, coarsening or imposed deformations constantly renew the contact network between bubbles: some bubbles come into contact while others separates from each other. These topological transformations are known as T1 events (Weaire & Hutzler 2000). When two bubbles come into contact, a new film is produced, with an initial thickness of the order of the meniscus size. This film thickness is generally much larger than the range of the disjoining pressure. Its subsequent drainage is therefore initially

†Email address for correspondence: isabelle.cantat@univ-rennes1.fr

<footnote n="1">¹Université de Rennes, CNRS, IPR (Institut de Physique de Rennes), UMR 6251, 35000 Rennes, France</footnote>

<footnote n="2">²MIE, CBI, ESPCI Paris, Université PSL, CNRS, 75005 Paris, France</footnote>
I. Cantat and C. Trégouët

governed solely by hydrodynamical laws (Denkov et al. 2008; Petit et al. 2015), until most of its volume has flowed in the surrounding menisci and the interfaces start to repel each other. After this hydrodynamical regime, the disjoining pressure regime begins: the drainage is governed by short-range forces, leading the film towards its equilibrium thickness (Stubenrauch & Klitzing 2003; Zhang & Sharma 2018; Chatzigiannakis & Vermant 2020), until its rejuvenation by a bubble motion, or its collapse. This cycle of thick-film production, film ageing and thin-film stabilisation has recently been shown to govern the average thickness in sheared foam (Saint-Jalmes & Trégouët 2023). It therefore influences foam coarsening, drainage and lifetime.

Since the seminal work of Mysels, Shinoda & Frankel (1959), and despite recent improvements (Chan, Klaseboer & Manica 2011; Chatzigiannakis, Jaensson & Vermant 2021), some aspects of hydrodynamic drainage have yet to be modelled and understood. A key result in the field was proposed by Barigou & Davidson (1994) and established by Aradian, Raphaël & de Gennes (2001), who predicted the asymptotic behaviour of a horizontal film of uniform thickness $h_∞$ in contact with a meniscus. Capillary suction leads to the formation of a groove along the meniscus, which becomes deeper with time and spreads into the film until the disjoining pressure limits the film thinning. This is the process at the origin of the formation of dimples in thin films. The model of Aradian et al. provides an exact solution to the hydrodynamical problem of film capillary drainage, under the important assumption that the flow satisfies the symmetry of the problem, i.e. an invariance under translation along the direction of the meniscus axis (referred to as invariant solutions hereafter). The stability of this solution is an open question and is clearly called into question by experimental observations.

Experimentally, the contact between a film and meniscus has mainly been investigated in two configurations: in vertical films, under the name of marginal regeneration (Stein 1991; Nierstrasz & Frens 1998, 1999), and in bubbles floating at a fluid interface (Lhuissier & Villermaux 2012; Frostad et al. 2016; Miguet et al. 2021). In both geometries, instabilities are observed at the bottom meniscus and the unstable thin groove is thus below the remaining part of the film. In such cases gravity has a destabilising effect, by a Rayleigh-Taylor mechanism (Shabalina et al. 2019). However, a similar instability is observed in horizontal films (Joye, Hirasaki & Miller 1994; Velev et al. 1995; Trégouët & Cantat 2021) indicating the presence of other destabilising processes.

The theoretical determination of the stability of the invariant solution is difficult: this invariant profile is indeed evolving in time at a rate controlled by the thin-film equation (Oron, Davis & Bankoff 1997; Chomaz 2001), which involves nonlinearities and high-order derivatives. As a consequence, previous analytical approaches are based on strong geometrical approximations. More specifically, Bruinsma (1995) uses the same equation set as we do, but assumes a quasi-flat invariant profile and shows that, under this assumption, the groove destabilisation requires gravity. A more refined model is numerically solved by Shi, Fuller & Shaqfeh (2022), and reproduced the destabilisation observed by Frostad et al. (2016). Here again, gravity is present, so the simulation does not allow conclusions to be drawn about the stability of the groove in a horizontal film.

In the hope of identifying the minimal ingredients responsible for the destabilisation of such an invariant groove in a soap film, we propose in this paper to address the case of an invariant groove or bump embedded in an otherwise flat film, so without meniscus. In that situation, the problem becomes analytically solvable and a destabilisation process is clearly identified: the invariant groove or bump is destabilised by capillary effects only. The analytical solution is made possible because of the self-similar properties of the invariant thickness profile established by Benzaquen, Salez & Raphaël (2013). As this reference
Stability of elongated thickness fluctuation in a soap film

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic representation of the film and of its non-dimensional dimensions. The reference bump is modulated with a wavelength $\lambda$. A film element is shown in blue, with the expression of the force exerted by the external film on its purple face (f). The black arrows on the left represent the parabolic velocity field, with the interfacial velocity $v$ in red.</figcaption></figure>

solution is unsteady, its destabilisation is not exponential in time and we show that the amplitude of a small perturbation breaking the invariance of the reference solution grows as $t^{7/4}$ at large time and large wavelength. The destabilisation mechanism we identify should still be present in the case of a groove between a thin film and a meniscus, so we believe that this study sheds new light on the old and important problem of thin-film drainage.

Section 2 is dedicated to the physical assumptions underlying the model and to the resulting thin-film equations, used in the limit of large Gibbs elasticity; § 3 discusses the properties of the reference thickness profile, invariant along one direction. It is the self-similar solution of the linearised thin-film equation whose analytical expression was established by Benzaquen et al. (2013). This solution is then perturbed by a sinusoidal variation of its height and width of wavelength $\lambda$, and a linear stability analysis is performed in § 4, leading to a linear, time-dependent, set of equations giving the growth rate of the perturbation. This equation set is eventually solved in § 5, where exact polynomial expressions for a stable and an unstable mode are determined.

## 2. General assumptions and equations of motion in a soap film

### 2.1. *Physical properties of the film*

We consider a film of thickness $2h_\infty$ with a bump localised along the y axis. The reference bump is invariant along the y direction and symmetric with respect to the planes (0, x, y) and (0, y, z). The aim of our work is to establish that an arbitrarily small modulation of the bump height in the y direction (satisfying the two mirror symmetries; see figure 1) is amplified with time, meaning that the invariance under translation of the reference bump is spontaneously broken. Note that the film tends to minimise its interface area and thus relaxes to a flat final state. The question is therefore to determine if its transient shape, during the relaxation process, remains invariant under translation or not.

The surface tension is $\gamma_0$ at large $|x|$, far from the bump, and slightly varies close to the bump. In order to establish the fact that capillary effects alone are responsible for
I. Cantat and C. Trégouët

the spontaneous symmetry breaking, we neglect all other potential destabilising physical processes: disjoining pressure, gravity and inertial effects are not taken into account in our model. For the sake of simplicity, we also neglect the air friction: the damping forces are due to the bulk viscosity η and to the interface viscosity η_s, with the additional assumption that η_s ≫ h_∞η. We assume that the interface Gibbs elasticity is large enough so that the interfaces can be considered as incompressible (Mysels et al. 1959), which is the relevant limit for soap films with small (Seiwert, Dollet & Cantat 2014; Champougny et al. 2015) or vanishing (Lenavetier et al. 2024) external forcing. This condition imposes that

$$ \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} = \mathrm{div}^{2D}(\mathbf{v}) = 0, \qquad (2.1) $$

with $\mathbf{v}(x, y, t)$ the interfacial velocity, assumed to be identical on both interfaces of the film. In this limit the surface tension plays the role of the Lagrange multiplier that prevents the surface from compressing/dilating and the hydrodynamical problem can thus be solved without an additional equation of state for the interface.

The problem is made non-dimensional using $h_∞$ as the length unit, $τ = 3ηh_∞/γ_0$ as the time unit and $γ_0$ as the tension unit. In the following, we use the bar notation for dimensional variables. The equations of motion are determined from the mass and force balances on a film element, defined as the liquid trapped between two pieces of interface dx dy, chosen symmetrically in the top and bottom interfaces and followed along their trajectories (see figure 1). The bump width w in the x direction is assumed to be much larger than the film thickness. Under this assumption, the lubrication equations can be used, which are Taylor expansions of the Stokes equations, with the thickness gradient as a small parameter.

2.2. Poiseuille flow

The evolution of the film thickness $2h(x, y, t)$ is determined from the mass conservation of a film element (Bruinsma 1995):

$$ \partial_t h + \mathbf{v} \cdot \nabla h = -\nabla \cdot (h^3 \nabla \Delta h), \qquad (2.2) $$

with $\nabla$ and $\Delta$ the gradient and Laplacian two-dimensional operators in the (x, y) plane. The left-hand side is the Lagrangian time derivative of the half-thickness of a film element followed along its trajectory. The right-hand side quantifies the flux of liquid leaving the film element due to the relative motion of the liquid with respect to the interfaces (see figure 1). This motion is a Poiseuille flow driven by the Laplace pressure $\Delta h$, with a mean velocity $h^2\nabla\Delta h$ and a flux $h^3\nabla\Delta h$. Note that, without this Poiseuille flow, the film element becomes a closed system of constant volume $2h \, dx \, dy$. In that case, the thickness is conserved, as the condition $\mathrm{div}^{2D}(\mathbf{v}) = 0$ imposes the conservation of the film element area $dx \, dy$.

2.3. Capillary and viscous stress tensors

The second relationship between $h$ and $\mathbf{v}$ is obtained from a force balance on a film element, which also involves the non-dimensional surface tension $\gamma = 1 + \delta\gamma$. An example of a film element of volume $2h \, dx \, dy$ is shown in blue in figure 1. The surface tension and pressure forces acting on one of its lateral faces, e.g. its purple face, can be expressed as a function of a two-dimensional stress tensor $\sigma_{cap}$. This tensor is defined so that the force $dF^f$ acting on a face (f) of area $2h \, d\ell^f$ and of normal $\mathbf{n}^f$ is
$d\boldsymbol{F}^f = \boldsymbol{\sigma}_{\text{cap}} \cdot \boldsymbol{n}^f d\ell^f$. In the case of the purple face shown in figure 1, the area is $2h$ dy and the normal $\boldsymbol{n}^f = \boldsymbol{e}_x$.

This stress depends on the local tension value, on the Laplace pressure $\Delta h$ integrated over the film thickness $2h$ and on the thickness gradient $\nabla h$. Indeed the interface slope leads to quadratic corrections on the capillary force. The expression of the stress tensor, established in Lenavetier et al. (2024), is written below, both in the reference basis $\mathcal{B}_0 = (\boldsymbol{e}_x, \boldsymbol{e}_y)$, with $\boldsymbol{e}_x$ and $\boldsymbol{e}_y$ the unit vectors associated with the directions x and y shown in figure 1, and in a local basis $\mathcal{B}_e = (\boldsymbol{n}, \boldsymbol{t})$ associated with the local film geometry. The unit vector $\boldsymbol{n}$ is oriented in the direction of the thickness gradient, so that $\nabla h = |\nabla h|\boldsymbol{n}$ and $\boldsymbol{t}$ is a unit vector oriented along the lines of constant thickness, $\boldsymbol{t} = \boldsymbol{e}_z \wedge \boldsymbol{n}$. With these conventions, we get

$$ \boldsymbol{\sigma}_{\text{cap}} = \boldsymbol{\sigma}_{\text{cap}}^* + \sigma^f \boldsymbol{I}, \qquad (2.3) $$

with $\boldsymbol{I}$ the identity matrix. The isotropic term

$$ \sigma^f = 2(1 + \delta\gamma) + 2h\Delta h \qquad (2.4) $$

comes from the Laplace pressure and tension forces. The deviatoric part is

$$ \boldsymbol{\sigma}_{\text{cap}}^* = \begin{pmatrix} -|\nabla h|^2 & 0 \\ 0 & |\nabla h|^2 \end{pmatrix}_{\mathcal{B}_e} = - \begin{pmatrix} (\partial_x h)^2 - (\partial_y h)^2 & 2\partial_x h \partial_y h \\ 2\partial_x h \partial_y h & (\partial_y h)^2 - (\partial_x h)^2 \end{pmatrix}_{\mathcal{B}_0} \qquad (2.5) $$

It arises from the correction of the tension forces, due to their projection in the $(x, y)$ plane. A film element is thus submitted to an anisotropic stress, which is higher in the direction perpendicular to the thickness gradient.

The resulting force on the film element is

$$ F_{\text{cap}} = \nabla \cdot \boldsymbol{\sigma}_{\text{cap}} = 2h\nabla(\Delta h) + 2\nabla\delta\gamma. \qquad (2.6) $$

For a film of arbitrary thickness, the first term $2h\nabla(\Delta h)$ is a priori not the gradient of a function and thus cannot be balanced by the second term $2\nabla\delta\gamma$, whatever the value of the surface tension. The force balance thus requires an additional term, which is here assumed to be the interfacial viscous friction. The viscous stress in a single interface is given by

$$ \boldsymbol{\sigma}_{\text{vis}} = \zeta(\nabla\boldsymbol{v} + (\nabla\boldsymbol{v})^T), \qquad (2.7) $$

with $\zeta = \eta_s/(\tau\gamma_0) = \eta_s/(3\eta h_\infty)$ a dimensionless friction coefficient.

With the interfacial incompressibility condition, the resulting viscous force in a single interface is $F_{\text{vis}} = \zeta \Delta \boldsymbol{v}$ and the force balance becomes (Bruinsma 1995)

$$ \zeta \Delta \boldsymbol{v} + \nabla \delta \gamma + h \nabla (\Delta h) = 0. \qquad (2.8) $$

Equations (2.1), (2.2), (2.8) constitute a closed set of equations for the unknown functions $h$, $\delta\gamma$ and $\boldsymbol{v}$ and determine the film evolution from any initial state. This equation set corresponds to the limit of high Marangoni number and low Bond number of the equation set numerically solved in Shi et al. (2022), with $\zeta$ based on the interfacial viscosity instead of the bulk viscosity.

Our initial condition is shown in figure 1: a flat film with a bump whose height and width are modulated with a wavelength $\lambda$ in the y direction. We restrict our stability analysis to the case $\lambda \gg w$. This assumption allows us to define a near-field domain where $x$ is of the order of $w$ and a far-field domain where $x$ is of the order of $\lambda$. The thickness only varies at the small scale, whereas the velocity is shown to vary at the large scale. We take advantage of this scale separation to simplify further (2.8).
I. Cantat and C. Trégouët

<figure><img src="image_4.png" /></figure>

Figure 2. Scheme of the capillary and viscous forces in the vicinity of the bump. The light-blue domain represents the bump and the grey rectangle the film element of interest. The stress on the red edge, spanning from one side of the bump to the other, is dominated by the capillary stress. Its integral along the edge is, by definition, the line tension. On the black edges, the capillary stress tensor is zero, and the viscous forces dominate.

2.4. Line tension

As the thickness gradients are localised in the near-field domain, the line tension formalism is used to express the capillary forces at the large length scale, as proposed in Lenavetier et al. (2024). In our geometry, the line tension is defined, for a single interface, as

$$
T = \frac{1}{2} \int_{-\infty}^{\infty} \mathbf{e}_y \cdot (\boldsymbol{\sigma}_{\text{cap}} - 2\mathbf{l}) \cdot \mathbf{e}_y \, dx. \qquad (2.9)
$$

It represents the force excess in the e_y direction exerted on a line oriented in the x direction, crossing the bump (see figure 2). In the limit where ∂_x h ≫ ∂_y h, its expression, established in Lenavetier et al. (2024), is

$$
T = \int_{-\infty}^{\infty} \left( \frac{\partial h}{\partial x} \right)^2 dx. \qquad (2.10)
$$

2.5. Far-field expression of the force balance

We anticipate that the thickness modulation of wavelength $\lambda$ induces a motion in a domain
of width $\lambda$ on both sides of the bump. In this far-field domain, the thickness gradients are
negligible and (2.8) becomes, for $x \sim \lambda \gg w$,

$$
\zeta \Delta \mathbf{v} + \nabla \delta \gamma = 0. \tag{2.11}
$$

The motion is then driven by the boundary conditions at $x \sim w \ll \lambda$. The force balance on the piece of film shown in figure 2, of characteristic size $w$ dy, leads to the condition

$$
2 \frac{\partial T}{\partial y} = -4 \mathbf{e}_y \cdot \boldsymbol{\sigma}_{\text{vis}} \cdot \mathbf{e}_x (x = w). \qquad (2.12)
$$

As we restrict our stability analysis to degrees of freedom that keep the mirror symmetry
$x \to -x$, the $x$ component of the velocity tends to 0 at small scale and we thus have
$\partial v_x / \partial y \ll \partial v_y / \partial x$ at $x \sim w$. Using (2.7), the condition (2.12) becomes

$$
\frac{\partial T}{\partial y} = -2\zeta \frac{\partial v_y}{\partial x} (x = 0), \qquad (2.13)
$$

where we used $\partial v_y / \partial x (x = 0) = \partial v_y / \partial x (x = w)$ at first order in $w / \lambda$.

To solve the whole stability problem, the far-field velocity field is determined from (2.11) solved for $x > 0$, with the boundary condition (2.13). The behaviour of this far-field solution at small $x$ is used to determine the deformation of the thickness pattern due to convection in (2.2).
### 3. Choice and properties of the reference bump

The investigation of the bump stability first requires characterisation of the evolution of a bump invariant under translation, which is our reference state.

#### 3.1. Self-similar solution

We first impose that the dynamics of this reference bump remains invariant under translation along the y direction. In this case, the condition of incompressible interfaces $div^{2D}v = 0$ becomes $\partial v_x/\partial x = 0$. Associated with the condition $v_x(0) = 0$ given by the symmetry, it imposes the condition of vanishing velocity at the interface. The set of equations of motion is strongly simplified: (2.2) becomes

$$ \partial_t h = -\partial_x(h^3 \partial_{xxx} h). \qquad (3.1) $$

This is a very classical equation which also governs thin films deposited on a solid (Landau & Levich 1942; Hammond 1983; Oron et al. 1997; Chomaz 2001; Lister et al. 2006).

Equation (2.8) is no longer useful, but consistently admits a solution with zero velocity for any y-invariant thickness field. Indeed it leads to

$$ \delta\gamma = -\int_{-\infty}^{x} h \partial_{xxx} h \, dx', \qquad (3.2) $$

which provides the tension value if needed.

Equation (3.1) is analytically solved in Benzaquen et al. (2013) in the linear regime. For a small bump, it becomes

$$ \partial_t h = -\partial_{xxxx} h. \qquad (3.3) $$

The authors show that the profile converges at long time towards an asymptotic self-similar profile:

$$ h^{\text{as}}(x, t) = 1 + \delta h_0(t) \Phi \left( \frac{x}{w_0(t)} \right), \qquad (3.4) $$

where $\Phi$ is the function plotted in figure 3. It is an even function, decreasing exponentially fast at large $x$ and normalised so that $\int_{-\infty}^{\infty} \Phi(U) dU = 1$. The functions $\delta h_0$ and $w_0$ are

$$ \delta h_0(t) = At^{-1/4} \quad \text{and} \quad w_0(t) = t^{1/4}, \qquad (3.5a,b) $$

with $A$ the area of the bump section, which is conserved by conservation of the liquid volume. The height of the bump is $\Phi(0)\delta h_0$ and decreases with time, whereas the bump width $w_0$ increases: the bump spreads on the soap film until it disappears.

The explicit expression of $\Phi$ is given in Benzaquen et al. (2013) and will not be useful in the following. However, we use a few relationships involving this function and its derivative $\Phi'$, and more generally its derivative $\Phi^{(n)}$ of order $n$, which are established below. To shorten the notations, we introduce the auxiliary functions $\Phi_0(x, t) = \Phi(x/w_0(t))$ and $\Phi'_0(x, t) = \Phi'(x/w_0(t))$.
<figure><img src="image_5.png" /><figcaption>Figure 3. Graph of the function $\Phi$ given by equation (19) of Benzaquen et al. (2013) (solid line). The function $-U\Phi'(U)$ (dashed line) is also used in § 4.</figcaption></figure>

Using the expressions (3.4) and (3.5a,b), we get the time and space derivatives of $h^{as}$:

$$ \partial_t h^{as} = -\frac{A}{4}t^{-5/4}(\Phi_0 + xt^{-1/4}\Phi'_0) = -\frac{1}{4}\frac{\delta h_0}{w_0^4}\left(\Phi_0 + \frac{x}{w_0}\Phi'_0\right), \quad (3.6) $$

$$ \partial_{xxxx} h^{as} = \frac{\delta h_0}{w_0^4} \Phi_0^{(4)}. \quad (3.7) $$

Using the expressions (3.6) and (3.7) in (3.3), we obtain the differential equation

$$ 4\Phi^{(4)}(U) = \Phi(U) + U\Phi^{(1)}(U). \quad (3.8) $$

These properties are used below in the following form:

$$ \partial_t(\delta h_0 \Phi_0) = -\frac{\delta h_0}{w_0^4} \Phi_0^{(4)} = -\frac{\delta h_0}{4w_0^4} \left( \Phi_0 + \frac{x}{w_0} \Phi_0^{(1)} \right), \quad (3.9) $$

directly deduced from (3.3), (3.7) and

$$ \partial_t \left( \frac{\delta h_0}{w_0} \Phi'_0 \right) = \partial_{t,x} (\delta h_0 \Phi_0) = -\partial_x \left( \frac{\delta h_0}{w_0^4} \Phi_0^{(4)} \right) = -\frac{\delta h_0}{w_0^5} \Phi_0^{(5)}. \quad (3.10) $$

In this equation, the first and last equalities are obtained by spatial integration and differentiation, respectively, and the second equality is obtained by switching the time and space derivatives and using (3.9).

## 3.2. Line tension of the reference solution

This bump generates a line tension given by (2.10), the expression for which becomes

$$ T_0 = \int_{-\infty}^{\infty} \left( \frac{\delta h_0}{w_0} \Phi'(U) \right)^2 w_0 \, dU = I_1 \frac{\delta h_0^2}{w_0}, \quad (3.11) $$

with

$$ I_1 = \int_{-\infty}^{\infty} [\Phi'(U)]^2 dU \approx 0.058. \quad (3.12) $$
Stability of elongated thickness fluctuation in a soap film

<figure><img src="image_6.png" /><figcaption>Figure 4. Scheme of the thickness field given by (4.1). The crest line of the bump is at height $1 + \delta h(y, t)\Phi(0)$, and the characteristic width of the bump is $w(y, t)$.</figcaption></figure>

## 4. Governing equations for a small non-invariant fluctuation

In §§ 4 and 5, we determine an explicit unstable mode of wavelength $\lambda$ for this y-invariant self-similar solution, and show that this unstable mode grows faster than the reference bump decreases. It is built using two simple degrees of freedom, the width and the height of the self-similar solution.

### 4.1. Initial perturbation

We assume that, at a given time $t_0$, an arbitrarily small fluctuation of wavelength $\lambda = 2\pi/k$ is imposed to the y-invariant self-similar profile $h^{as}$ of section area A. Note that the choice of the initial time $t_0$ imposes the initial width and height of the reference bump through the relations (3.5a,b). The solution thus depends on the choice of this initial time. For the evolution of the modulated bump, at $t > t_0$, the function is sought in the form

$$h(x, y, t) = 1 + \delta h(y, t) \Phi \left( \frac{x}{w(y, t)} \right), \quad (4.1)$$

with

$$\delta h = \delta h_0(t)(1 + \varepsilon_h(y, t)) \quad \text{and} \quad w = w_0(t)(1 + \varepsilon_w(y, t)), \quad (4.2a,b)$$

as shown in figure 4. This form presumes that, at each y and t, the thickness profile keeps its self-similar shape $\Phi$, which is a strong assumption. The existence of a solution satisfying the form of (4.1) is thus not ensured at this stage, but will be proved a posteriori by exhibiting a solution.

The functions $\delta h_0(t)$ and $w_0(t)$ are given by (3.5a,b), and the small parameters controlling the fluctuation amplitude are $\varepsilon_h = \hat{\varepsilon}_h(t)e^{iky}$ for the height and $\varepsilon_w = \hat{\varepsilon}_w(t)e^{iky}$ for the width of the bump profile.

The functions $\hat{\varepsilon}_h$ and $\hat{\varepsilon}_w$ are assumed to depend on time only. We assume that, for $t > t_0$, $\delta h_0 \ll 1$, $w_0 \ll 1/k$, $\varepsilon_h \ll 1$ and $\varepsilon_w \ll 1$. The equation of motion (2.2) is linearised with respect to $\delta h_0$, $\varepsilon_h$ and $\varepsilon_w$. Equations (2.11) and (2.13) are of quadratic order in $\delta h_0$ and are linearised in $\varepsilon_h$ and $\varepsilon_w$ only.

### 4.2. Qualitative analysis

Before tackling the full analytical problem, we first qualitatively identify the physical processes involved.
<figure><img src="image_7.png" /><figcaption>Figure 5. Schematic top view of the film, with the colour indicating the domain of various thicknesses: the flat film (light blue), the bump (blue) and a positive height fluctuation (darker blue). (a) The local width $w^{loc}$ is larger than the average width $w^{av}$. (b) The local height $h^{loc}$ is larger than the average height $h^{av}$. The amplitude and orientation of the Poiseuille flow, relative to the interfaces, are given by the blue arrows, and the velocity of the interfaces is represented by the red arrows.</figcaption></figure>

The figure 5 provides a schematic basis to discuss the influence of the height and width perturbations on the time evolution of the bump. The Poiseuille flow, localised at $|x| \sim w$ and represented with blue arrows, is the motion of the fluid with respect to the interfaces. It occurs perpendicularly to the bump and is oriented towards the film. It tends to spread the bump by increasing its width and decreasing its height.

Let us first assume that the local width $w^{loc}$ is larger than the average width $w^{av}$, the local bump height $\delta h^{loc}$ remaining at the average height $\delta h^{av}$, as shown in figure 5(a). Locally, the interface curvature and thus the Laplace pressure below the bump are smaller than their average values and the Poiseuille flow (P) decreases. As a consequence, the width increases slower than in the rest of the bump ($0 < dw^{loc}/dt|_P < dw^{av}/dt|_P$) and the height decreases slower ($\delta h^{av}/dt|_P < \delta h^{loc}/dt|_P < 0$). The Poiseuille flow is thus stabilising the width: as $w^{loc}$ grows slower than the width of the remaining part of the bump, the width is uniform again after some time.

However, as the flow no longer needs to remain invariant under translation along the y direction, the interfacial velocity $v$, which was zero in § 3, now plays an important role: the flow is the superposition of P discussed above and a plug flow (i.e. a flow that is uniform through the entire thickness of the film) at a velocity equal to the interfacial velocity $v$. The driving force for this second part of the flow results from the gradients of line tension (T). The line tension scales as $\delta h/w^2$ (see (3.11)) and a locally larger width thus induces a smaller line tension. The resulting line-tension gradient along the bump axis y drives a circulation in the film, schematised by the red arrows and circles in figure 5(a): at the point of low line tension, the bump is stretched along its axis y and, by interface incompressibility, it is compressed in the transverse direction x. This part of the flow leaves the thickness modulation unchanged ($\delta h^{av}/dt|_T = \delta h^{loc}/dt|_T = 0$) but decreases the width ($dw^{loc}/dt|_T < dw^{av}/dt|_T = 0$): this is again a stabilising ingredient for the width.

The total influence of the width fluctuation is obtained by adding the effects of the pressure-induced Poiseuille flow (P) and that of the tension-induced recirculation flow (T). It can be written in the linear regime as follows:

$$ \frac{d(\delta h^{loc} - \delta h^{av})}{dt} = a_{hw}(w^{loc} - w^{av}) \quad \text{and} \quad \frac{d(w^{loc} - w^{av})}{dt} = -a_{ww}(w^{loc} - w^{av}), \qquad (4.3a,b) $$

with $a_{wh}$ and $a_{ww}$ positive numbers.
The same reasoning can be done with an initial fluctuation of the bump height, as represented in figure 5(b). At a point of larger thickness ($\delta h^{loc} > \delta h^{av}$), the line tension and Laplace pressure are larger, so the opposite scenario occurs: the Poiseuille contribution leads to a faster height decrease $d\delta h^{loc}/dt|_P < d\delta h^{av}/dt|_P < 0$ and width increase $0 < dw^{lav}/dt|_P < dw^{loc}/dt|_P$, whereas the tension contribution does not modify the thickness but imposes $0 = dw^{av}/dt|_T < dw^{loc}/dt|_T$. The sum of the two contributions leads to

$$ \frac{d(\delta h^{loc} - \delta h^{av})}{dt} = -a_{hh}(\delta h^{loc} - \delta h^{av}) \quad \text{and} \quad \frac{d(w^{loc} - w^{av})}{dt} = a_{wh}(\delta h^{loc} - \delta h^{av}), \qquad (4.4a,b) $$

with $a_{hh}$ and $a_{hw}$ positive numbers.

Finally, considering a single degree of freedom (either the height only or the width only) leads to the conclusion that the invariant solution is stable. However, considering together these two degrees of freedom allows us to built a coupled system which can, in contrast, lead to an instability. Indeed the eigenvalues of the system (4.3a,b), (4.4a,b) satisfy

$$ \omega^2 + \omega(a_{hh} + a_{ww}) + a_{hh}a_{ww} - a_{hw}a_{wh} = 0. \qquad (4.5) $$

If $a_{hh}a_{ww} < a_{hw}a_{wh}$, one solution is positive and the perturbation grows.

In the following, we determine the coupling coefficients associated with the growth of the unknown functions $\varepsilon_h$ and $\varepsilon_w$ and we show that, indeed, the invariant solution is unstable. The whole destabilisation process can be summarised by the following positive feedback: a positive height fluctuation leads to a higher tension, thus to a local axial compression, thus to a bump width increase, thus to a decrease of the Poiseuille flow, thus to a slowdown of the height decrease and finally to an increase of the initial fluctuation.

### 4.3. First linearisations

The different equations established in § 2 need to be expanded to first order in $\varepsilon_h$ and $\varepsilon_w$. For any function $F(x/w)$ implicitly written $F$, we define the function $F_0 = F(x/w_0)$. In the following, we use that, at first order in $\varepsilon_w$, we have

$$ F = F_0 - \frac{x}{w_0} \varepsilon_w F'_0. \qquad (4.6) $$

#### 4.3.1. Thickness profile derivatives

The thickness profile given by (4.1) is, at first order,

$$ h = 1 + \delta h_0 (1 + \varepsilon_h) \Phi \left( \frac{x}{w_0 (1 + \varepsilon_w)} \right) = 1 + \delta h_0 \Phi_0 + \varepsilon_h \delta h_0 \Phi_0 - x \varepsilon_w \frac{\delta h_0}{w_0} \Phi'_0, \qquad (4.7) $$

and its time derivative is

$$ \partial_t h = \partial_t (\delta h_0 \Phi_0) + \varepsilon_h \partial_t (\delta h_0 \Phi_0) - x \varepsilon_w \partial_t \left( \frac{\delta h_0}{w_0} \Phi'_0 \right) + \dot{\varepsilon}_h \delta h_0 \Phi_0 - x \dot{\varepsilon}_w \frac{\delta h_0}{w_0} \Phi'_0. \qquad (4.8) $$
### 4.3.2. Poiseuille flow

We first linearise the Poiseuille term of (2.2), leading to

$$-\nabla \cdot (h^3 \nabla \Delta h) = -\nabla \cdot (\nabla \Delta h). \quad (4.9)$$

As $1/k \ll w$, the thickness variations are dominated by the $x$ derivative, so

$$\Delta h = \frac{\delta h}{w^2} \Phi'' \quad \text{and} \quad -\nabla \cdot \nabla \Delta h = -\frac{\delta h}{w^4} \Phi^{(4)}. \quad (4.10a,b)$$

This latter quantity must now be linearised with respect to $\epsilon_w$ and $\epsilon_h$, and using (3.9), (3.10) and (4.6), we get

$$\left.
\begin{aligned}
-\nabla \cdot \nabla \Delta h &= -\frac{\delta h_0}{w_0^4} \Phi^{(4)}(1 + \epsilon_h - 4\epsilon_w), \\
-\nabla \cdot \nabla \Delta h &= -\frac{\delta h_0}{w_0^4} \Phi_0^{(4)}(1 + \epsilon_h - 4\epsilon_w) - x \frac{\delta h_0}{w_0^5} \Phi_0^{(5)} \epsilon_w, \\
-\nabla \cdot \nabla \Delta h &= \partial_t (\delta h_0 \Phi_0)(1 + \epsilon_h - 4\epsilon_w) - x \partial_t \left(\frac{\delta h_0}{w_0} \Phi_0'\right) \epsilon_w.
\end{aligned}
\right\} \quad (4.11)$$

### 4.3.3. Tension

From (2.10), the tension is

$$T = \int_{-\infty}^{\infty} \left( \frac{\delta h}{w} \Phi'(U) \right)^2 w \, dU = \frac{\delta h_0^2}{w_0} (1 + 2\epsilon_h - \epsilon_w) \int_{-\infty}^{\infty} [\Phi'(U)]^2 \, dU, \quad (4.12)$$

and using (3.11), we get

$$T = (1 + 2\epsilon_h - \epsilon_w)T_0 \quad \text{and} \quad \partial_y T = ik(2\epsilon_h - \epsilon_w)T_0. \quad (4.13a,b)$$

## 4.4. Far-field recirculation

As discussed in § 2.5, the recirculation flow at the length scale $\lambda$ is governed by (2.11) and the bump properties only appear through the boundary condition (2.13). The velocity field is thus determined in the half-plane $x > 0$ from these equations.

As $\operatorname{div}^{2D} v = 0$, we introduce the stream function $\phi$ defined as

$$v_x = \partial_y \phi \quad \text{and} \quad v_y = -\partial_x \phi, \quad (4.14a,b)$$

with $\phi = 0$ at large $x$ and

$$\phi = f(x, t) e^{iky}. \quad (4.15)$$

The two components of (2.11) become

$$\left.
\begin{aligned}
\zeta(\partial_{xx} + \partial_{yy})\partial_y\phi + \partial_x\delta\gamma &= 0, \\
\zeta(\partial_{xx} + \partial_{yy})\partial_x\phi - \partial_y\delta\gamma &= 0.
\end{aligned}
\right\} \quad (4.16)$$

The unknown tension variation $\delta\gamma$ can be removed by taking the cross-derivatives and combining both equations to obtain

$$(\partial_{xx} + \partial_{yy})(\partial_{xx} + \partial_{yy})\phi = 0. \quad (4.17)$$

Using (4.15), we get a fourth-order equation of the function $f$:

$$f^{(4)} - 2k^2 f^{(2)} + k^4 f = 0, \quad (4.18)$$

which has $k$ and $-k$ as double roots.
<figure><img src="image_8.png" /><figcaption>Stability of elongated thickness fluctuation in a soap film</figcaption></figure>

Figure 6. Velocity field ($v_x$, $v_y$) given by (4.24), (4.25) as a function of $x/\lambda$ and $y/\lambda$, for $2\hat{\varepsilon}_h - \hat{\varepsilon}_w > 0$. The bump is schematically represented by the red rectangle.

In the domain $x > 0$, only the negative value leads to a non-divergent solution (the other half-plane being obtained by symmetry). The condition $v_x(x = 0, y) = 0$ imposes $f(x = 0, t) = 0$, so $f = -i\phi_1 x e^{-kx}$, with $\phi_1$ a constant determined below. It leads to

$$v_x = k\phi_1 x e^{-kx} e^{iky}, \quad (4.19)$$

$$v_y = -i\phi_1(kx - 1) e^{-kx} e^{iky}, \quad (4.20)$$

$$\partial_x v_y = -i\phi_1 [-k(kx - 1) + k] e^{-kx} e^{iky}, \quad (4.21)$$

$$\partial_x v_y(x = 0, y) = -2ik\phi_1 e^{iky}. \quad (4.22)$$

The boundary condition (2.13) imposes, using (4.13b),

$$\phi_1 = (2\hat{\varepsilon}_h - \hat{\varepsilon}_w) \frac{T_0}{4\zeta}. \quad (4.23)$$

So, far from the bump, from (4.19), (4.20) and (4.23), we get

$$v_x = (2\hat{\varepsilon}_h - \hat{\varepsilon}_w) \frac{T_0}{4\zeta} kx e^{-kx} e^{iky}, \quad (4.24)$$

$$v_y = -i(2\hat{\varepsilon}_h - \hat{\varepsilon}_w) \frac{T_0}{4\zeta} (kx - 1) e^{-kx} e^{iky}. \quad (4.25)$$

This velocity field is shown in figure 6.

The convective term involved in (2.2) is $-v \cdot \nabla h$. The term $-v_y \partial_y h$ is of second order in $\varepsilon$ so the convective term becomes $-v_x \partial_x h$, with $\partial_x h = (\delta h_0 / w_0) \Phi'_0$ at the zeroth order in $\varepsilon$. Extrapolating the far-field velocity at the bump, in the small-$x$ domain where $x e^{-kx} \sim x$, we obtain the convective term involved in the thickness profile evolution:

$$-v_x \partial_x h = -x(2\varepsilon_h - \varepsilon_w) \frac{kT_0}{4\zeta} \frac{\delta h_0}{w_0} \Phi'_0. \quad (4.26)$$
# I. Cantat and C. Trégouët

## 4.5. Equation of stability

The Taylor expansion of the equation of motion (2.2) can now be performed using (4.8), (4.11) and (4.26). At the zeroth order in $\epsilon_h$ and $\epsilon_w$ we consistently get $\partial_t(\delta h_0 \Phi_0) = \partial_t(\delta h_0 \Phi_0)$.

At first order in $\epsilon$, we get

$$ \begin{aligned} & \epsilon_h \partial_t(\delta h_0 \Phi_0) - x \epsilon_w \partial_t \left( \frac{\delta h_0}{w_0} \Phi'_0 \right) + \dot{\epsilon}_h \delta h_0 \Phi_0 - x \dot{\epsilon}_w \frac{\delta h_0}{w_0} \Phi'_0 \\ & = \partial_t(\delta h_0 \Phi_0)(\epsilon_h - 4\epsilon_w) - x \partial_t \left( \frac{\delta h_0}{w_0} \Phi'_0 \right) \epsilon_w - x \frac{\delta h_0}{w_0} \Phi'_0 (2\epsilon_h - \epsilon_w) \frac{kT_0}{4\zeta}, \end{aligned} \quad (4.27) $$

which is reorganised using (3.9) into

$$ \Phi_0 \delta h_0 \left( \dot{\epsilon}_h - \epsilon_w \frac{1}{w_0^4} \right) + x \Phi'_0 \frac{\delta h_0}{w_0} \left( -\dot{\epsilon}_w - \epsilon_w \frac{1}{w_0^4} + (2\epsilon_h - \epsilon_w) \frac{kT_0}{4\zeta} \right) = 0. \quad (4.28) $$

We observe only two kinds of x-dependencies in the different terms: either functions proportional to $\Phi_0 = \Phi(x/w_0(t))$ or functions proportional to $x\Phi'_0 = x\Phi'(x/w_0(t))$. As these two functions are not proportional (see figure 3), the equality (4.28) requires that the prefactor of each one is zero.

We thus get the conditions

$$ \left. \begin{aligned} \dot{\epsilon}_h &= \epsilon_w \frac{1}{w_0^4}, \\ \dot{\epsilon}_w &= -\epsilon_w \frac{1}{w_0^4} + (2\epsilon_h - \epsilon_w) \frac{kT_0}{4\zeta}. \end{aligned} \right\} \quad (4.29) $$

Note that this set of two linear equations no longer involves any dependency in $x$ and admits solutions that are discussed below. This justifies a *posteriori* the functional choice made for $\delta h$ in (4.1) and (4.2a,b), and the assumption that $\epsilon_h$ and $\epsilon_w$ do not depend on $x$.

## 5. Stability analysis

The linear system (4.29) depends on time through $w_0$ and $T_0$ and the stability analysis requires one to solve this time-dependent equation set. However, the instantaneous growth rate discussed below provides a first intuition of the physical process.

### 5.1. Eigenvalues of the system

In the equation set (4.29), the Poiseuille flow and the flow induced by the line tension are governed, respectively, by the characteristic rates $\omega_P = 1/w_0^4$ and $\omega_T = kT_0/(4\zeta) = (3I_1/4)(k\delta h_0^2/w_0)(\eta h_\infty/\eta_s)$. Note that the time scale is $\tau = 3\eta h_\infty/\gamma_0$, so the associated physical quantities satisfy $\tilde{\omega}_P \sim (\delta h_0/w_0)^4 \gamma_0/(\eta h_\infty)$ and $\tilde{\omega}_T \sim \delta h_0^2/(\lambda w_0) \gamma_0/\eta_s$. As anticipated in § 4.2, the driving force is always the minimisation of the interface energy, proportional to $\gamma_0$. Moreover, the damping force is either the bulk friction proportional to $\eta$, due to a motion of the liquid phase with respect to the interfaces, or the interface friction proportional to $\eta_s$, due to in-plane recirculations in the film.
<figure><img src="image_9.png" /><figcaption>Stability of elongated thickness fluctuation in a soap film</figcaption></figure>

Figure 7. Growth rate $\omega^+$ of the instability as a function of $\omega_T$ and $\omega_P$, as given by (5.2). Some level curves are shown in black, for integer values of $\omega^+$.

With these notations, the system (4.29) becomes

$$ \frac{d}{dt} \begin{bmatrix} \varepsilon_h \\ \varepsilon_w \end{bmatrix} = \begin{bmatrix} 0 & \omega_P \\ 2\omega_T & -\omega_T - \omega_P \end{bmatrix} \begin{bmatrix} \varepsilon_h \\ \varepsilon_w \end{bmatrix}. \quad (5.1) $$

Coming back to our qualitative analysis in §4.2, we get the expected signs in the matrix. Note, however, that one coupling coefficient is zero. This comes from the fact that we work here with the relative height perturbation $\varepsilon_h$ and not with the local bump height $\delta h_0(1 + \varepsilon_h)$ as previously. The Poiseuille flow is proportional to the bump height (see (4.11)), leading to the cancellation of the two terms proportional to $\varepsilon_h$ in (4.27).

The eigenvalues of (5.1) are

$$ \omega^{\pm} = \frac{1}{2} \left( -(\omega_T + \omega_P) \pm \sqrt{(\omega_T + \omega_P)^2 + 8\omega_T\omega_P} \right), \quad (5.2) $$

and the associated eigenvectors $(\varepsilon_h, \varepsilon_w)^{\pm} = (1, \omega^{\pm}/\omega_P)$. The growth rate $\omega^+$ is strictly positive if both $\omega_T$ and $\omega_P$ are strictly positive, so the corresponding initial fluctuation can a priori grow.

The asymptotic behaviour of $\omega^+$ is $2\omega_P$ at large $\omega_T/\omega_P$ and $2\omega_T$ at small $\omega_T/\omega_P$, as shown in figure 7. The growth rate is thus dominated by the slowest process, either the Poiseuille flow, driven by $\omega_P$, or the tension-induced recirculation, driven by $\omega_T$. As expected from §4.2, the instability disappears if one of the two processes becomes infinitely slow, i.e. if the corresponding characteristic rate ($\omega_P$ or $\omega_T$) vanishes. If both processes evolve on the same time scale, $\omega_T = \omega_P = \omega_{T,P}$ and the growth rate becomes $\omega^+ = (\sqrt{3} - 1)\omega_{T,P} \approx 0.73\omega_{T,P}$.

The negative eigenvalue $\omega^-$ scales as $-\omega_P$ and as $-\omega_T$, respectively, at small and large $\omega_T/\omega_P$ and is thus dominated by the fastest process.
# I. Cantat and C. Trégouët

## 5.2. System solution

### 5.2.1. Explicit time-dependent equations

Using expression (3.5a,b) for the reference bump and expression (3.11) for the line tension, (4.29) becomes

$$ \dot{\epsilon}_h = t^{-1} \epsilon_w, \qquad (5.3) $$

$$ \dot{\epsilon}_w = -t^{-1}\epsilon_w + t^{-3/4}K(2\epsilon_h - \epsilon_w), \qquad (5.4) $$

with $K = kI_1A^2/(4\zeta)$, $I_1$ a number defined in (3.11) and $A$ the invariant bump area, as defined in (3.5a,b).

To simplify the time dependency, we define a new time variable $\Theta = 4Kt^{1/4}$ so $t = (\Theta/(4K))^4$, and the associated unknown functions $X(\Theta) = \hat{\epsilon}_h(t)$ and $Y(\Theta) = \hat{\epsilon}_w(t)$ for the height and the width.

After the change of variable, (5.3) and (5.4) become

$$ \frac{d}{d\Theta} \begin{bmatrix} X \\ Y \end{bmatrix} = \begin{bmatrix} 0 & \Omega_p \\ 2\Omega_T & -\Omega_T - \Omega_p \end{bmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} \quad \text{with} \quad \Omega_p = \frac{4}{\Theta} \quad \text{and} \quad \Omega_T = 1. \qquad (5.5) $$

This system is formally identical to (5.1) and its two eigenvalues $\Omega^+$ and $\Omega^-$ are thus deduced from (5.2). At large $\Theta$, $\Omega^+$ scales as $2\Omega_p = 8/\Theta$ and $\Omega^-$ tends to $-\Omega_T = -1$.

### 5.2.2. A polynomial solution

As the positive eigenvalue tends to 0 at large $\Theta$, we guess the existence of a growing solution $(X^+, Y^+)$ with a power-law behaviour and we try to build a solution of the form

$$ X^+ = \Theta^N \sum_{j=0}^{\infty} a_j \Theta^{-j} \quad \text{and} \quad Y^+ = \Theta^N \sum_{j=0}^{\infty} b_j \Theta^{-j}, \qquad (5.6a,b) $$

satisfying, from (5.5),

$$ \frac{dX^+}{d\Theta} = \frac{4}{\Theta}Y^+ \quad \text{and} \quad \frac{dY^+}{d\Theta} = 2X^+ - \left(\frac{4}{\Theta} + 1\right)Y^+. \qquad (5.7a,b) $$

These conditions lead, for the terms of highest order ($j = 0$), to

$$ Na_0 = 4b_0 \quad \text{and} \quad 2a_0 - b_0 = 0, \qquad (5.8a,b) $$

which imposes $N = 8$.

For the generic terms $j \ge 0$, we thus get the conditions

$$ (8-j)a_j = 4b_j \quad \text{and} \quad 2a_{j+1} - b_{j+1} = b_j(12-j), \qquad (5.9a,b) $$

and after substitution of the first relation in the second, we get

$$ a_{j+1} = a_j \frac{(8-j)(12-j)}{1+j}. \qquad (5.10) $$

The coefficients $a_j$ are zero for $j > 8$ and a solution of (5.5) is the polynomial

$$ X^+ = \sum_{j=0}^{8} a_j \Theta^{8-j}; \quad Y^+ = \sum_{j=0}^{7} a_j \frac{8-j}{4} \Theta^{8-j}, \qquad (5.11a,b) $$
<figure><img src="image_10.png" /><figcaption>Stability of elongated thickness fluctuation in a soap film</figcaption></figure>

Figure 8. Numerical trajectories in the (X, Y) plane, with an initial condition ($X_0, Y_0$) = 0.1(cos($\theta$), sin($\theta$)) at $\Theta$ = 1, obtained by solving (5.5). All trajectories diverge at long time, except that obtained with $\theta$ $\approx$ 2.12[$\pi$] rad, corresponding to the stable solution (dashed lines). The arrows indicate the direction of the trajectories and the time $\Theta$ = 2.5 (arbitrarily chosen but visible on most of the trajectories) is indicated by a circle of the same colour as the curve. Note that the yellow circles are hidden by the blue ones, and green ones are not shown as they would be outside the figure.

with $a_0$ = 1 arbitrarily chosen and, from (5.10),

$$a_j = \frac{8!12!}{(8-j)!(12-j)!j!} \quad \text{and} \quad b_j = \frac{(8-j)8!12!}{4(8-j)!(12-j)!j!} \quad (5.12a,b)$$

This explicit solution of the problem diverges at large time and thus provides evidence that the invariant bump is unstable.

Another solution ($X^-, Y^-$), decreasing faster than exponentially, can also be built (see (A9) in Appendix A). All the solutions can be written as $(X, Y) = a^+(X^+, Y^+) + a^-(X^-, Y^-)$ with $a^+$ and $a^-$ two prefactors determined from the initial conditions at $\Theta_0 = 4Kt_0^{1/4}$, with $t_0$ the arbitrary initial time defined in § 4.1.

For the purpose of illustration, and to check our computations, the numerical integration of (5.5) has been performed with Matlab. The trajectories obtained for a set of initial conditions satisfying $X^2 + Y^2 = 0.01$ at the arbitrary time $\Theta_0 = 1$ are plotted in figure 8. It shows a divergence at large time, except for a particular initial condition, corresponding to the stable solution. The functions $X(\Theta)$ and $Y(\Theta)$ are shown in figure 9 as a function of time, for the same initial conditions as in figure 8. They have been divided by $X(6)$ to provide evidence that the large-$\Theta$ behaviour does not depend on the initial condition. Indeed, at large time, all the solutions are proportional to the polynomial solution $(X^+, Y^+)$ as expected.

### 5.3. Physical growth rate of the instability

In the following we consider only the contribution of the unstable mode (denoted with '+'), for $t > t_0$, and we come back to the physical variables, using the bar notation when needed. The explicit expression of the relative amplitude fluctuation of the unstable
I. Cantat and C. Trégouët

<figure><img src="image_11.png" /><figcaption>Figure 9. Numerical values of $X(\Theta)/X(6)$ (a) and $Y(\Theta)/X(6)$ (b), obtained by solving (5.5). The colour code is the same as in figure 8 and thus indicates the initial condition. The red circles represent the analytical solutions ($X^+/X^+(6)$, $Y^+/X^+(6)$), given by (5.11a,b).</figcaption></figure>

mode is directly deduced from the solution (5.11a,b), (5.12a,b), using the definition of the function $(X, Y)$ given in § 5.2.1:

$$ \hat{\varepsilon}_h^+ = \alpha \sum_{j=0}^{8} \frac{8!12!}{(8-j)!(12-j)!j!} \left( \frac{\bar{t}}{\tau(4K)^{-4}} \right)^{(8-j)/4}, \quad (5.13) $$

$$ \hat{\varepsilon}_w^+ = \alpha \sum_{j=0}^{7} \frac{8!12!}{(8-j)!(12-j)!j!} \frac{8-j}{4} \left( \frac{\bar{t}}{\tau(4K)^{-4}} \right)^{(8-j)/4}, \quad (5.14) $$

with $K = kI_1A^2/(4\zeta)$ a parameter proportional to the wavenumber $k$ and $\alpha$ a prefactor determined by the initial amplitude fluctuation, small enough to ensure that $\hat{\varepsilon}_h^+ \ll 1$. Note that, as we consider the unstable mode only, the initial values of $\hat{\varepsilon}_h^+$ and $\hat{\varepsilon}_w^+$ are not independent.

At large time, the solution behaves as

$$ \hat{\varepsilon}_h^+ \sim \alpha \left( \frac{\bar{t}}{\tau(4K)^{-4}} \right)^2, \quad \hat{\varepsilon}_w^+ \sim 2\alpha \left( \frac{\bar{t}}{\tau(4K)^{-4}} \right)^2. \quad (5.15a,b) $$

In order to get some orders of magnitude, we set the physical parameters to the values presented in table 1 (left). The chosen tension and viscosity correspond to usual surfactant solutions. In order to keep the viscous parameter $\zeta = \eta_s/(3\eta h_\infty)$ close to unity, we choose to consider an interface shear viscosity of $10^{-8}$ kg s$^{-1}$, which should be reached for highly soluble surfactants (Zell et al. 2014). Note that the friction of air may play a role for such a low interface viscosity, which has not been considered in the model (Lenavetier et al. 2024). With these values, the time scale and the viscous parameter are calculated and presented in table 1 (left). As initial bump shape, we choose the self-similar shape of height $\bar{h}^b(\bar{t}_0) = h_\infty\Phi(0)\delta h_0(t_0) = 200$ nm and of width $\bar{w}(\bar{t}_0) = h_\infty w_0(t_0) = 30$ µm. As the aspect ratio of the bump evolves with time following (3.5a,b), these choices impose the value of the initial time $t_0$ and of the geometrical parameter $A$, which are specified in table 1 (centre). Finally, the initial height perturbation is 1% of $\delta\bar{h}_0$.

The value of $\hat{\varepsilon}_h^+$ is plotted in figure 10(a) for the wavelengths $\bar{\lambda} = [1, 2, 5, 10]\bar{w}_0$, for the time range satisfying the constraint $\hat{\varepsilon}_h^+ < 0.3$. Note that only the last case satisfies the
$$
\begin{align*}
h_{\infty} &= 1 \, \mu\text{m} & \bar{h}^b(\bar{t}_0) &= 0.2 \, h_{\infty} = 200 \, \text{nm} & \hat{\epsilon}_h^+(\bar{t}_0) &= 0.01 \\
\gamma_0 &= 30 \times 10^{-3} \, \text{N m}^{-1} & \bar{w}(\bar{t}_0) &= 30 \, \mu\text{m} & & & \\
\eta &= 10^{-3} \, \text{kg m}^{-1} \, \text{s}^{-1} & t_0 &= (\bar{w}(\bar{t}_0)/h_{\infty})^4 & & & \\
\eta_s &= 10^{-8} \, \text{kg s}^{-1} & \delta h_0(t_0) &= (\bar{h}^b(\bar{t}_0)/h_{\infty})/\Phi(0) & & & \\
\tau &= \frac{3\eta h_{\infty}}{\gamma_0} = 10^{-7} \, \text{s} & w_0(t_0) &= 30 & & & \\
\zeta &= \frac{\eta_s}{\tau \gamma_0} = 3.33 & A &= w_0(t_0)\delta h_0(t_0) = 6/\Phi(0) = 20.8 & & & \\
& & K &= 1.88\bar{k}h_{\infty} = 11.81h_{\infty}/\bar{\lambda} & & & \\
& & \bar{t}_0 &= \tau w_0(t_0)^4 = 81 \, \text{ms} & & &
\end{align*}
$$

<figure><img src="image_12.png" /><figcaption>Figure 10. (a) Fluctuation amplitude divided by the bump height $\hat{\epsilon}_h(\bar{t})$, predicted by (5.13), for the physical parameters given in table 1. The wavelength is $\lambda = [1, 2, 5, 10]w_0$, respectively, for the magenta, red, green and blue curves. (b) Value of the different monomials in (5.13) for $\lambda = 2w_0$, from the $\bar{t}^{3/8}$ term (yellow) to the $\bar{t}^2$ term (dark brown). The monomials of lower orders are indistinguishable from zero. The dominant terms are the monomial $\bar{t}^{3/2}$ and $\bar{t}^{7/4}$, respectively, before and after $\bar{t} = 0.57$ s. The asymptotic power law, dominated by the $\bar{t}^2$ term, would only be reached after $\bar{t} = 22$ s.</figcaption></figure>

condition $\lambda \gg w_0$ assumed in the model, and that the other cases should be considered as extrapolated results. However, as the short-wavelength cases are the most unstable modes, the physical destabilisation process is expected to be close to these extrapolated modes. The contribution of each term of the polynomial (5.13) is plotted in figure 10(b). Low-order monomials dominate at short time, while high-order monomials only make a significant contribution at long time. For the example shown in figure 10, at the end of the linear regime at $\bar{t} - \bar{t}_0 \sim 0.5$ s, the largest contribution comes from the $\bar{t}^{7/8}$ term so the power law (5.15a,b) is never observed.

In order to characterise more precisely the bump shape, we define its maximal height $\bar{h}^{\text{max}}(\bar{t}) = h_{\infty}\Phi(0)\delta h_0(\bar{t})(1 + \hat{\epsilon}_h^+)$, reached at the point $(\bar{x} = 0, \bar{y} = 0)$, and the height at the saddle points along the bump crest $\bar{h}^{\text{sad}}(\bar{t}) = h_{\infty}\Phi(0)\delta h_0(\bar{t})(1 - \hat{\epsilon}_h^+)$, reached at the point $(\bar{x} = 0, \bar{y} = \bar{\lambda}/2)$. These characteristic heights are plotted in figure 11 as a function of time, for the same parameter and wavelength values as in figure 10. In this representation, the fate of the initially invariant bump clearly appears. The saddle-point height decreases faster than the average bump height, and the height of the maximum also decreases with
I. Cantat and C. Trégouët

<figure><img src="image_13.png" /><figcaption>Figure 11. Geometrical properties of the bump as a function of time in physical units. The wavelength is $\lambda = [1, 2, 5, 10]w_0$, respectively, for the magenta, red, green and blue curves. (a) The top curves are the maximal height $\bar{h}^{max}$ and the bottom curves the saddle-point height $\bar{h}^{sad}$. The black dashed line is the average bump height. (b) Similarly, the top and bottom curves are the maximal and minimal values of the bump width $\bar{w}$ and the black dashed line is the average bump width.</figcaption></figure>

time, but slower than the average bump height. Note that the thickness is only modified by the Poiseuille flow, which leads to a thickness decrease at every point along the bump crest: the increase of $\bar{h}^{max}$ observed at long time for the shortest wavelength is thus unphysical and comes from the fact that the system has left the linear regime. The difference between the maximum height and the saddle-point height $\delta\bar{h}^* = (\bar{h}^{max} - \bar{h}^{sad})/2 = h_\infty\Phi(0)\delta h_0 \hat{\epsilon}_h^+$ quantifies the absolute amplitude of the perturbation. As $\delta h_0 \sim \bar{t}^{-1/4}$ and $\hat{\epsilon}_h^+$ grows with an apparent power law in the range $\bar{t}^{[1.5-2]}$, the perturbation amplitude grows slightly faster than linearly.

The minimum and maximum widths of the perturbation are shown in figure 11. Importantly $\hat{\epsilon}_h^+$ and $\hat{\epsilon}_w^+$ have the same sign, so the width increases at the places of maximum height. A three-dimensional profile built from (5.13), (5.14) with parameters of table 1 is shown in figure 12(a,b) and shows this transition from the initial invariant bump towards well-separated round hills expected at long time (Trégouët & Cantat 2021). The shape evolution is also plotted for a groove in figure 12(c,d) to underline that the whole study remains valid for a negative value of $\delta h_0$, so a groove instead of a bump, as potentially relevant for the marginal regeneration problem.

Another important point is that the growth is faster at small wavelength. The model assumes $\lambda \gg w_0$, so the small-wavelength limit cannot be discussed in this frame. As a consequence, this model does not predict the fastest wavelength. A reasonable guess is that this fastest wavelength is of the order of the bump width.

6. Conclusion

This paper demonstrates that an elongated perturbation in the thickness of a soap film (whether a bump or a groove) destabilises into circular spots before eventually disappearing. More specifically, we show that the self-similar shape constructed by Benzaquen et al. (2013) is linearly unstable, and we provide the analytical expression of the unstable mode, whose growth over time follows a polynomial law. Our analytical predictions are restricted to the long-wavelength regime and therefore do not predict the fastest-growing wavelength. Determining this directly observable quantity would require numerical exploration. Experimental validation of this theory is highly appealing,
<figure><img src="image_14.png" /><figcaption>Figure 12. Examples of film-thickness profiles calculated with the parameters of table 1 (a,b) or with the same parameters but for a groove of initial depth $\bar{h}^b(\bar{t}_0) = -0.2h_\infty$ (c,d). The initial perturbations are shown in (a,c) and the profiles after 500 ms are shown in (b,d), for $\lambda = 20\bar{w}(\bar{t}_0)$. The colour bar represents $|\bar{h} - h_\infty|$ in micrometres.</figcaption></figure>

although preparing a controlled initial shape with the geometric properties discussed here may prove challenging. The most exciting prospect lies in extending the theory to scenarios involving a film in contact with a meniscus. Since all the relevant physical processes contributing to destabilisation remain unchanged in this slightly modified configuration, we assume that the instability causing marginal regeneration shares a common nature with the phenomenon discussed here. However, due to the non-self-similarity of the corresponding reference solution constructed by Aradian et al. (2001), the analytical treatment becomes potentially more complex. Thus, a numerical simulation, guided by the theoretical framework established in this work, appears to be a promising prospect.

**Acknowledgements.** We thank T. Salez, J. Eggers and J. Lister for enlightening discussions and S. Cantat for determining the stable solution and providing general mathematical advice.

**Funding.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 725094) and the ANR (DRAINFILM). It was partly carried out at the DAMTP, thanks to the hospitality of R. Goldstein.

**Declaration of interests.** The authors report no conflict of interest.

Author ORCIDs.

Isabelle Cantat https://orcid.org/0000-0002-3256-1254;

Corentin Trégouët https://orcid.org/0000-0002-2300-486X.

Appendix A

In order to determine the decreasing solution $(X^-, Y^-)$ of (5.5), we first use

$$Y = \frac{\Theta}{4}X' \quad \text{and} \quad Y' = \frac{\Theta}{4}X'' + \frac{X'}{4} \qquad (A1a,b)$$
to transform

$$Y' = 2X - \left(\frac{4}{\Theta} + 1\right) Y \qquad (A2)$$

into

$$X'' + \left(\frac{5}{\Theta} + 1\right) X' - \frac{8X}{\Theta} = 0. \qquad (A3)$$

We then introduce the auxiliary function $W$ defined so that $X^{-} = WX^{+}$, with $X^{+}$ the polynomial solution (5.11a,b). Equation (A3) becomes

$$X^{+}W'' + \left[ \left( \frac{5}{\Theta} + 1 \right) X^{+} + 2(X^{+})' \right] W' = 0. \qquad (A4)$$

Setting $Z = W'$ we obtain

$$-\frac{Z'}{Z} = \frac{5}{\Theta} + 1 + 2\frac{(X^{+})'}{X^{+}}, \qquad (A5)$$

leading to

$$\left.
\begin{aligned}
-\ln\left(\frac{Z}{Z^*}\right) &= \ln(\Theta^5) + \Theta + \ln((X^+)^2), \\
\frac{Z^*}{Z} &= \Theta^5 (X^+)^2 e^\Theta, \\
Z &= Z^* \frac{e^{-\Theta}}{\Theta^5 (X^+)^2},
\end{aligned}
\right\} \qquad (A6)$$

with $Z^*$ an arbitrary constant.

The functions $W$ and $X^{-}$ are thus

$$\left.
\begin{aligned}
W &= Z^* \int_{\Theta_0}^{\Theta} \frac{e^{-\Theta}}{\Theta^5 (X^+)^2} d\Theta + W(\Theta_0), \\
X^{-} &= X^+ \left( Z^* \int_{\Theta_0}^{\Theta} \frac{e^{-\Theta}}{\Theta^5 (X^+)^2} d\Theta + W(\Theta_0) \right).
\end{aligned}
\right\} \qquad (A7)$$

The constant $W(\Theta_0)$ can be determined from the constraint that $X^{-}$ tends to 0 at large $\Theta$. This imposes

$$W(\Theta_0) = -Z^* \int_{\Theta_0}^{\infty} \frac{e^{-\Theta}}{\Theta^5 (X^+)^2} d\Theta. \qquad (A8)$$

Finally we get, choosing $Z^* = -1$,

$$\left.
\begin{aligned}
X^{-} &= X^{+} \int_{\Theta}^{\infty} \frac{e^{-\Theta}}{\Theta^5 (X^{+})^2} d\Theta, \\
Y^{-} &= \frac{\Theta}{4} (X^{-})' = Y^{+} \int_{\Theta}^{\infty} \frac{e^{-\Theta}}{\Theta^5 (X^{+})^2} d\Theta - \frac{1}{4} \frac{e^{-\Theta}}{\Theta^4 X^{+}}.
\end{aligned}
\right\} \qquad (A9)$$
Stability of elongated thickness fluctuation in a soap film

REFERENCES

ARADIAN, A., RAPHAËL, E. & DE GENNES, P.-G. 2001 "Marginal pinching" in soap films. *Europhys. Lett.* **55**(6), 834.

BARIGOU, M. & DAVIDSON, J.F. 1994 Soap film drainage: theory of experiment. *Chem. Engng Sci.* **49** (11), 1807–1819.

BENZAQUEN, M., SALEZ, T. & RAPHAËL, E. 2013 Intermediate asymptotics of the capillary-driven thin-film equation. *Eur. Phys. J. E* **36**, 1–7.

BRUINSMA, R. 1995 Theory of hydrodynamic convection in soap films. *Phys. A* **216** (1–2), 59–76.

CHAMPOUGNY, L., SCHEID, B., RESTAGNO, F., VERMANT, J. & RIO, E. 2015 Surfactant-induced rigidity of interfaces: a unified approach to free and dip-coated films. *Soft Matt.* **11** (14), 2758–2770.

CHAN, D.Y.C., KLASEBOER, E. & MANICA, R. 2011 Film drainage and coalescence between deformable drops and bubbles. *Soft Matt.* **7**, 2235–2264.

CHATZIGIANNAKIS, E., JAENSSON, N. & VERMANT, J. 2021 Thin liquid films: where hydrodynamics, capillarity, surface stresses and intermolecular forces meet. *Curr. Opin. Colloid Interface Sci.* **53**, 101441.

CHATZIGIANNAKIS, E. & VERMANT, J. 2020 Breakup of thin liquid films: from stochastic to deterministic. *Phys. Rev. Lett.* **125** (15), 158001.

CHOMAZ, J.M. 2001 The dynamics of a viscous soap film with soluble surfactant. *J. Fluid Mech.* **442**, 387.

DENKOV, N.D., TCHOLAKOVA, S., GOLEMANOV, K., ANANTHAPADMANABHAN, K.P. & LIPS, A. 2008
Viscous friction in foams and concentrated emulsions under steady shear. *Phys. Rev. Lett.* **100**, 138301.

FROSTAD, J.M., TAMMARO, D., SANTOLLANI, L., BOCHNER DE ARAUJO, S. & FULLER, G.G. 2016 Dynamic fluid-film interferometry as a predictor of bulk foam properties. *Soft Matt.* **12**, 9266–9279.

HAMMOND, P.S. 1983 Nonlinear adjustment of a thin annular film of viscous fluid surrounding a thread of another within a circular cylindrical pipe. *J. Fluid Mech.* **137**, 363–384.

JOYE, J.-L., HIRASAKI, G.J. & MILLER, C.A. 1994 Asymmetric drainage in foam films. *Langmuir* **10** (9), 3174–3179.

LANDAU, L. & LEVICH, B. 1942 Dragging of a liquid by a moving plate. *Acta Physicochim. USSR* **17**, 42.

LENAVETIER, T., AUDÉOUD, G., BERRY, M., GAUTHIER, A., PORYLES, R., TRÉGOUËT, C. & CANTAT, I. 2024 Line tension in a thick soap film. *Phys. Rev. Lett.* **132** (5), 054001.

LHUISSIER, H. & VILLERMAUX, E. 2012 Bursting bubble aerosols. *J. Fluid Mech.* **696**, 5–44.

LISTER, J.R., RALLISON, J.M., KING, A.A., CUMMINGS, L.J. & JENSEN, O.E. 2006 Capillary drainage of an annular film: the dynamics of collars and lobes. *J. Fluid Mech.* **552**, 311–343.

MIGUET, J., PASQUET, M., ROUYER, F., FANG, Y. & RIO, E. 2021 Marginal regeneration-induced drainage of surface bubbles. *Phys. Rev. Fluids* **6** (10), L101601.

MYSELS, K.J., SHINODA, K. & FRANKEL, S. 1959 *Soap Films: Study of their Thinning and a Bibliography*. Pergamon.

NIERSTRASZ, V.A. & FRENS, G. 1998 Marginal regeneration in thin vertical liquid films. *J. Colloid Interface Sci.* **207** (2), 209–217.

NIERSTRASZ, V.A. & FRENS, G. 1999 Marginal regeneration and the marangoni effect. *J. Colloid Interface Sci.* **215**, 28–35.

ORON, A., DAVIS, S.H. & BANKOFF, S.G. 1997 Long-scale evolution of thin liquid films. *Rev. Mod. Phys.* **69** (3), 931.

PETIT, P., SEIWERT, J., CANTAT, I. & BIANCE, A.-L. 2015 On the generation of a foam film during a topological rearrangement. *J. Fluid Mech.* **763**, 286–301.

SAINT-JALMES, A. & TRÉGOUËT, C. 2023 Foam coarsening under a steady shear: interplay between bubble rearrangement and film thinning dynamics. *Soft Matt.* **19**, 2090–2098.

SEIWERT, J., DOLLET, B. & CANTAT, I. 2014 Theoretical study of the generation of soap films: role of interfacial visco-elasticity. *J. Fluid Mech.* **739**, 124–142.

SHABALINA, E., BÉRUT, A., CAVELIER, M., SAINT-JALMES, A. & CANTAT, I. 2019 Rayleigh–Taylor-like instability in a foam film. *Phys. Rev. Fluids* **4**, 124001.

SHI, X., FULLER, G.G. & SHAQFEH, E.S.G. 2022 Instability and symmetry breaking of surfactant films over an air bubble. *J. Fluid Mech.* **953**, A26.

STEIN, H.N. 1991 On marginal regeneration. *Adv. Colloid Interface Sci.* **34**, 175–190.

STUBENRAUCH, C. & KLITZING, R.V. 2003 Disjoining pressure isotherms of thin liquid films - new concepts and perspectives. *J. Phys.: Condens. Matter* **15**, R1197–R1232.

TRÉGOUËT, C. & CANTAT, I. 2021 Instability of the one-dimensional thickness profile at the edge of a horizontal foam film and its plateau border. *Phys. Rev. Fluids* **6** (11), 114005.

VELEV, O.D., CONSTANTINIDES, G.N., AVRAAM, D.G., PAYATAKES, A.C. & BORWANKAR, R.P. 1995 Investigation of thin liquid films of small diameters and high capillary pressures by a miniaturized cell. *J. Colloid Interface Sci.* **175** (1), 68–76.
I. Cantat and C. Trégouët

WEAIRE, D. & HUTZLER, S. 2000 *The Physics of Foams*. Oxford University Press.

ZELL, Z.A., NOWBAHAR, A., MANSARD, V., LEAL, L.G., DESHMUKH, S.S., MECCA, J.M., TUCKER, C.J. & SQUIRES, T.M. 2014 Surface shear inviscidity of soluble surfactants. *Proc. Natl Acad. Sci. 111* (10), 3677–3682.

ZHANG, Y. & SHARMA, V. 2018 Nanoridge formation and dynamics of stratification in micellar freestanding films. *Langmuir* **34** (3), 1208–1217.