# Surfactant-induced dissipation in sheared foams: mechanics and thermodynamics

Yedhir Mezache, François Detcheverry, Bastien Di Pierro, Peter D M Spelt,

Anne-Laure Biance, Marie Le Merrer

► To cite this version:

Yedhir Mezache, François Detcheverry, Bastien Di Pierro, Peter D M Spelt, Anne-Laure Biance, et al. Surfactant-induced dissipation in sheared foams: mechanics and thermodynamics. Physical Review Fluids, 2024, 9 (094008), 10.1103/PhysRevFluids.9.094008. hal-04249279v2

HAL Id: hal-04249279

https://hal.science/hal-04249279v2

Submitted on 9 Oct 2024

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
Surfactant-induced dissipation in sheared foams: mechanics and thermodynamics.

Yedhir Mezache,<fnref n="1" /> François Detcheverry,<fnref n="1" /> Bastien Di Pierro,<fnref n="2" /> Peter

D.M. Spelt,<fnref n="2" />, * Anne-Laure Biance,<fnref n="1" />, † and Marie Le Merrer<fnref n="1" />

Though the influence of surfactant type on foam rheological properties is well established exper-
imentally, the underlying physical mechanisms are far from understood. Here, using fully resolved
numerical simulation of an elementary T1 event taking into account both flow and surfactant dynam-
ics, we unveil the origin of surfactant-induced dissipation in sheared foams. To model the diffusive
and exchange contributions, we revisit the classical Lucassen model by switching from a mechanical
to a thermodynamic perspective and find that in spite of its extreme simplicity, it captures well some
behaviors of our numerical foam. Our approach should be useful in the class of surfactant-controlled
systems, from soap films and emulsions to coating films and spreading droplets.

<footnote n="1">¹Universite Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, UMR5306, F-69100, Villeurbanne, France.</footnote>

<footnote n="2">²University of Lyon, Université Claude Bernard Lyon 1, Ecole Centrale de Lyon, CNRS, Laboratoire de Mécanique des Fluides et d'Acoustique, F-69622 Villeurbanne, France.</footnote>

<footnote>* Deceased 2020</footnote>

<footnote>† anne-laure.biance@univ-lyon1.fr</footnote>
# I. INTRODUCTION

Foams are soft granular materials made of bubbles in a soapy solution. They are ubiquitous in everyday life and nature, from shampoos to sea foams and frog nests [1, 2]. Properties such as large specific area, low-cost and low-weight have also made them instrumental in a variety of applications such as cosmetics, fire extinguishers, templates for insulating building materials, or food texture modifiers among others [3]. Understanding the rheology of this class of materials is a long-standing quest [4, 5]. Just as for gels and suspensions, the origin of this complex behavior lies in their microstructure. Yet, foams are unique in two aspects. First, the bubbles are soft and highly deformable objects, a key property for their organization [6]. Second, these out-of-equilibrium materials involve surface-active species which adsorb at interfaces and induce repulsion forces that stabilize liquid films [7]. These surfactant molecules act as a third component that is essential to the properties of those biphasic systems.

It is now well-established that surfactants play a key role in the rheology of foams [8, 9]. A first line of evidence lies in interfacial rheology, where dedicated models and techniques, such as the oscillating bubble method [10, 11], could relate the surfactant mobility in the liquid to the interfacial complex modulus [12, 13], as first proposed by Lucassen [14]. More complex situations include slightly deformed two-bubble assembly or foam [15, 16], or other small deformation modes with viscous dissipation [17, 18]. The effect of the surfactant nature is also apparent in the case of large deformation, such as bubble neighbor switching as evidenced at the bubble scale [19–21] and in 3D foams [22]. The recurring observation is that the type of surfactant, though a component of molecular size, can deeply impact the flow behavior at a much larger scale. Yet, the underlying mechanisms have remained somewhat elusive because experiments typically probe only a global observable such as the loss modulus, and not the "hidden variables" [23] of complex surfactant dynamics.

In this work, we resolve the origin of surfactant-induced dissipation in a flowing foam and its dependency on surfactant properties. Using level-set simulations [24] to probe quantities usually inaccessible in experiments, we present a complete view for the processes at work in the elementary event of a sheared liquid foam. In particular, adopting a thermodynamic point of view, we have been able to decipher the origins of surfactant dissipation, whether it comes from surfactant diffusion or exchanges between the bulk and the interface. In addition, still relying on our thermodynamic approach, we theoretically extended the linear Lucassen model [14] to propose a tractable prediction for the surfactant-induced dissipation. Though it is a crude representation of a flowing foam, the model is able to rationalize several trends revealed by the simulations. Our approach is applicable to other surfactant-rich systems, from thin films to droplets or bubbly flows.

# II. DISSIPATION IN A NUMERICAL SHEARED FOAM.

## A. Numerical implementation: Governing equations and level-set method

The elementary event in a flowing foam is a so-called T1 event, wherein a bubble switches neighbors (Fig. 1a). We consider a minimal two-dimensional system of a few bubbles within a cell of height $H$ (Fig. 1b). Shear is imposed through pinned contact lines by walls moving at a velocity $U$. The liquid and gas flows are governed by the Navier-Stokes equation for incompressible fluids:

$$ \nabla \cdot \mathbf{u} = 0, \quad \rho \left[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right] = \nabla \cdot \bar{\sigma}, \quad (1) $$

where $\mathbf{u}$ is the fluid velocity field and $\bar{\sigma}$ is the stress tensor for an incompressible Newtonian fluid: $\bar{\sigma} = -p\bar{\bar{I}} + \eta(\nabla\mathbf{u} + (\nabla\mathbf{u})^T)$. The jump in the stress tensor $\bar{\sigma}$ across gas-liquid interfaces determines the boundary condition:

$$ [\bar{\sigma} \cdot \mathbf{n}] = -\gamma C \mathbf{n} - \nabla_s \gamma. \quad (2) $$

Here, the brackets denote a jump in quantity across the interface, $[x] = (x_{\text{liq}} - x_{\text{gas}})$, and $\mathbf{n}$ is the unit normal vector at interfaces, which points outwards from the gas to the liquid. $\gamma$ is the local surface tension and $C = -\nabla \cdot \mathbf{n}$ is the local interface curvature as seen from the liquid phase. $\nabla_s = \bar{I}_s \cdot \nabla$ is the surface gradient tensor, where $\bar{I}_s = \bar{\bar{I}} - \mathbf{n} \otimes \mathbf{n}$. The jump in normal stress is balanced by the normal component, $-\gamma C \mathbf{n}$, as given by the Young-Laplace equation. The tangential component, $\nabla_s \gamma$, balances the jump in traction across interfaces due to Marangoni stresses that arise from variations in surface tension along the interfaces.

The stress-jump boundary condition at the interface thus involves a Marangoni stress which depends on the surfactant nature and local concentration gradient. The concentrations of surfactants in the liquid and at the interface, $F$ and $f$ respectively, evolve through flow advection and through diffusion in the bulk and along the interface with
<figure><img src="image_2.png" /><figcaption>FIG. 1. Description of a sheared foam. (a) Schematic of a T1 event: upon shearing, bubble 1 and 4 become neighbors. (b) Numerical model of a sheared foam: the 4-bubble assembly sheared by plate motion is simulated with the level-set method. The dashed-line rectangle shows the minimal system of interest, for which periodic boundary conditions apply in the longitudinal direction. The liquid fraction is set to $\varphi = 0.3$.</figcaption></figure>

respective coefficients $D$ and $D_f$. Finally, the exchange of surfactants between the liquid and the interface is described by a Langmuir equation [23], with a flux $J = r_a F_s (f_\infty - f) - r_d f$; such that, $r_a$ and $r_d$ denote the adsorption and desorption rate coefficients, $f_\infty$ is the interfacial concentration at saturation, and $F_s$ the subsurface concentration.

At equilibrium where $J = 0$, the Langmuir adsorption isotherm is

$$ \frac{f_e}{f_\infty} = \frac{kF_e}{1 + kF_e}, \qquad (3) $$

where $k \equiv r_a/r_d$, and $kF_e$ is the so-called adsorption number [25]. The corresponding equation of state for the surface tension as a function of surfactant concentration is then

$$ \gamma(f) = \gamma_0 \left[ 1 + \frac{RTf_\infty}{\gamma_0} \ln \left( 1 - \frac{f}{f_\infty} \right) \right], \qquad (4) $$

where $R$, $T$, and $\gamma_0$ are respectively the ideal gas constant, the temperature, and the surface tension for clean interfaces ($f = 0$).

The surfactant properties are described by three dimensionless numbers: the bulk Péclet number $Pe \equiv UH/D$, the Biot number $Bi \equiv r_d H/U$ and the normalized adsorption length $h \equiv f_e/F_e H$, with $F_e$ and $f_e$ the values at equilibrium. These parameters may be thought of as characterizing respectively the surfactant diffusion, its sorption dynamics and the relative amount of surfactant lying at the interface with respect to the bulk. As a first approach, we use the size domain $H$ as a characteristic length of the system. Other possible choices for this characteristic length will be discussed in details below in Sec. IV.

The coupling between liquid flow, surfactant transport, interface boundary condition and the large interface deformation makes for an intricate problem that can only be solved numerically. We use a level-set method [26, 27] extended to account for the surfactant dynamics [28], which allows access to the fluid flow field and local surfactant distribution. The level-set scheme introduces a distance function $\phi$ that is advected by the flow:

$$ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = 0, \qquad (5) $$

where the interfaces are captured at each time step wherever $\phi = 0$. All governing equations, definitions of dimensionless quantities and default parameters are further detailed in the Supplementary Material (SM [29], section I).

## B. Numerical results: Local and global dissipation

Fig.2 shows a typical view of the various stages of the T1 rearrangement, once an oscillatory state is reached. Indeed, as shown in Fig. 3, various computed quantities such as the injected power or the viscous dissipation reach a steady state after two oscillations. The corresponding time period $H/U\sqrt{3}$ is determined by the hexagonal bubble arrangement and the relative velocity between facing bubbles.

When the bubble assembly is sheared, the interface may be locally compressed or dilated, inducing enrichment or depletion in surfactants, as can be seen in Fig.2a. These periodic changes in interfacial area are illustrated in Fig. 3a, where we report the quantity $\int_{\Gamma} (\nabla_s \cdot \mathbf{u})^2 d\Gamma$, with $\Gamma$ the interface contour and $\nabla_s$ the surface divergence, which accounts for the local compression or dilation of the surface. The local variation in interfacial area induces
<figure><img src="image_3.png" /><figcaption>FIG. 2. Locally resolved T1 event. Time increases from left to right, corresponding to $Ut/H = 1.8, 2.1, 2.4$ and $2.7$ successively. Top: surfactant concentration in the liquid and at interface ($F$ and $f$ resp.). Bottom: local rate of dissipation due to surfactant diffusion ($D_{\text{dif}}^{\text{loc}}$) and to exchange processes ($D_{\text{ex}}^{\text{loc}}$). The simulation parameters are $Pe = 1, Bi = 10$ and $h = 1$. The capillary number $\eta U/\gamma_e$, with $\eta_l$ the liquid viscosity and $\gamma_e$ the surface tension at equilibrium, is 0.1 in all simulations.</figcaption></figure>

exchanges of surfactant with the adjacent liquid resulting in concentration heterogeneities in the bulk, which are counteracted by the smoothing effect of diffusion and possibly by flow advection.

To identify dissipation sources during the T1 process, we use an energy balance. Fig. 3b shows the temporal evolution of the injected power, evaluated from forces applied on the wall [28], from which we calculate the time-average $P_{\text{inj}}$, taken during the last oscillation, once the steady-state is reached. The local viscous dissipation reads $D_v^{\text{loc}} \equiv \eta \nabla u : (\nabla u + \nabla u^T)$, and we report in Fig. 3b this spatially integrated viscous dissipation as a function of time. We denote as $D_v \equiv \langle D_v^{\text{loc}} \rangle$ its time-averaged value after the steady-state is reached, with $\langle \cdot \rangle$ indicating both a spatial integration over the whole system and a time average over the period of one T1 rearrangement. As expected, $D_v$ is found to evolve with the surfactant-related parameters, due to surfactant-induced Marangoni stresses. We also observe that $\langle D_v \rangle$ is systematically smaller than $P_{\text{inj}}$, with deviations of up to 30%. Since we are in steady-state, the rate of change of kinetic energy is zero when averaged over one oscillation (see one example in Fig. S3 in SM [29]). The difference $\langle D_s \rangle \equiv P_{\text{inj}} - \langle D_v \rangle$ is thus the surfactant-induced dissipation [28] and is the focus of this work.

To understand its microscopic origins, we switch from the classical mechanical approach to a thermodynamic framework. According to non-equilibrium thermodynamics [30], the dissipation rate is proportional to the rate of entropy

<figure><img src="image_4.png" /><figcaption>FIG. 3. Temporal evolution of global quantities. The time unit is $H/U$. (a) Variations in $\int_{\Gamma} (\nabla_s \cdot \mathbf{u})^2 d\Gamma$ (the time-average of which gives the parameter $C$). (b) Variations of injected power, total viscous dissipation, surfactant-induced dissipation and surface work $\int_{\Gamma} \gamma \nabla_s \cdot \mathbf{u} d\Gamma$. The horizontal lines show the values averaged over the last oscillation. The simulation parameters are $Pe = 1, Bi = 10$ and $h = 0.1$ (case 1 in Fig. 4). See also Fig. S2 in SM [29].</figcaption></figure>
production, where each source contributes through the product of a thermodynamic force and an associated flux. In a sheared foam, the viscous dissipation in the flow is supplemented by two contributions coming from the surfactants. The first is the diffusion contribution, present whenever there are gradients in surfactant concentration [31]; the second is the exchange contribution deriving from adsorption/desorption processes at the interface. Their local dissipation rates, per unit volume and per unit interfacial area respectively, are (see SM [29], section II.B for details):

$$ \mathcal{D}_{\text{dif}}^{\text{loc}} = \nabla \mu_F (D \nabla F) = DRT (\nabla F)^2 / F, \qquad (6a) $$

$$ \mathcal{D}_{\text{ex}}^{\text{loc}} = (\mu_{F_s} - \mu_f)J \approx RTJ^2/(r_d f_e), \qquad (6b) $$

with R the ideal gas constant and T the temperature. In Eq. (6a), the last equality assumes that the chemical potential of the surfactants in the liquid is $\mu_F = \mu_0 + RT \ln F/F_e$. In Eq. (6b), the thermodynamic force is the difference in chemical potentials between the surface ($\mu_f = \mu_0 + RT \ln [f/(f_\infty - f)]$) and the adjacent liquid ($\mu_{F_s} = \mu_F(F = F_s)$). The last equality is valid close to equilibrium. An example map of these local dissipation rates is shown in Fig.2b.

We first check the consistency of our description. Besides its definition as $\langle D_s \rangle \equiv P_{\text{inj}} - \langle D_v \rangle$, $\langle D_s \rangle$ can also be expressed thermodynamically as $\langle D_{\text{dif}}^{\text{loc}} \rangle + \langle D_{\text{ex}}^{\text{loc}} \rangle$. A third alternative formulation from purely mechanical considerations [28] is $\langle \gamma \nabla_s \cdot \mathbf{u} \rangle$ (where the spatial integration is made on the interface contour). Note that $\langle \gamma \nabla_s \cdot \mathbf{u} \rangle$ represents the rate of change of surface energy, which is zero on average for surfactant-free interfaces, but positive in presence of surfactants [28]. We find in Fig. 3b (see also Fig. S2 in SM [29]) that the three numerical estimates for $\langle D_s \rangle$ match each other, typically within 10% due to numerical discretization errors reflected in the error bars of Fig. 4a. This agreement confirms that no significant source of surfactant dissipation has been left out.

As expected, the surfactant properties do significantly influence the dissipation. The total surfactant dissipation is shown in Fig. 4 for various combinations of Pe, Bi and h, which characterize the surfactant dynamics. We observe a highly non-trivial behaviour, with a maximum dissipation that appears at intermediate values of h. The rationale for this behavior is difficult to pinpoint at this stage.

Whether surfactant dissipation is dominated by diffusion or exchange processes depends on the surfactant properties. The ratio $\langle D_{\text{dif}} \rangle / \langle D_{\text{ex}} \rangle$ obtained numerically for a variety of Pe, Bi and h values varies by several orders of magnitude (Fig. 4a). Is the difference reflected in the maps of surfactant distribution and dissipation? Consider first the exchange dominated regime, as in case 1 of Fig. 4. Here, the surfactant distribution is very homogeneous both in the bulk and at the interface. Consequently, dissipation due to diffusion across the liquid bulk is weak compared to dissipation generated by surfactant exchanges with the interface. The situation is reversed upon increasing Péclet to Pe = 100, i.e. taking a surfactant with much lower diffusivity (case 2 of Fig. 4). The diffusion contribution becomes dominant because the surfactants in the bulk and at the interface are heterogeneously distributed. Dissipation due to diffusion is then more pronounced and localized in the vicinity of the interfaces. Finally, a different dissipation pattern can be observed by varying h while keeping the same ratio $\langle D_{\text{dif}} \rangle / \langle D_{\text{ex}} \rangle$ (case 3 of Fig. 4). Here the concentration and dissipation heterogeneities are parallel to the interface, instead of being normal as in case 2. Overall, we observe that how and where dissipation occurs depends sensitively on surfactant properties.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Simulation results. (a) Surfactant dissipation $\langle D_s \rangle$ as a function of h for various Pe and Bi. The simulation unit for the dissipation is $\rho_l U^3 H$ in our 2D geometry, with $\rho_l$ the liquid density. (b) Ratio $\langle D_{\text{dif}} \rangle / \langle D_{\text{ex}} \rangle$ as a function of h for various Pe and Bi, as obtained from simulations. (c) Maps of local concentration (top) and dissipation (bottom) for cases 1-2 and 3 circled in (b). The color scales are the same as in Fig. 2.</figcaption></figure>
### III. EXTENDED LUCASSEN MODEL

#### A. Model and approximations

In our simulations, the bubble motion generates local variations of the interfacial area, which vary periodically in time, as observed from the evolution of the parameter $\int_{\Gamma}(\nabla_s \cdot \mathbf{u})^2 d\Gamma$ in Fig. 3a. To shed light on our numerical observations, we revisit the 60-year-old model of Lucassen and van den Tempel [14] who first rationalized how surfactant dynamics could induce an effective surface viscosity [32]. The model focuses on the surfactant-induced dissipation that arises when imposing periodic variations in interfacial area. Specifically, for the simplified geometry of Fig. 5, we study how a liquid film of finite thickness $2W$ and unit depth that contains surfactants responds to an imposed area variation $\epsilon = \delta A/A = \epsilon_0 \exp(i\omega t)$, where $\epsilon_0$ is a small parameter and $\omega$ the forcing frequency.

The model is approximate in several respects. First, it focuses only on surfactant dissipation, and therefore neglects many ingredients present in the simulations, such as: the liquid flow and viscous dissipation, advection of surfactants by this liquid flow, as well as the complex geometry and the heterogeneity of surface compression and dilation. Second, the model only considers the linear regime, meaning that the interfacial area variations are small compared to the total area of the system. In our simulations, the amplitude of the area variations can be estimated by the typical values of $\int_{\Gamma} \text{div}_s^2 u \, d\Gamma$. In the Lucassen configuration, this would correspond to $A\epsilon^2\omega^2 \sin^2\omega t$. The typical relative deformation is then such that $\epsilon^2 \sim \int_{\Gamma} \text{div}_s^2 u \, d\Gamma / (A\omega^2)$. Taking the values of Fig. 3, the initial total area (per unit depth) $A \approx 3.2H$ and the dimensionless frequency $\tilde{\omega} = 2\pi/\sqrt{3}$, this yields $\epsilon \sim 0.1$, which suggests a linear model is appropriate.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Extended Lucassen model. A soap film of thickness $2W$ has its area $A$ subject to sinusoidal relative variations $\epsilon = \delta A/A = \epsilon_0 \exp(i\omega t)$. For symmetric forcing, the area variations of top and bottom interfaces are identical, while for antisymmetric forcing, they are opposite.</figcaption></figure>

Our approach differs from the original model in two ways: (i) Instead of instantaneous exchanges and an infinite reservoir, we assume both a finite adsorption/desorption kinetics and a finite thickness for the film. Though these extensions were considered separately [13, 33], we treat the general case here. The thin film system considered here has two interfaces, which may be excited in a different manner. We will consider the two cases of a symmetric forcing where both interfaces are compressed and dilated simultaneously, and an antisymmetric forcing, with one interface that is dilated while the other is compressed. (ii) More importantly, we adopt a complementary thermodynamic view where we compute not only the loss modulus that encapsulates all dissipative relaxation processes, but also the separate contributions from each dissipation mechanism, allowing a measure of their relative importance. The calculations for this extended Lucassen model are fully detailed in section II of SM [29].

#### B. Results: diffusion and exchange dissipations

The main result is a set of explicit formulas for dissipation induced by diffusion and exchange of surfactants between liquid and interface, together with the definition of an elastic complex modulus of the interface $E^*$. In what follows, we use characteristic length and time scales $W$ and $1/\omega$ to introduce the dimensionless numbers $Pe \equiv \omega W^2/D$, $Bi \equiv r_d/\omega$ and $h = f_e/F_e W$. Note that, even though these dimensionless numbers differ from those used above for numerical simulations, we use the same symbols for convenience.
## 1. Symmetric forcing

Let us first consider the case of a symmetric forcing. Focusing on time and spatial average for simplicity, we find

$$ \frac{\langle \mathcal{D}_{\text{dif}} \rangle}{\mathcal{D}_{\text{un}}} = \frac{\bar{\chi} h \sqrt{8\text{Pe}}}{|\zeta^{-1} + 1 - i|^2} \mathcal{G}(\sqrt{2\text{Pe}}), \qquad \mathcal{G}(u) \equiv \frac{\sinh(u) - \sin(u)}{\cosh(u) - \cos(u)}, \qquad (7a) $$

$$ \frac{\langle \mathcal{D}_{\text{ex}} \rangle}{\mathcal{D}_{\text{un}}} = \frac{4\bar{\chi}}{\text{Bi}|\zeta^{-1} + 1 - i|^2}, \qquad \frac{\zeta^{-1}}{\sqrt{2}\bar{\chi}} \equiv h\sqrt{\text{Pe}} \coth\left(j\sqrt{\text{Pe}}\right) + \frac{j}{\text{Bi}}, \quad j \equiv \frac{1+i}{\sqrt{2}}, \qquad (7b) $$

with $1 - \bar{\chi} \equiv \chi \equiv f_e/f_\infty$ a parameter for surface coverage that remains constant. The left-hand sides are made dimensionless with $\mathcal{D}_{\text{un}} \equiv E_{\text{GM}}\omega A\epsilon_0^2/2$, where $E_{\text{GM}}$ is the Gibbs-Marangoni elasticity modulus defined as $-d\gamma/d(\ln f)|_e$ evaluated at thermodynamic equilibrium [12]. For the Langmuir equation of state that we use here (Eq. (4)), this gives $E_{\text{GM}} = RTf_e/\bar{\chi}$. Thus, from the physico-chemical properties of the surfactants as specified through the dimensionless numbers, Eqs. (7a) and (7b) can predict the origin and magnitudes of surfactant dissipation induced by deformations of the interface. In particular, we can identify which mechanism dominates the surfactant dissipation. In this regard, the prediction of the model is very simple:

$$ \frac{\langle \mathcal{D}_{\text{dif}} \rangle}{\langle \mathcal{D}_{\text{ex}} \rangle} = \text{Da} \mathcal{G}(\sqrt{2\text{Pe}}), \qquad (8) $$

with $\mathcal{G}$ defined in Eq. (7a), and $\text{Da} \equiv \bar{\chi} h \text{Bi} \sqrt{\text{Pe}/2}$ the Damköhler number, which compares sorption and diffusion timescales.

The original Lucassen model adopts a mechanical perspective that focuses on the complex modulus. We derive also the complex viscoelastic modulus of the film in our extended Lucassen approach (section II.C in the SM [29]). It is defined as

$$ E \equiv E' + iE'' = \frac{2d\gamma}{d(\ln A)} = \frac{-2d\gamma}{d(\ln f)} \times \frac{-d(\ln f)}{d(\ln A)}. \qquad (9) $$

The factor 2 here accounts for the two interfaces of the film. When using the $\zeta$ variable defined in Eq. (7b), the complex modulus take a form similar to that obtained in the original Lucassen model for a single interface, namely

$$ \frac{E}{E_{\text{GM}}} = \frac{1}{1 + (1 - i)\zeta} = \frac{1 + \zeta + i\zeta}{1 + 2\zeta + 2\zeta^2}. \qquad (10) $$

Once rewritten in terms of the dimensionless numbers Pe, Bi, $h$ and $\chi$, and accounting for the two interfaces of the film, the complex modulus of the film in the extended Lucassen model reads as:

$$ \frac{E}{E_{\text{GM}}} = 2 \left[ 1 + \frac{1}{\bar{\chi} \left[ h(1+i)\sqrt{\text{Pe}/2} \coth\left( (1+i)\sqrt{\text{Pe}/2} \right) + i/\text{Bi} \right]} \right]^{-1}. \qquad (11) $$

Since the loss modulus $E''$, the imaginary part of $E$, is a measure of the total energy dissipated per cycle of sinusoidal deformation, it should satisfy

$$ \langle \mathcal{D}_{\text{ex}} \rangle + \langle \mathcal{D}_{\text{dif}} \rangle = \frac{\epsilon_0^2 \omega E'' A}{2} \qquad (12) $$

We have checked analytically that Eq. (12) is indeed satisfied, which confirms the consistency between mechanical and thermodynamical approaches.

## 2. Antisymmetric forcing

The case of an antisymmetric forcing, depicted in Fig. 5b, can be treated in a similar way. Here the interfaces are not compressed and dilated simultaneously but one interface is compressed while the other is dilated. Antisymmetric forcing is considered for two reasons. First, we want to assess the impact of the forcing type on the resulting dissipation.
Second, when inspecting the simulation results, such as in Fig. 4c for instance, interfaces facing each other might be subject to antisymmetric dilation and compression. We obtain (section II.E in SM [29]):

$$ \frac{\langle \mathcal{D}_{\text{dif}} \rangle}{\mathcal{D}_{\text{un}}} = \frac{\bar{\chi} h \sqrt{8\text{Pe}}}{|\bar{\zeta}^{-1} + 1 - i|^2} \bar{\mathcal{G}}(\sqrt{2\text{Pe}}), \qquad \bar{\mathcal{G}}(u) \equiv \frac{\sinh(u) + \sin(u)}{\cosh(u) + \cos(u)}, \qquad (13a) $$

$$ \frac{\langle \mathcal{D}_{\text{ex}} \rangle}{\mathcal{D}_{\text{un}}} = \frac{4\bar{\chi}}{\text{Bi}|\bar{\zeta}^{-1} + 1 - i|^2}, \qquad \frac{\bar{\zeta}^{-1}}{\sqrt{2}\bar{\chi}} = h\sqrt{\text{Pe}}\tanh(j\sqrt{\text{Pe}}) + \frac{j}{\text{Bi}}, \quad j \equiv \frac{1+i}{\sqrt{2}}. \qquad (13b) $$

Note that the formulas are almost identical to the symmetric case, except for the tanh function in $\bar{\zeta}$ and the plus signs in $\bar{\mathcal{G}}$. The ratio between diffusion and exchange dissipation now reads:

$$ \frac{\langle \mathcal{D}_{\text{dif}} \rangle}{\langle \mathcal{D}_{\text{ex}} \rangle} = h \text{Bi} \sqrt{\text{Pe}/2} \bar{\mathcal{G}}(\sqrt{2\text{Pe}}). \qquad (14) $$

Finally, the complex modulus can also be computed from the mechanical route, giving $E/E_{\text{GM}} = 2/(1 + (1 - i)\bar{\zeta})$, in complete analogy with the symmetric case. We checked again the consistency between the thermodynamic and mechanical approaches.

### C. Discussion

We now discuss some specific features of the dissipations in the two cases. We first consider the total dissipation and then discuss the dominant mechanism of dissipation.

#### 1. Total surfactant dissipation and h dependence

For convenience, we consider here the dimensionless dissipation defined as $\langle \tilde{\mathcal{D}}_s \rangle = E''/E_{\text{GM}}$. The dependence of dissipation on normalized adsorption length $h$ is illustrated in Fig. 6 for the two cases considered and for some values of Bi and Pe. Two notable features are visible: the first is a plateau at low $h$, whose magnitude is controlled by Bi, the second is a maximum at a position $h_{\text{max}}$ controlled by Pe. To better understand such a behavior, we consider $\langle \mathcal{D}_s \rangle$ in the $(h, \text{Pe})$ plane, as plotted in Fig. 7 for various Bi numbers and identify three simple features.

(i) Plateau at low $h$ and Pe. For small $h$, the dissipation approaches the constant

$$ h \to 0, \qquad \langle \tilde{\mathcal{D}}_s \rangle = \frac{2\bar{\chi}\text{Bi}}{\bar{\chi}^2 + \text{Bi}^2}. \qquad (15) $$

In this case, most surfactant lays within the liquid film and its distribution is unaffected by driving. Dissipation originates only in exchange processes and is therefore identical for symmetric and anti-symmetric driving. For $\text{Pe} \to 0$, Eq. (15) still applies in the antisymmetric case, while the symmetric case leads to a slightly different value [34].

(ii) Decay at high $h$ and Pe. Dissipation vanishes for large $h$ and Pe numbers. Specifically, $\langle \tilde{\mathcal{D}}_s \rangle \sim h^{-1}$ for $h \to \infty$ and $\langle \tilde{\mathcal{D}}_s \rangle \sim \text{Pe}^{-1/2}$ for $\text{Pe} \to \infty$.

(iii) Maximum. Let us focus first on the high Pe regime. For sufficiently high Bi and whatever the type of driving, a maximum in $\langle \tilde{\mathcal{D}}_s \rangle(h)$ is observed, whose position and value are given by

$$ h_{\text{max}} = \frac{1}{\bar{\chi}\sqrt{\text{Pe}}} \left[ 1 - (1 + \sqrt{2})\bar{\chi}/\text{Bi} \right], \qquad \langle \tilde{\mathcal{D}}_s \rangle_{\text{max}} = \frac{1}{(1 + \sqrt{2})(1 - \bar{\chi}/\text{Bi})}. \qquad (16) $$

The position of this maximum is shown with solid line in Fig. 7. Now what is the physical origin of this maximum? It is useful at this point to introduce the dimensional adsorption length $l_a \equiv \bar{\chi}f_e/F_e$ and the thickness of the diffusive boundary layer $|\kappa^{-1}| \equiv \sqrt{D/\bar{\omega}} = W/\sqrt{\text{Pe}}$. Taking for simplicity the high Bi limit, the dissipation maximum occurs when $|\kappa^{-1}| = l_a$, that is when the diffusion acts over a typical depth that is precisely the length characterizing surfactant distribution between film and interface.

Consider next the low Pe regime. It is apparent from Fig. 7 that the maximum may be absent (top middle figure) in the symmetric driving but that it survives in the antisymmetric case (dashed line). Specifically, at lowest order in Pe, one finds

$$ h_{\text{max}} = \frac{1}{\bar{\chi}\text{Pe}}, \qquad \langle \tilde{\mathcal{D}}_s \rangle_{\text{max}} = 1 - \frac{\text{Pe}}{3}. \qquad (17) $$
The condition for $h_{max}$ can be rewritten as $|\kappa^{-1}| \equiv \sqrt{l_a W}$. The surfactant dissipation is maximal when the thickness of diffusive layer is the geometric average of adsorption length and film thickness. Presumably, the later becomes important in antisymmetric forcing because diffusion across the film is now relevant.

Finally, we have taken so far the three dimensionless parameters governing our problem (h, Bi, Pe) as independent quantities, because this is natural from a theoretical point of view. However, in a typical experiment, the parameter that would most easily be varied while keeping all others fixed would be the driving frequency $\omega$, which would modifies Bi and Pe simultaneously. In this case, we find again a maximum in $\langle \tilde{D}_s \rangle(\omega)$. Indeed, $\langle \tilde{D}_s \rangle(\omega)$ is increasing at low $\omega$ and increasing at high $\omega$, specifically:

$$ \omega \to 0, \quad \langle \tilde{D}_s \rangle = \frac{2l_a(1 + k_F W/3D)}{k_F(1 + l_a/W)^2} \omega, \qquad \omega \to \infty, \quad \langle \tilde{D}_s \rangle = \frac{2k_F}{l_a \omega}, \qquad (18) $$

where $k_F \equiv r_a f_e(\chi^{-1} - 1)$. The exact position of the maximum cannot be written explicitly but can be approximated from the crossing of the two limiting regimes above, giving

$$ \omega_{\max} = \frac{k_F(l_a^{-1} + W^{-1})}{\sqrt{1 + k_F W/3D}}, \qquad (19) $$

The pulsation $\omega_{\max}$ is controlled either by the film half-thickness $W$ or adsorption length $l_a$, whichever is smallest, with a prefactor that is simply $k_F$ in the limit of low $k_F$ and $\sqrt{3k_F D/W}$ in the limit of high $k_F$.

## 2. Origin of dissipation

We now discuss how the dissipation induced by exchanges $\langle D_{ex} \rangle$ compares with the one induced by diffusion $\langle D_{dif} \rangle$. We find from Eqs. (8) and (13a)-(13b) that, for both types of forcing, the ratio $\langle D_{dif} \rangle / \langle D_{ex} \rangle$ scales with the Damköhler number $Da = \bar{\chi} h \sqrt{Pe/2}$, with a prefactor that only depends on the Péclet number: $G(\sqrt{2Pe})$ (resp. $\bar{G}(\sqrt{2Pe})$) for the symmetric (resp. antisymmetric) forcing. These functions are plotted in Fig. 8 with Pe as the argument. We observe that both $G$ and $\bar{G}$ are approximately unity (within 10%) provided that $Pe \ge 5$. If we now consider a mixed model, which averages the contributions from symmetric and antisymmetric forcing, we recover the same scaling in Da, but the prefactor is now a linear combination of $G$ and $\bar{G}$, namely

$$ \mathcal{G}_{\text{eff}} = \Delta\mathcal{G}(\sqrt{2\text{Pe}}) + (1 - \Delta)\bar{\mathcal{G}}(\sqrt{2\text{Pe}}), \qquad (20) $$

with $\Delta$ a function of Pe, Bi and $h$ [35]. Figure 8 shows how $\mathcal{G}_{\text{eff}}$ changes over the range of dimensionless parameters of interest in this work. In the mixed model, we observe that the effective ratio is close to unity, with deviation below 10%, as soon as $Pe \ge 1$. Overall, whatever the exact configuration (symmetric, antisymmetric or mixed forcing), we find that the ratio $\langle D_{dif} \rangle / \langle D_{ex} \rangle$ is well predicted by the Damköhler number $Da = \bar{\chi} h \sqrt{Pe/2}$, provided the Péclet number is larger than unity.

Finally, for Pe below unity, $G(\sqrt{2Pe})$ and $\bar{G}(\sqrt{2Pe})$ are proportional to $\sqrt{2Pe}$, with a prefactor 1/3 and 1 respectively. The ratio $\langle D_{dif} \rangle / \langle D_{ex} \rangle$ is then controlled by the product $\bar{\chi} h \text{BiPe}$.

<figure><img src="image_7.png" /><figcaption>FIG. 6. Prediction from extended Lucassen model. Normalized surfactant dissipation for the film $\langle \tilde{D}_s \rangle(h)$ for Pe = 0.1, 1, 100 and Bi = 1, 10 in the (a) symmetric and (b) antisymmetric forcing.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. Density plot of the surfactant dissipation ($\langle D_s \rangle$)($h$, Pe) for Bi = 1, 10, 100 (left to right) in the symmetric and antisymmetric forcing (top and bottom resp.). The black line indicates the approximate expressions for maximum position $h_{max}$ at high Pe (solid line, Eq. (16)) and at low Pe (dashed line, Eq. (17)).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Prefactor $\langle D_{dif} \rangle / \langle D_{ex} \rangle$ / Da arising in the symmetric model ($G$, black), antisymmetric model ($\bar{G}$, grey) and mixed model ($G_{eff}$, colors), as defined from Eq. (20) and for various values of $h$ and Bi.</figcaption></figure>

## IV. DISCUSSION: COMPARISON OF EXTENDED LUCASSEN MODEL TO NUMERICAL RESULTS

To compare the surfactant dissipation in the model and in simulations, we first need to match the dimensionless parameters from the two situations. To do so, we introduce the dimensionless pulsation $\tilde{\omega} = \omega H/U$ and equivalent film half-thickness $\tilde{W} = W/H$. The correspondence between the relevant dimensionless parameters is then:

$$Pe_L = \tilde{W}^2 \tilde{\omega} Pe, \quad Bi_L = \frac{Bi}{\tilde{\omega}}, \quad h_L = \frac{h}{\tilde{W}}, \qquad (21)$$

where the simulation (resp. model) dimensionless parameters are $Pe$, $Bi$ and $h$ (resp. $Pe_L$, $Bi_L$ and $h_L$).
<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison between simulation results and extended Lucassen model. Rescaled surfactant dissipation $\langle \tilde{D}_s \rangle$, equivalent to the loss modulus or effective viscosity of the surface, as a function of $h$. (a) Simulation results rescaled with $\tilde{\omega} = 3$. (b) Mixed extended Lucassen model (average of symmetric and antisymmetric models) with $\tilde{\omega} = 3$ and $\tilde{W} = 0.5$. facteur 2 du modèle à expliquer</figcaption></figure>

Our first comparison between the simulations and model focuses on the total surfactant dissipation. For the comparison to be meaningful, the use of a re-scaled dissipation $\langle \tilde{D}_s \rangle$ is required. For the model, we consider $\langle \tilde{D}_s \rangle \equiv \langle D_s \rangle / 2D_{un}$, where $D_{un}$ was introduced in Eq. (7a). We consider here the dissipation per interface (and not per film area), hence the factor 1/2. For the simulations, $\langle \tilde{D}_s \rangle \equiv \xi \tilde{\omega} \langle D_s \rangle / C$, where $\omega = \tilde{\omega} U/H$, $\xi$ is a parameter that is kept constant in our simulations, and $C \equiv \langle (\nabla_s \cdot \mathbf{u})^2 \rangle$ (where the spatial integration is along interfaces). On the one hand, we define $\xi = \text{Re} C a_0 \bar{\chi} / \beta \chi$, where the Reynolds number $\text{Re}$, capillary number $C a_0$ and $\beta$ parameters are defined in section I.C of SM [29]. It accounts for the fact that the simulation unit for energy is $\rho_l U^2 H^3$, with $\rho_l$ the liquid density. On the other hand, the parameter $C$ accounts for a heterogeneous dilatation along the interface and varies with the surfactant-related parameters $h$, Pe and Bi (section I.F in SM [29]). Both renormalizations are equivalent in the extended Lucassen geometry where the interfacial deformation is homogeneous (section III in SM [29]). The re-scaled dissipation therefore corresponds to the dimensionless loss modulus or effective viscosity of the surface.

We now discuss the choice of $\tilde{\omega}$ and $\tilde{W}$ prefactors. The simplest expectation would be plain geometric matching. As observed above in Fig. 3, the rearrangement period is $H/(U\sqrt{3})$, yielding $\tilde{\omega} = 2\pi/\sqrt{3}$. The choice of $\tilde{W}$ is less obvious. Indeed, in the simulations, the film thickness is an ill-defined quantity, which varies with space, time and dimensionless parameter values. One could however assume that the simulated sheared foam can be reduced to a single film whose constant thickness is fixed by the liquid fraction $\varphi$, giving $\tilde{W} = \varphi/2 = 0.15$. As detailed in SM [29], this choice of parameters only provides a poor or partial agreement between the simulations and the model. The discrepancy indicates that the complexity of the former is not reducible to the latter with plain geometric matching but that effective parameters are needed. Now, as illustrated in Fig. 9, taking $\tilde{\omega} = 3$ and $\tilde{W} = 0.5$ yields for the rescaled dissipations $\langle \tilde{D}_s \rangle$ an agreement that is at least qualitative. One may interpret the large value of $\tilde{W}$ as indicating that the complex geometry, with three films meeting at a junction cannot be assimilated to a unique film of similar thickness. Figure 9 shows $\langle \tilde{D}_s \rangle$ as a function of $h$ for various Pe and Bi, for the simulations and a mixed model, which takes a mean between symmetric and antisymmetric forcings. We observe that the Pe-independent plateaus at low $h$ match quantitatively, whereas the existence and location of the maximum agree at least qualitatively.

As a second comparison, we now identify which mechanism dominates the surfactant dissipation. As discussed in Sec. III C 2 above, the prediction from the extended Lucassen model, provided Pe is above unity, is remarkably simple and reads

$$ \frac{\langle D_{\text{dif}} \rangle}{\langle D_{\text{ex}} \rangle} \approx D_a, \qquad (22) $$

with $D_a = \bar{\chi} h \text{Bi} \sqrt{\text{Pe}/2}$ the Damköhler number, which compares sorption and diffusion timescales. Plotting the simulation data for $\langle D_{\text{dif}} \rangle / \langle D_{\text{ex}} \rangle$ versus $D_a$, we observe in Fig. 10 a good collapse of all curves, whatever the parameters considered. Moreover, for all low or intermediate $h$ values, the data and prediction are proportional to each other, with a prefactor $\approx 0.4$ of order unity, showing the ability of the simple model to capture a complex system. To explain the deviations visible at large $h$, we recall that the model assumes lateral invariance of surfactant concentration, which is not true in the simulations where the bubbles are deformed heterogeneously. This has two consequences. First, the contribution from surfactant diffusion along the interface is neglected. This is actually legitimate since this term is
<figure><img src="image_11.png" /><figcaption>FIG. 10. Ratio $\langle D_{\text{dif}} \rangle / \langle D_{\text{ex}} \rangle$ as a function of the Damköhler number $Da = hBi\sqrt{Pe}/2$ for the same simulations as in Fig. 4. The circles correspond to the three cases shown in Fig. 4c. The dashed line corresponds to $y = 0.4x$.</figcaption></figure>

numerically negligible [36]. Second, concentration gradients parallel to the interface are forbidden. This is different from simulations (cases 2-3 Fig. 10), pointing to a possible cause for discrepancy at large ratios values.

Overall, although its geometry is drastically simplified compared to a sheared foam, the model is able to capture some of the generic features of the dissipation mechanisms, provided we focus on time-averaged global quantities.

V. CONCLUSION

To summarize, we investigated numerically the surfactant dissipation in a sheared foam and found that a linear theory for the surfactant-induced dissipation can capture some of the trends observed in T1 simulation, as regards magnitude and origin. In both simulations and in a simplified model, we used a thermodynamic approach to identify which transport mechanisms –diffusion or exchange– dominates surface dissipation. Our findings may help to rationalize the influence of surfactants on foam behavior. Strong effects of surfactant physico-chemistry are reported [8, 20–22], when insoluble surfactants (e.g. dodecanol, fatty acids) are added to solutions of soluble ones (e.g. sodium dodecyl sulfate, SLES-CAPB mixtures). In terms of dimensionless numbers, the insoluble and soluble molecules are similar in size, hence similar in Pe values. Also, the adsorption lengths vary from 1 to 5 µm [18], leading to $h$ values with comparable order of magnitude. By contrast, exchanges can be slowed down by three orders of magnitude in the presence of insoluble species [37], corresponding to a decrease in Biot number by a factor $10^3$. Our results (Fig. 9) show that lowering Bi can increase the total surface dissipation. Such a trend is qualitatively consistent with the observed increase of foam viscosity [8] and slowing down of bubble rearrangements [20–22], and suggests that surfactant dissipation might be one key factor controlling the impact of surfactant composition on foam rheology.

This work should be pursued in several directions. First, it motivates further experimental work aimed at correlating specific microscopic surfactant properties and rheological measurements. Second, future efforts must be made towards extending our model to account for the viscous dissipation in the fluid, so as to better understand its coupling with surfactant properties and overall foam rheology. Indeed, viscous dissipation has been shown to play a key role when assemblies of large bubbles are deformed [18, 21], and is also affected by surfactant dynamics. Third, simulations of more realistic dry foams are as of yet prohibited by the issue of numerical coalescence [38]. New numerical developments are needed to handle this point. Finally, from a wider perspective, the methodology presented here is not specific to foams but is applicable to many surfactant-controlled systems, including emulsions [9], drop and bubble dynamics [39] and thin films [40].
*Acknowledgements* This research was funded by a grant from ANR (project Surfbreak, ANR-18-CE08-0013, 2018). We thank Isabelle Cantat and Christophe Ybert for inspiring discussions.

[1] Cantat, I., Cohen-Addad, S., Elias, F., Graner, F., Hohler, R., Pitois, O., Rouyer, F., Saint-James, A., *Foams - Structure and Dynamics*, (Oxford 2013).

[2] Hill, C., Eastoe, J., Foams: From nature to industry. *Adv. Colloid Interface Sci.* **247**, 469 (2017).

[3] Stevenson, P., editor, *Foam Engineering*, (Wiley-Blackwell 2012).

[4] Weitz, D., Foam flow by stick and slip. *Nature* **381**, 475 (1996).

[5] Höhler, R., Cohen-Addad, S., Rheology of liquid foam. *J. Phys.: Condens. Matter* **17**, R1041 (2005).

[6] Drenckhan, W., Hutzler, S., Structure and energy of liquid foams (2015).

[7] Bergeron, V., Forces and structure in thin liquid soap films. *J. Phys.: Condens. Matter* **11**, R215 (1999).

[8] Denkov, N. D., Tcholakova, S., Golemanov, K., Ananthpadmanabhan, K. P., Lips, A., The role of surfactant type and bubble surface mobility in foam rheology. *Soft Matter* **5**, 3389 (2009).

[9] Cohen-Addad, S., Höhler, R., Pitois, O., Flow in Foams and Flowing Foams. *Annu. Rev. Fluid Mech.* **45**, 241 (2013).

[10] Johnson, D. O., Stebe, K. J., Oscillating bubble tensiometry: A method for measuring the surfactant adsorptive-desorptive kinetics and the surface dilatational viscosity. *J. Colloid Interface Sci.* **168**, 21 (1994).

[11] Wantke, K.-D., Fruhner, H., The oscillating bubble method, in *Studies in Interface Science*, volume 6, 327–365, (Elsevier 1998).

[12] Langevin, D., Rheology of Adsorbed Surfactant Monolayers at Fluid Surfaces. *Annu. Rev. Fluid Mech.* **46**, 47 (2014).

[13] Christov, C., Ting, L., Wasan, D., The apparent dilational viscoelastic properties of fluid interfaces. *J. Colloid Interface Sci.* **85**, 363 (1982).

[14] Lucassen, J., Van Den Tempel, M., Dynamic measurements of dilational properties of a liquid interface. *Chem. Eng. Sci.* **27**, 1283 (1972).

[15] Besson, S., Debrégeas, G., Statics and dynamics of adhesion between two soap bubbles. *Eur. Phys. J. E* **24**, 109 (2007).

[16] Besson, S., Debrégeas, G., Cohen-Addad, S., Höhler, R., Dissipation in a sheared foam: From bubble adhesion to foam rheology. *Phys. Rev. Lett.* **101**, 214504 (2008).

[17] Bussonnière, A., Shabalina, E., Ah-Thon, X., Le Fur, M., Cantat, I., Dynamical Coupling between Connected Foam Films: Interface Transfer across the Menisci. *Phys. Rev. Lett.* **124**, 018001 (2020).

[18] Bussonnière, A., Cantat, I., Local origin of the visco-elasticity of a millimetric elementary foam. *J. Fluid Mech.* **922**, A25 (2021).

[19] Durand, M., Stone, H. A., Relaxation Time of the Topological T1 Process in a Two-Dimensional Foam. *Phys. Rev. Lett.* **97**, 226101 (2006).

[20] Biance, A., Cohen-Addad, S., Höhler, R., Topological transition dynamics in a strained bubble cluster. *Soft Matter* **5**, 4672 (2009).

[21] Petit, P., Seiwert, J., Cantat, I., Biance, A., On the generation of a foam film during a topological rearrangement. *J. Fluid Mech.* **763**, 286 (2015).

[22] Le Merrer, M., Cohen-Addad, S., Höhler, R., Duration of bubble rearrangements in a coarsening foam probed by time-resolved diffusing-wave spectroscopy: Impact of interfacial rigidity. *Phys. Rev. E* **88**, 022303 (2013).

[23] Manikantan, H., Squires, T. M., Surfactant dynamics: Hidden variables controlling fluid flows. *J. Fluid Mech.* **892**, 1 (2020).

[24] Though the numerical method of level-set is not strictly speaking a simulation, we use the term “simulation” for convenience and brevity.

[25] Erik Teigen, K., Song, P., Lowengrub, J., Voigt, A., A diffuse-interface method for two-phase flows with soluble surfactants. *J. Comput. Phys.* **230**, 375 (2011).

[26] Osher, S., Fedkiw, R., *Level Set Methods and Dynamic Implicit Surfaces*, volume 153, (Springer 2003).

[27] Sethian, J. A., Smereka, P., Level set methods for fluid interfaces. *Annu. Rev. Fluid Mech.* **35**, 341 (2003).

[28] Titta, A., Le Merrer, M., Detcheverry, F., Spelt, P., Biance, A.-L., Level-set simulations of a 2D topological rearrangement in a bubble assembly: effects of surfactant properties. *J. Fluid Mech.* **838**, 222 (2018).

[29] See Supplemental Material at ... for some details on the numerical simulation and full calculation of the extended Lucassen model. It includes Refs. [12,18,23,25,28,30].

[30] de Groot, S. R., Mazur, P., *Non-Equilibrium Thermodynamics*, (Dover, New York 1985).

[31] Buzza, D. M. A., Lu, C.-Y. D., Cates, M. E., Linear Shear Rheology of Incompressible Foams. *J. Physique II* **5**, 37 (1995).

[32] Levich, V. G., *Physicochemical Hydrodynamics*, (Prentice-Hall, Englewood Cliffs 1962).

[33] van den Tempel, M., Lucassen, J., Lucassen-Reynders, E. H., Application of surface thermodynamics to Gibbs elasticity. *J. Phys. Chem.* **69**, 1798 (1965).

[34] The Bi² term in the denominator of Eq.(15) is replaced with Bi²(1 + χh)².

[35] The explicit formula is $\Delta = |\bar{\zeta}^{-1} + 1 - i|^2 / (|\bar{\zeta}^{-1} + 1 - i|^2 + |\zeta^{-1} + 1 - i|^2)$.

[36] Because of the high surface Péclet number Pe$_f$ ≡ HU/D$_f$ = 10$^3$, this term contributes to less than 2% of the total surfactant dissipation in 21 out of 22 simulations. The only exception is a contribution of 10% when Pe = 100, Bi = 10, h = 1.
[37] Minkov, I. L., Arabadzhieva, D., Salama, I. E., Mileva, E., Slavchov, R. I., Barrier kinetics of adsorption-desorption of alcohol monolayers on water under constant surface tension. *Soft Matter* **15**, 1730 (2019).

[38] Innocenti, A., Jaccod, A., Popinet, S., Chibbaro, S., Direct numerical simulation of bubble-induced turbulence. *J. Fluid Mech.* **918**, A23 (2021).

[39] Anna, S. L., Droplets and Bubbles in Microfluidic Devices. *Annu. Rev. Fluid Mech.* **48**, 285 (2016).

[40] Ou Ramdane, O., Quéré, D., Thickening Factor in Marangoni Coating. *Langmuir* **3**, 2911 (1997).