Mean-field model for the bubble size distribution in coarsening wet foams

Abstract

Aqueous foams are subject to coarsening, whereby gas from the bubbles diffuses through the liquid phase. Gas is preferentially transported from small to large bubbles, resulting in a gradual decrease of the number of bubbles and an increase in the average bubble size. Coarsening foams are expected to approach a scaling state at late times in which their statistical properties are invariant. However, a model predicting the experimentally observed bubble-size distribution in the scaling state of foams with moderate liquid content, as a function of the liquid fraction ϕ, has not yet been developed. To this end, we propose a three-dimensional mean-field bubble growth law for foams without inter-bubble adhesion, validated against bubble-scale simulations, and use it to derive a prediction of the scaling-state bubble-size distribution for any ϕ from zero up to the unjamming transition ϕ_c ≈ 36%. We verify that the derived scaling state is approached from a variety of initial conditions using mean-field simulations implementing the proposed growth law. Comparing our predicted bubble-size distribution with previous simulations and experimental results, we likewise find a large population of small bubbles when ϕ > 0, but there are qualitative differences from prior results which we attribute to the absence of rattlers, i.e. bubbles not pressed into contact with their neighbours, in our model.

---

1 Introduction

Aqueous foams—packings of gas bubbles in liquid^1,2—are of theoretical importance as complex fluids whose properties arise from their well-characterised small-scale structure.³ Foams also have many industrial applications, including fire suppression, ore separation, foods, and drinks.^1,4,5 However, a foam's lifetime is limited by film rupture, gravitational drainage, and coarsening.² Our focus is on the last of these instabilities, which occurs due to solubility of the gas phase in the liquid: dissolved gas diffuses between bubbles, preferentially transferring from high-pressure bubbles to those at lower pressure.^6,7 Small bubbles tend to have higher pressure (the pressure of isolated spherical bubbles is inversely proportional to their radius by the Young–Laplace law¹), and thus shrink until they disappear, while large bubbles grow and the mean bubble size increases.^2,6 We neglect film rupture and drainage henceforth, noting that they can be mitigated in experiments, albeit with difficulty in the latter case.^8–10

The effects of coarsening are fairly well characterised in both the dry and dilute limits of a foam's liquid content,^2,11–17 which is measured by the liquid fraction ϕ; i.e. the ratio of liquid volume to the foam's total volume. At late times, the foam approaches a scaling state in which its statistical properties do not change with time, except for a scaling arising from growth of the average bubble size.^2,11,13,17 Let R denote the effective radius of a bubble; i.e. the radius of a sphere with the same volume (or that of a circle of the same area in 2D).¹ Then the growth exponent α with which the mean radius 〈R〉 (along with other averages¹⁸ of R) increases with time t in the scaling state, via 〈R〉 ∼ t^α, is 1/2 in the dry limit ϕ → 0 and 1/3 in the dilute limit ϕ → 1.^2,11,17 However, real foams lie between these limits, and moderately-wet foams arise in applications such as fire suppression and ore separation.^4,5

Coarsening at moderate ϕ is not yet well characterised. Experiments and simulations indicate that the foam still approaches a scaling state,^19–22 although they are not all consistent regarding the form of the crossover in α from 1/2 to 1/3 as ϕ increases,^8,18–21 and unexpected bubble-size distributions have been observed in simulations²¹ and in recent experiments on the International Space Station (ISS).²² The observed distributions exhibit a large population of small bubbles at moderate ϕ, which has been interpreted as resulting from the small shrinkage rate of rattlers (also termed roamers); i.e. small bubbles, which have been observed directly, that are either out of contact with their neighbours, or that have contacts due only to inter-bubble adhesion.^21–23 However, rigorous theoretical predictions of the variation of α with ϕ, and of the scaling-state bubble-size distribution, have yet to be developed.

Progress on predicting α has recently been made by Durian,²⁴ who derived an approximate growth law for the average bubble radius, with contributions from gas transfer through thin films separating the bubbles, and through the bulk liquid in the foam. The growth exponent α was taken to be a weighted average of these contributions, with the typical film area obtained by fitting to experimental data.²⁴ Encouraging agreement with experiments was found,^8,18,24 but there is considerable scope for improving the argument's rigour, and possibly for avoiding the use of free parameters.

Our focus in the present work is instead to derive an approximation of the scaling-state distribution as a function of ϕ. We apply standard techniques^11,25 to derive the distribution predicted by the three-dimensional (3D) version of a mean-field bubble growth law we proposed in prior work.²⁶ This is a border-blocking growth law; i.e. it accounts only for gas transfer through a bubble's contact films, omitting that through the bulk liquid,^9,27 and hence corresponds to the limit in which the ratio of film thickness h to bubble size tends to zero.^7,18,24 This would appear to be a natural first case to consider: as a foam coarsens, 〈R〉 grows while h is usually expected to remain approximately constant,^1,9,18 and so the limit h/〈R〉 → 0 should be approached at late times. A border-blocking model might therefore be expected to describe the eventual coarsening dynamics of a real foam, and thus its scaling state (analogous behaviour is understood to occur in alloy coarsening, where the dynamics is eventually dominated by bulk diffusion instead^28,29). This argument assumes that there are films in the foam,† and thus that the foam is either jammed, or flocculated due to bubble adhesion.^18,31,32 For simplicity, and as the natural starting point for developing theory, we consider only foams with zero adhesion in the present work; i.e. with a contact angle θ = 0 between film and bulk-liquid interfaces.³³ Therefore, we restrict our attention to jammed foams, with ϕ < ϕ_c, where ϕ_c ≈ 36% is the unjamming transition in 3D.²

Mean-field models are widely used to provide tractable approximations of coarsening systems.^6,18,34,35 In the dry limit ϕ = 0, Lemlich's law^6,12,36 approximates the growth rate of a 3D bubble as follows. Let the bubble have effective radius R, and let γ be the liquid/gas surface tension, D a gas diffusion coefficient, H Henry's constant² (for brevity, we use the definition employed by Schimming and Durian,⁷ which incorporates the molar gas volume²), and h the film thickness. Let R_c denote the critical bubble radius, which is the radius of bubbles which neither grow nor shrink under coarsening,¹¹ and which is set by the condition that the total gas volume is conserved.⁶ Under the assumption that the equivalent volume of dissolved gas²⁵ is negligible compared to the bubble volumes, this condition requires that the instantaneous mean volumetric bubble growth rate satisfies^2,6 〈V̇〉 ≡ 〈dV/dt〉 = 0. Lemlich's law is then^6,12,36

where^2,6R_c = 〈R²〉/〈R〉 ≡ R₂₁. This law is approximate due to mean-field assumptions made during its derivation, and we consider it as a prediction for dry foams (ϕ = 0) due to its assumption that the whole surface of a bubble is available for gas transfer through films.^6,35 In the dilute limit ϕ → 1, for which all gas flow is through the bulk liquid (Ostwald ripening),² the Lifshitz–Slyozov–Wagner (LSW) law gives^11,12in 3D, where imposing 〈V̇〉 = 4π〈R²Ṙ〉 = 0 gives¹²R_c = 〈R〉.

Much prior work^29,34,37,38 has been done to develop mean-field growth laws for intermediate values of ϕ. However, to our knowledge, most of these studies are focussed on more general coarsening systems, rather than foams specifically, and thus do not typically model the films between contacting bubbles. Gas transfer in foams is most efficient through films due to their small thickness,² and so it seems necessary to approximate film sizes to obtain an effective foam coarsening model. Progress has been made by Pasquet et al.,¹⁸ who obtain a correction to the prefactor in eqn (1) for ϕ > 0 using an approximation for film area as a function of ϕ derived by Höhler et al.,³⁹ as we shall discuss in Section 2. While their model predicts the variation of the coarsening rate with ϕ, it is not able to predict the effect of ϕ on the bubble size distribution, since their growth law is a constant multiple of Lemlich's law (1) at fixed ϕ (see Section 2), thus justifying development of the more detailed model that we define below.

In Section 2, we define our proposed growth law, and compare its predictions with our previous bubble-scale simulations.^26,40 Then, in Section 3, we adapt standard techniques^11,25 to derive the scaling-state bubble-size distribution predicted by this growth law. Supporting derivations and properties of the distribution are given in Appendix A. Next, in Section 4, we give the results of mean-field simulations using our proposed growth law to check that the derived scaling state is approached from various initial conditions. Our simulation methods are adapted from those of De Smet et al.⁴¹ We then compare the derived distribution with previous experiments^22,23 and simulations²¹ in Section 5. We conclude in Section 6, and discuss simulation convergence and apparent exceptions to the universality of our derived scaling state (comparable to apparent exceptions previously observed for other mean-field laws^29,42,43) in Appendix B.

2 Mean-field growth law for wet foams

The border-blocking growth rate of a bubble in a 3D wet foam is known exactly under standard assumptions (including that the thin film thickness h is the same for all contact films), although we will make mean-field approximations below in order to develop a tractable model. Let the considered bubble have effective radius R (i.e. the radius of a sphere with equal volume¹), pressure p (we define all pressures relative to that of the liquid, which is constant), and n contacting neighbours (and hence n films). Let the bubble's kth contact film have area A_k and adjoin a bubble of pressure p_k. Then the bubble's border-blocking growth rate is exactly²

In two dimensions (2D), the equivalent growth law is²⁷

where L_k is the length of the film between the bubble and its kth neighbour. Due to the dependence of eqn (3) on multiple bubble properties (not just R, for example), its implications for the evolution of a coarsening foam as a whole are not clear, inviting further approximations. Pasquet et al.¹⁸ proposed a mean-field border-blocking growth law which generalises eqn (1)⁶ to wet foams. It is based on the following approximation, which was derived by Höhler et al.,³⁹ for the ratio of total contact film area A to surface area S of bubbles in monodisperse 3D foams. Let the foam's osmotic pressure^33,44 be Π_O, and define the scaled osmotic pressure [capital Pi, Greek, circumflex] = Π_O(R₃₂/γ)/(1 − ϕ) (this definition differs from that typically used in other studies^39,44 by the factor of 1/(1 − ϕ), which we incorporate for brevity in the results given below), where²R₃₂ ≡ 〈R³〉/〈R²〉. Then³⁹

Following Pasquet et al.,¹⁸ let us adapt the derivation⁶ of Lemlich's law (1). We simplify eqn (3) by making the mean-field approximation that all neighbours and films of the considered bubble are identical, and take the bubble to have pressure p = 2γ/R; i.e. that of an isolated spherical bubble,⁶ by the Young–Laplace law.¹ The equal neighbour pressures are written p_k = 2γ/R_c (anticipating that R_c will serve as the critical radius), and the bubble's total film area is approximated using eqn (5) (taking S ≈ 4πR², noting that the error in this approximation is about 10% for a dry 3D foam⁴⁵). Pasquet et al.¹⁸ thereby obtained an approximate growth law of the form

where R_c = R₂₁ since the law has the same functional form as eqn (1). While this growth law predicts the rate at which the foam coarsens¹⁸ (see Section 4.2 below), it is not suitable for approximating the variation of the scaling-state bubble-size distribution with ϕ—the right-most factor, approximating the bubble's film coverage, scales all bubble growth rates equally, and so the same scaling-state distribution, that of Lemlich's law (1),^12,36,46 is obtained for all ϕ < ϕ_c (beyond which [capital Pi, Greek, circumflex] = 0 and there is no coarsening in the border-blocking model^7,39). This suggests the development of growth laws which account for the correlation between bubble pressure and film area, an approach which was proposed by Pasquet et al.¹⁸

We now describe a 2D mean-field growth law that we derived previously by following this approach,²⁶ before introducing the corresponding growth law obtained by adapting the derivation to 3D. We validate both growth laws against our bubble-scale simulations.^26,40

2.1 Two-dimensional mean-field growth law

In ref. 26, we developed mean-field approximations for bubble pressure p and the ratio of a bubble's total film length L_f to its perimeter P, which depend only on bubble radius R (for the non-adhesive foams we consider here). Thus, the correlation between p and L_f/P is approximated through their mutual dependence on R. These approximations for p and L_f/P are generalisations of two results from Höhler et al.³⁹—a relation between the foam's capillary and osmotic pressures, and eqn (5)—to individual bubbles in polydisperse foams, and are derived using similar techniques (but are subject to additional simplifying assumptions). In ref. 26, we used the pressure and film length approximations, in a similar manner to the derivation⁶ of eqn (1) and (6), to derive the following approximate mean-field border-blocking growth law, which incorporates the variation of film coverage with bubble size. Let [capital Pi, Greek, macron] = Π_O(R₂₁/γ)/(1 − ϕ) be a 2D equivalent of [capital Pi, Greek, circumflex], and let R̄ = R/R₂₁. Then²⁶where R_c is the critical radius.

In ref. 26, we compared eqn (7) with the growth rates predicted by bubble-scale finite-element simulations of 2D foams. These simulations were implemented in the Surface Evolver,^47,48 and adapted the approaches of Kähärä et al.⁴⁹ and Boromand et al.:⁵⁰ every liquid/gas interface in the foam was resolved (requiring a film thickness considerably larger than expected in real static foams), and made to interact with other interfaces through a disjoining pressure. Hence, the simulated bubbles were deformable. We found that eqn (7) appears to approximate the trend in the simulated border-blocking growth rates (when corrected for variations in simulated bubble film thickness h; see below), albeit with a large degree of scatter among individual bubbles.

Using the data⁵¹ from ref. 26, we see in Fig. 1 that eqn (7) is actually in good agreement with the binned mean growth rate. Here, we have used eqn (4) to calculate the border-blocking growth rates of the simulated bubbles (expressed in terms of their effective neighbour number^19,26). Eqn (4) assumes that all films have equal thickness h, which is a standard approximation for real foams,^1,9,18 whereas film thickness varies considerably between bubbles in our simulations for numerical reasons, an effect which would tend to enhance the shrinkage rates of small bubbles (see ref. 26 for details). However, we do not expect the thickness variations to affect the bubble pressures and film sizes substantially, and so our use of eqn (4) should correct for these variations.

Fig. 1 Border-blocking growth rate (expressed as the rate of change in bubble area A_b = πR²) versus effective radius for 1024 bubbles in 2D simulated foams at (a) ϕ = 2% and (b) ϕ = 10% without inter-bubble adhesion, the simulation data having been taken from our previous study.^26,51 The bubble-area distribution is a compressed exponential fitted to experimental data by Roth et al.⁹ The individual bubble data is shown alongside its binned mean (20 equal bins are used; the error bars give the standard deviation within each bin). Comparison is made with eqn (7) (—) and the 2D version of eqn (6) (⋯), in which R_c is obtained by numerically solving 〈dA_b/dt〉 = 0 as described in the text and [capital Pi, Greek, macron] is measured in the simulations.²⁶ The values of R_c are shown on the panels. The calculation of dA_b/dt for the simulated bubbles is explained in the text.

We calculate R_c in Fig. 1 by summing the bubble area growth rates dA_b/dt = 2πR dR/dt predicted by the considered mean-field growth law (eqn (7) or the 2D equivalent of eqn (6)) for each bubble in the finite-element simulations. We then use a numerical rootfinder⁵² to obtain R_c such that this sum is zero. Hence, we apply the 2D equivalent 〈dA_b/dt〉 = 0 of the gas conservation condition^2,6 introduced in Section 1 to the simulated foams. The osmotic pressure is measured directly in the simulations by evaluating its definition in terms of the bubble interface geometry³⁹ (see ref. 26 for details).

In Fig. 1, we also see improved agreement over the 2D equivalent^26,36 of eqn (6) for small bubbles, further justifying the more complicated growth law (7).

2.2 Three-dimensional mean-field growth law

A 3D equivalent of eqn (7) can be derived using very similar arguments, with the required approximations for bubble pressure and film area being obtained by slightly adapting the 2D derivations.²⁶ The arguments in the 3D case are given in ref. 40, and are summarised within Appendix A. The resulting mean-field approximation for an individual bubble's pressure p (relative to that of the liquid; i.e. the bubble's capillary pressure³⁹) is

The corresponding approximation for the ratio of a bubble's total film area A to its surface area S (i.e. the fraction of its surface in contact with other bubbles) is

where R̂ = R/R₃₂. Eqn (9) predicts that smaller bubbles are less deformed by their contacting neighbours^27,53 (i.e. they have smaller A/S), and that A/S is larger in drier foams (i.e. for larger [capital Pi, Greek, circumflex]) as expected.³⁹ The scaling of radius by R₃₂ arises from an approximation that all neighbouring bubbles have pressure equal to the capillary pressure Π_c of the foam³⁹ (see Appendix A). Other pressures could have been chosen, but we shall see below that the resulting growth law is insensitive to a replacement of R₃₂ by R_c, which is another average radius. Eqn (9) becomes eqn (5)³⁹ in a monodisperse foam (where R̂ = 1 for all bubbles).

An approximate 3D growth law is then obtained, like eqn (1) and (6), by assuming identical films and neighbours in eqn (3). Eqn (8) is used for the bubble's pressure (instead of the Young–Laplace law for spherical bubbles), and the equal pressures of its neighbours (taken to have radius R_c), while eqn (9) gives the total film area A (again approximating⁴⁵S ≈ 4πR²). Hence, we obtain the approximate 3D mean-field border-blocking growth law

where R_c is the critical radius. Comparing with eqn (6), we see that there is an additional dependence on bubble size R, which is intended to approximate the variation of film size with R (see below), and thus more accurately describe a bubble's growth rate (still within a mean-field model). As noted above, this additional complexity is required to predict the variation of scaling-state bubble-size distribution with ϕ (upon which it is known to depend^22,54).

We provide a test of eqn (10) in Fig. 2, similar to Fig. 1, using 3D versions of the finite-element simulations discussed in Section 2.1. These 3D simulations are described in detail in ref. 40, and use the same approach as in 2D, that of Kähärä et al.⁴⁹ and Boromand et al.,⁵⁰ and thus are similar to prior simulations of deformable 3D frictional particles⁵⁵ and biological cells.⁵⁶ The simulations are again performed with the Surface Evolver,⁴⁸ and the initial dry-foam geometries are generated using the Neper software.^57,58 The border-blocking growth rates are calculated using eqn (3) (thus correcting for film thickness variations between bubbles as in the 2D simulations). The critical radius R_c in each of the plotted mean-field growth laws is separately obtained by numerically solving the gas conservation condition^2,6 similarly to Fig. 1: the considered growth law's prediction of dV/dt = 4πR²dR/dt for each bubble in the simulations is summed, and a rootfinder⁵² is used to give R_c such that 〈dV/dt〉 = 0 for the growth law. As in Fig. 1, the osmotic pressure is measured in the simulations using its definition in terms of the bubble interface geometry³⁹ (see ref. 40 for details).

Fig. 2 Border-blocking relative growth rate² (where V = 4πR³/3 is the bubble volume) versus effective radius for 256 bubbles (aggregated from four 64-bubble simulations at equal ϕ) in 3D foams, using finite-element simulations described in ref. 40. There is no inter-bubble adhesion, and the liquid fraction is (a) ϕ = 10% and (b) ϕ = 30%. The bubble-size distribution¹⁷ is lognormal with standard deviation 0.4 with respect to R/〈R〉. The simulation data is plotted as in Fig. 1, except 5 equal bins are used due to the smaller number of bubbles, and is compared with eqn (10) (—) and (6)¹⁸ (⋯); where [capital Pi, Greek, circumflex] and R₃₂ are the mean values measured in the aggregated simulations, and R_c is obtained⁶ by numerically solving 〈dV/dt〉 = 0 as explained in the text (the values of R_c are shown on the panels). The means of calculating dV/dt for the simulated bubbles is also given in the text. Comparison is additionally made to eqn (10) with R₃₂ replaced by R_c throughout (- -), as described later in the text.

Computational resources limit our 3D simulations to much smaller foam sizes than in 2D: in Fig. 2, we aggregate bubble data from four distinct 64-bubble foams with periodic boundary conditions. The individual plotted bubbles have their radii scaled by R₃₂ from the foam sample in which they originate (rather than the mean value of R₃₂ used in the plotted growth laws). While comparison to larger foams would be highly desirable to reduce finite-size effects and statistical fluctuations, Fig. 2 does indicate fairly good agreement in the mean between the simulations and eqn (10), and gives some evidence to prefer the latter growth law over eqn (6), at least at higher ϕ. However, we note that the level of agreement with eqn (10) for small bubbles does worsen if the number of bins is doubled to 10 (although the bin of smallest radius then contains few bubbles).

Fig. 2 shows that the main difference between our proposed growth law (10) and eqn (6)¹⁸ is that the former predicts suppressed shrinkage rates for small bubbles^2,59 (we note that the radius corresponding to the local minimum predicted by eqn (10) depends on the measure of growth rate selected; there is no minimum for dR/dt). This suppression follows from eqn (9), which implies that the fraction of a bubble's surface in contact with other bubbles decreases with decreasing bubble radius.²⁷ Thus, smaller bubbles have relatively-less film area through which gas may transfer in the border-blocking model we use.^9,27 On the other hand, eqn (6) assumes that the relative film areas are independent of bubble size, through its use of eqn (5).³⁹

As discussed later in Section 5, the scaling-state distribution is expected to be sensitive to the form of the growth law for small bubbles.^21,22 In 2D, Fig. 1 suggests that eqn (7) applies over a wide range of bubble sizes, although it would be desirable to test the growth law against simulations of more polydisperse foams, with more small bubbles. However, Fig. 2 suggests that, in 3D, eqn (10) may overestimate the shrinkage rate of small bubbles at higher ϕ, perhaps due to the presence of rattlers; i.e. bubbles without contact films.^21,22 Further 3D simulations of larger systems may resolve this—the discrepancy might be an artefact of aggregating data from several foams. The role of rattlers in mean-field models remains unclear, and will be discussed in Section 5. In spite of these unresolved matters, the level of agreement in the mean observed in Fig. 1 and 2 (particularly the former, noting that the same arguments are used in 2D and 3D to derive the growth law) encourages us to explore the consequences of eqn (10) in 3D, which is the case of primary practical interest. The 2D growth law (7) is used no further in the present work.

2.3 Application and simplification of the 3D growth law

We derive the scaling state of eqn (10) in Section 3. First, however, we note that a foam is more commonly parametrised by liquid fraction ϕ than by osmotic pressure Π_O (which appears in scaled form [capital Pi, Greek, circumflex] within the growth law). In a foam without bubble adhesion (such as we consider here), Π_O is expected^39,44 to have a one-to-one correspondence with ϕ up to the unjamming transition ϕ_c, beyond which Π_O = 0. We convert between ϕ and Π_O using the empirical law^39,60for critical liquid fraction ϕ_c = 36% and K = 3.2.¹⁸ The extra factor of 1 − ϕ compared to Maestro et al.⁶⁰ is due to our different definition of the scaled osmotic pressure [capital Pi, Greek, circumflex] (above). In real foams, the critical liquid fraction depends on the bubble-size distribution, and is likely to be smaller in the polydisperse scaling-state distributions we derive below.^21,22 For example, Galvani et al.²² found ϕ_c = 31% for a simulated foam with polydispersity equal to the foams in the ISS experiments. However, for concreteness, we select the standard value ϕ_c = 36% which applies to monodisperse disordered 3D foams.² We emphasise that eqn (11) was not used in Fig. 2, where, as in Fig. 1, the osmotic pressure was measured directly from the finite-element simulations.^2,26,40

We also simplify the growth law (10) by replacing all instances of R₃₂ therein by the critical radius R_c. Both measures of the average bubble size are expected to be fairly close in value. This approximation affects the derived scaling-state distributions very slightly (not shown), but gives the great advantage that the scaling-state properties can be derived analytically in closed form. This replacement is achieved by redefining R̂ = R/R_c and [capital Pi, Greek, circumflex] = Π_O(R_c/γ)/(1 − ϕ), and its very small effect on the predicted growth rates is shown in Fig. 2. These redefinitions are retained for the remainder of this article. We note that a similar replacement was made by Ardell³⁷ in a generalisation of the LSW law (2) to ϕ < 1.

3 Derivation of scaling state

We now derive ρ(R̂), the scaling-state bubble-size probability distribution predicted by the 3D border-blocking growth law (10). We follow the approach of Marqusee and Ross,²⁵ which is among the methods commonly used to derive scaling-state distributions in Ostwald ripening and coarsening theory.^{11,36,38,61,62} The derivation will also give us the growth exponent α (see Section 1).

Using a similar notation to Marqusee and Ross,²⁵ let r = R/L and τ = t/T be dimensionless bubble radius and time variables, with T = L²h/(2γDH) selected to simplify eqn (10), and with L chosen such that the total (conserved) gas volume is unity.‡ Let the dimensionless critical radius be r_c(τ) = R_c(τ)/L, which is expected to increase with time like the average bubble size³⁵ (with the same growth exponent¹⁸α) and define the dimensionless bubble radius r̂ = r/r_c (= R̂) for consistency of notation. The growth law, eqn (10), becomes

Let n(r,τ) be the bubble number distribution, such that n(r,τ)dr is the number of bubbles with radius between r and r + dr at time τ, and is the total number of bubbles. Since bubble radii vary continuously, the distribution n satisfies the continuity equation¹¹

Following Binder,⁶¹ let x = rτ^−α be the scaled bubble size, recalling that α is the growth exponent (and so fixed x corresponds to fixed R/R_c, for example). Define^11,25n̄(x,τ) as the number distribution in the variables (x,τ). Since n̄dx = ndr, we have n̄ = τ^αn. Substituting this definition into eqn (13), and expressing in terms of the derivatives of the variables (x,τ) (noting that this change of variables alters ∂/∂τ) gives the continuity equation^11,25 for n̄,

where²⁵is the velocity of a bubble's scaled radius x (i.e. dx/dτ).¹¹ We seek a scaling-state solution n̄(x,τ) of eqn (14). Recalling Section 1, this is a solution for which the distribution of bubble sizes is independent of time except for an overall scaling,² which is incorporated into the definition of x. Since bubbles shrink and disappear during coarsening,¹ the bubble number density n̄ will decrease with time while the x distribution remains constant. Hence, a scaling state takes the separable form^25,63n̄(x,τ) = [small rho, Greek, macron](x)[small sigma, Greek, macron](τ), where [small rho, Greek, macron] is proportional to the probability distribution in x (which is τ-independent in the scaling state), and [small sigma, Greek, macron] accounts for gradual disappearance of bubbles, and thus the decay in the number density n̄.

Following Binder,⁶¹ we assume [small sigma, Greek, macron] = τ^−β, and obtain the exponent β using conservation of the amount of gas. Marqusee and Ross²⁵ account for the amount of dissolved discontinuous phase when imposing its conservation. However, as in Section 1, we take the total bubble volume to account for all gas in the system (which is conserved),⁶¹ thus presuming a sufficiently low gas solubility in the foaming liquid. We recall that r is defined such that the total conserved gas volume is unity after nondimensionalisation; i.e.

Substituting r = xτ^α, and using the fact that n = τ^−αn̄ = τ^−α−β[small rho, Greek, macron](x) in the scaling state, gives

The right-hand side is evidently independent⁶¹ of time τ, and so β = 3α and n̄ = τ^−3α[small rho, Greek, macron](x) in the scaling state. Substituting this form for n̄ into eqn (14), we obtain²⁵

The exponent α is determined by noting that the right-hand side of this equation can have no explicit dependence^25,61 on time τ (since the left-hand side has none), once eqn (12) and (15) are substituted for u. We define x_c = r_cτ^−α, which is τ-independent in the scaling state since α is the growth exponent. Hence, x̂ = x/x_c (= r̂) also has no explicit dependence on τ in the scaling state. Furthermore, [capital Pi, Greek, circumflex] is τ-independent since we fix the liquid fraction ϕ. By eqn (12) and (15),

Substituting^25,61 this result into eqn (18) gives α = 1/2, as expected in experiments using dry foams,^14,17 and as in the small-ϕ model of Schimming and Durian.⁷ However, eqn (10) is intended as an approximate growth law for all ϕ at which bubbles are in contact, and for which h is vanishingly small, so that gas transfer via films dominates.^7,18 Therefore eqn (10) accords with the standard expectation that α = 1/2 in any such border-blocking case.¹⁸ Substituting α = 1/2 into eqn (18), and using that τu is then τ-independent, we obtain a linear, first-order ordinary differential equation for [small rho, Greek, macron], which is straightforwardly solved (by separation of variables) to give⁶¹

where C is a constant. Rewriting eqn (19), for α = 1/2, we obtain

The numerator is quadratic in x, so τu has at most two zeros. We now follow standard arguments^11,25,36,62 to constrain these zeros, and thus to determine x_c analytically.§ First, we argue that [small rho, Greek, macron](x) = 0 for all x ≥ x₀, where x₀ > 0 is a finite cutoff.^11,25 Following Marqusee and Ross,²⁵ we note that, at large x, eqn (18) becomes

on substituting α = 1/2 and eqn (19). Thus [small rho, Greek, macron] = k/x⁴ for constant k. But k = 0 is the only value for which the integral in eqn (17) is finite.²⁵ Therefore, [small rho, Greek, macron] must be zero beyond some cutoff,^11,25 ¶ denoted x₀ (which we define to take its smallest possible value).

Let us consider a bubble b₀ which initially has scaled radius x = x₀ in the scaling state, and which therefore is no smaller than any of the other bubbles. We recall¹¹ that u gives the rate of change of x for b₀. Thus, if u < 0 at x₀, then x will decrease. But u is a single-valued function of x (and τ) for all bubbles, by eqn (19), and so all bubbles initially with x smaller than x₀ will remain smaller than b₀ as it shrinks.¹¹ Therefore, [small rho, Greek, macron] will evolve,³⁶ with its cutoff equal to the decreasing scaled radius of b₀, contradicting the assumption that [small rho, Greek, macron] is a function of x only (and thus describes a scaling state). It follows that u ≥ 0 at x = x₀.

But u should be nowhere positive for x ≤ x₀ in the scaling state, because that would result in a class of bubbles which never vanish during coarsening,³⁶ contradicting the expected dynamics that all but one bubble should eventually vanish in a finite foam sample.^1,2 Therefore, we must have²⁵u = 0 at x = x₀; i.e. the scaled radius of the bubble b₀ is constant.

The quadratic numerator in eqn (21) thus has a zero at x = x₀. From this equation, u = −[capital Pi, Greek, circumflex]/(τx_c) at x = 0, and u ∼ −x/(2τ) as x → ∞. We also recall that [capital Pi, Greek, circumflex] > 0 in the jammed foams without bubble adhesion^33,44 that we consider. Therefore, if u has exactly a pair of distinct zeros, then both zeros are at x > 0, and u > 0 between these zeros. Hence, x₀ must equal the smaller of these zeros in such a case, since we have argued above that u must be nonpositive for all x ≤ x₀. But a stability problem then arises: if bubbles with x just below the cutoff x₀ are perturbed to lie beyond x₀, then they will continue to grow²⁵ (since u > 0), which will cause [small rho, Greek, macron] to evolve. This is not a contradiction (since we have perturbed the distribution to begin with), but suggests that the scaling state is unstable in this case^25,29,43 (see Appendix B for further discussion).

Therefore, by the above exclusion of other possibilities, the desired stable scaling state should be such that u has exactly one zero^11,25 at x = x₀; i.e. the quadratic in the numerator of eqn (21) has a double zero at x₀. Hence, the quadratic's discriminant gives

This quadratic equation can be solved for x_c, and the cutoff x₀ is then the double zero of the quadratic numerator in eqn (21). Therefore,||

In the dry limit as [capital Pi, Greek, circumflex] → ∞, for which eqn (10) reduces to Lemlich's law (1),^6,12,36 we have and , in agreement with Marqusee and Ross.²⁵ Furthermore, recalling that x₀ is the double zero of the numerator in eqn (21), we may rewrite this equation as

To obtain the scaling-state distribution, we substitute this into eqn (20), giving

for x < x₀. However, we desire the probability distribution ρ(R̂) of relative bubble size in the scaling state. We recall that R̂ ≡ R/R_c = x/x_c by definition (which has no explicit τ-dependence in the scaling state), while ρ(R̂)dR̂ = [small rho, Greek, macron](x)dx. Thus, ρ(R̂) = x_c[small rho, Greek, macron](x). We also define the dimensional cutoff radius R₀ = x₀Lτ^α with α = 1/2, analogous to R and R_c, and its relative value R̂₀ = R₀/R_c (= x₀/x_c). Therefore, eqn (26) givesas the probability distribution of relative bubble size in the scaling state, where C is chosen to normalise the distribution. Eqn (27) is valid for R̂ < R̂₀, whereas ρ(R̂) = 0 otherwise. From eqn (24), the relative cutoff is

In the dry limit [capital Pi, Greek, circumflex] → ∞, eqn (27) becomes the Lemlich–Markworth distribution^12,36,46

which is the scaling-state distribution of Lemlich's law (1), to which eqn (10) reduces in this limit. In the limit [capital Pi, Greek, circumflex] → 0, corresponding to the approach to the unjamming transition, eqn (27) becomes

We recall that gas transfer through the bulk liquid is neglected in our model, because we take h → 0. Hence, eqn (30) differs from the scaling-state distribution¹¹ predicted by the LSW law (2), for which all gas diffusion is via bulk liquid rather than films.²

We plot eqn (27) in Fig. 3 for a variety of liquid fractions ϕ, which determine [capital Pi, Greek, circumflex] viaeqn (11).⁶⁰ We also include the above limiting distributions. The normalisation constant C and the first few moments of ρ(R̂) can be obtained analytically, and we give these in Appendix A.

Fig. 3 Probability distribution ρ of the relative bubble size R̂ in the scaling state, predicted by eqn (27) (C is obtained here by numerical integration of the distribution), for various liquid fractions. We also plot ρ in the limiting cases of a dry foam and a foam at the unjamming transition, given respectively by the Lemlich–Markworth distribution (29)^12,46 and eqn (30).

Perhaps the most noticeable property of the probability distribution (27) is that, for all wet foams, we have ρ > 0 at R = 0, which is a qualitative difference from its form in the dry limit (29).^36,46 This may be interpreted using the approximation (9) for the relative film area A/S of bubbles: smaller bubbles have a smaller proportion of their surface in contact with neighbouring bubbles,²⁷ and thus relatively less area available for gas transfer under the assumption of border blocking.^9,27 In our mean-field model, A/S tends to zero as R → 0; i.e. bubbles of zero size are rattlers. Therefore, in spite of the divergence of bubble pressure as R → 0, a large population of small bubbles remains, in the scaling state, due to their disproportionately small films.** In the dry limit,⁶A/S = 1 independent of R, and so this effect does not arise. We discuss the predicted probability distribution (27) further in Section 5, where we compare it to prior simulations and experiments.^21–23 For now, we note that a large population of small bubbles was also found in these studies, but the form of the bubble size distribution is qualitatively different.^21–23,54

In Fig. 4, we give the variation with ϕ of R_c/R₂₁, R_c/R₃₂, R₀/R₂₁, polydispersity⁴⁵ [scr P, script letter P] = R₃₂/〈R³〉^1/3 − 1, the standard deviation σ_R of R/〈R〉, and geometric disorder² in the scaling state (where V = 4πR³/3 is the bubble volume). We recall from Section 1 that R_c = R₂₁ in the dry limit.⁶ All quantities plotted in Fig. 4 are measures of the distribution's width, and so they each increase with ϕ, consistent with the broadening of the distribution apparent from Fig. 3. In Fig. 4, R_c/R₃₂ ≈ 1, which suggests that our replacement of R₃₂ by R_c in eqn (10) (see Section 2) is self consistent. The poor agreement of the polydispersity in Fig. 4 with the ISS experiments (for which [scr P, script letter P] > 0.3 for ϕ ≲ 30%, and [scr P, script letter P] tends to decrease as ϕ increases²²) is interpreted in Section 5.

Fig. 4 Various properties of eqn (27) as functions of liquid fraction. The ratios of critical radius and cutoff radius to R₂₁ (and R₃₂ in the former case) are shown, alongside the polydispersity⁴⁵ [scr P, script letter P] and geometric disorder²σ_V, which are defined in the text, and the standard deviation σ_R of R/〈R〉. These are all evaluated numerically from eqn (27) and (28).

We recall, by the definition of x_c, that the critical radius satisfies r_c = x_cτ^α in the scaling state, where α = 1/2. But x_c is a constant, and therefore so is dr_c²/dτ = x_c², which provides a measure of the coarsening rate of the foam^18,64 (proportional to the rate at which a typical bubble's surface area increases). We plot this in Fig. 5 as a function of ϕ, and compare it with the corresponding rate predicted by Pasquet et al.¹⁸ using the simpler growth law (6).

Fig. 5 Coarsening rate (see text) versus liquid fraction ϕ, as predicted by the growth law (10)viaeqn (24), alongside the prediction of Pasquet et al.¹⁸ The predictions agree at ϕ = 0, where they equal 1/2, and they both reach zero at ϕ = ϕ_c. The variation of ρ(0) with ϕ, given by eqn (27) (with C again being obtained numerically), is also plotted (ρ(0) = 0 at ϕ = 0).

As expected, the coarsening rate decreases with increasing ϕ in Fig. 5 due to the decrease in film areas, reaching zero when the bubbles lose contact at the unjamming transition. Our prediction is close to that of Pasquet et al.¹⁸ (agreeing at ϕ = 0 since the two growth laws are equivalent in this case), although our rates are slightly higher for intermediate ϕ. This is perhaps counterintuitive, as our model has enhanced border-blocking for small bubbles compared to eqn (6), thus slowing their disappearance along, presumably, with the growth of the average bubble size. However, as expressed in more detail in Section 5, the rate of change of the bubble number dN/dt (and thus the coarsening rate) depends on ρ(0) in addition to the limiting growth rate as R → 0.^13,65

We note that the prediction¹⁸x_c² = [capital Pi, Greek, circumflex]/(4 + 2[capital Pi, Greek, circumflex]) plotted in Fig. 5, which originates from the film area approximation (5),³⁹ has been successfully fitted²³ to data from coarsening experiments at small ϕ ≤ 8%. However, its agreement with the ISS experiments at larger ϕ is poor, likely due to bubble adhesion.^18,64 Recalling that we consider only nonadhesive foams in the present work, we consider this discrepancy no further, except to note that an adaptation of eqn (5) to adhesive foams has been proposed and used by Galvani et al.²³

Having derived and analysed the scaling-state bubble size distribution (27) predicted by the border-blocking mean-field growth law (10), we next give results from mean-field simulations to verify that this distribution is approached from a range of initial conditions.^25,29

4 Simulations

4.1 Simulation methods

We now describe the approach we use to simulate foams that coarsen according to a mean-field bubble growth law of the LSW^11,12 or Lemlich^6,12,36 type; i.e. which depends only on bubble radius R, like eqn (10).

The numerical methods are adapted from those of De Smet et al.,⁴¹ and hence involve predicting the evolution of a large number of individual bubbles (rather than calculating the evolution of a distribution function^6,29,41). We select this approach for its directness, and we note that present computational resources allow considerably larger system sizes to be studied than were simulated by De Smet et al.⁴¹ To our understanding, a similar approach is used by Brown.⁴² We implement the simulations in Python, using the NumPy and SciPy packages.^52,66

We use the dimensionless bubble radius and time variables r and τ defined in Section 3. Let v = 4πr³/3 be the dimensionless volume of a bubble.

We begin the simulation with N = 3 × 10⁷ bubbles with different radii r_i (and volumes v_i). The initial bubble number N is chosen to be as large as is feasible (we use a PC with a 16-core Intel i7 processor from circa 2021, 16GB of RAM, and an SSD), noting that memory usage appears to be the limiting factor. We make no attempt to parallelise the code. For comparison, Thomas et al.¹⁶ use initial foams containing about 2 × 10⁶ bubbles in their 3D Potts-model simulations of dry-foam coarsening (which account for the detailed foam structure, unlike our mean-field simulations). We sample the relative radii r_i/〈r〉 from a specified distribution, usually a narrow²⁹ lognormal distribution with standard deviation σ_R = 1/10. The radii r_i themselves are determined by imposing the condition (see Section 3) that the total gas volume is . We select a lognormal distribution since this is frequently^17,36,67 used to approximate the scaling-state distribution of three-dimensional coarsening systems with ϕ = 0, although as an initial condition the choice is fairly arbitrary. See Section 4.2 for a discussion of the above choice of σ_R, along with the other initial distributions we use.

Next, the main loop is entered, each pass through which evolves the foam by a single time step Δτ (to be defined below). At the beginning of each step, the critical radius r_c is calculated numerically by solving 〈dv/dτ〉 = 0; i.e. by imposing that the total volume of gas is conserved.^2,6

Once r_c is obtained, the growth rate dr/dτ or dv/dτ of each bubble is known from eqn (12). The root-mean-square volumetric growth rate is used to set the time step Δτ, which we define such that m = 10³ steps (this parameter has been varied to check convergence; see Appendix B) would be required for a hypothetical bubble initially with the mean volume 〈v〉 to shrink to zero size if its growth rate were fixed at the current value of . We use this adaptive definition of Δτ for numerical efficiency, noting that the rate of coarsening (measured by dr_c/dτ for example) is expected to slow with time in a real foam (since r_c ∼ τ^α for α < 1 in the scaling state;^8,18 in our model, α = 1/2 by Section 3).

Statistics of the foam, including the mean bubble radius 〈r〉 and polydispersity⁴⁵ [scr P, script letter P], are then calculated and output. Every 300 time steps, we store the bubble radius distribution as a histogram in r/r_c, with 100 equal bins between r = 0 and r/r_c = 2.5, which we find sufficient to characterise the distribution (indeed, we undersample the histogram in Section 4.2).

A single Euler step⁶⁸ is then applied to evolve the bubble volumes; i.e. v_i(τ + Δτ) = v_i(τ) + (dv_i/dτ)Δτ. More efficient algorithms might allow fewer (larger) time steps to be used without increasing error in the bubble radii. However, the execution time of a simulation is fairly short (requiring about 4 hours on the PC described above). Furthermore, we believe that the main inaccuracy induced by the finite time steps is due to the size threshold for small bubbles, which we now discuss.

We recall that small bubbles shrink and eventually vanish due to coarsening.⁶ In the simulations, we define a minimum bubble volume allowed following the above Euler step, below which a bubble is deleted. In order to avoid bubbles of negative volume, we delete a bubble if its volume v is small enough that the bubble would vanish after two further time steps Δτ at its present growth rate dv/dτ (both Δτ and dv/dτ are fixed at their values calculated prior to the Euler step). The choice of two steps is arbitrary, but is intended to allow for the variation of Δτ in future steps (recalling that this is chosen adaptively). The total volume of all the deleted bubbles is reassigned to the bubbles that remain by multiplying their volumes by a constant factor such that the total remains unity.^6,69 The minimum bubble volume affects the bubble-size distribution at small radii, as discussed in Appendix B.

This completes the time step, following which the next pass of the main loop is performed. We continue until fewer than 10³ bubbles remain, at which point the simulation is halted since we expect smaller systems to exhibit excessive statistical fluctuations¹⁶ (see Fig. 9 and 10 below).

4.2 Simulation results

Using the methods⁴¹ described in Section 4.1, we perform coarsening simulations to first verify that the distribution (27) is reproducible at a variety of liquid fractions ϕ for a given initial bubble-size distribution. We then vary the initial distribution to provide evidence that the derived scaling state is universal, as it is expected to be in a real foam.² Finally, we compare the evolution of other foam properties (such as the critical radius R_c) with theory. In Appendix B, we describe exceptions to the universality of the scaling state,^29,42 and discuss convergence of our simulations with respect to the time-step size.

We emphasise that the simulations use the same growth law (10) as analysed in Section 3. Hence, the growth law itself and the modelling assumptions are untested by the simulation results. Rather, our aim here is to show that the derived scaling state is approached at late times over a range of ϕ and from a variety of initial conditions (which is not proven by the derivation).

In Fig. 6, we plot bubble-size histograms from our simulations for ϕ ∈ {1%,10%,35%}. The initial bubble sizes are taken from a narrow lognormal distribution with standard deviation σ_R = 1/10 with respect to the relative radius R/〈R〉. Data from five simulations is aggregated in each case (with distinct samples from the same initial distribution), and we plot histograms when about 10⁵ bubbles remain and when about 10³ bubbles remain (i.e. at the end of the simulations). We observe from Fig. 6 that agreement with eqn (27) is good at each liquid fraction (whereas the initial distribution is dissimilar from the predicted scaling state), and the nonzero value of ρ(R̂ = 0) is reproduced to a fair approximation.

Fig. 6 Simulated bubble-size histograms for (a) ϕ = 1%, (b) ϕ = 10%, and (c) ϕ = 35%. The mean over five simulation runs is shown, and error bars give the sample standard deviation (where larger than the marker) over the five runs. The distributions are plotted when roughly N = 10⁵ bubbles remain (with 20 bins), and when the runs halt (with 10 bins). The corresponding number τ_n of completed time steps is given. The final N value is not exactly 10³ because multiple bubbles can disappear in one time step. The uncertainties in the legend are due to differences between the five runs. Comparison is made to the predicted scaling state (27), and to that for the dry limit, the Lemlich–Markworth distribution (29).^12,36,46 The initial distribution (stated in the text) is also shown (not normalised, for clarity), as is the value of R_c/R₃₂ for each histogram, alongside the theoretical prediction obtained as in Fig. 4.

Therefore, the derived scaling state is approached from one initial condition for a range of liquid fractions. In Fig. 7, we vary the initial bubble-size distribution (defined in the caption) to show that the theoretical prediction is approached from a variety of initial conditions (similar plots are given by Li et al.²⁹ for a different growth law). Comparison should be made with Fig. 6(b). We note that the speed with which the system approaches the scaling state (judged by eye using the deviation from theory at the end of the simulations) varies with the initial distribution, and we see that wider distributions tend to result in a longer transient, as observed in other coarsening simulations with different growth laws (though exceptions have previously been observed for extremely narrow initial distributions).^29,41 This is the reason we select a narrow initial distribution in Section 4.1 (lognormal with standard deviation σ_R = 1/10) for use in most of our simulations.

Fig. 7 Simulated bubble-size histograms for ϕ = 10%. These plots differ from Fig. 6(b) only in the initial bubble-size distributions (which are again shown unnormalised for clarity). In (a), the initial distribution is lognormal with standard deviation σ_R = 1/5 in R/〈R〉. The distribution is lognormal with σ_R = 1/100 in (b), and is exponential in (c).

Having shown that the scaling-state distribution (27) is approached from a variety of initial conditions, we now plot the time-dependence of a few foam statistics as a further test of the derivation in Section 3 (again, we emphasise that this is not a test of the modelling assumptions). We start with the evolution of critical radius r_c for a variety of liquid fractions, in Fig. 8.^8,18,21 Good agreement is observed with the scaling-state prediction x_cτ^1/2 after a transient of the expected qualitative form;^35,69 where x_c is given by eqn (24). The offset of the curves corresponds to the decrease in coarsening rate with liquid fraction (cf.Fig. 2 of Pasquet et al.¹⁸). The relative difference of r_c from x_cτ^1/2 is plotted in Fig. 9 for various initial conditions at ϕ = 10% to reinforce this agreement. All plotted runs exhibit a relative difference of below 0.1 at their termination. Consistent with Fig. 7, narrower initial distributions tend to approach the theoretical scaling state more closely by the end of the simulations;^29,41i.e. they typically exhibit a smaller deviation between r_c and x_cτ^1/2. The two cusps in the curves for the narrowest initial condition correspond to changes in the sign of r_c − x_cτ^1/2: very narrow initial conditions can cause r_c to undershoot its scaling-state prediction, consistent with the theory of Chieco and Durian.³⁵

Fig. 8 Dimensionless critical radius r_cversus dimensionless time τ, on logarithmic scales, from simulations at different liquid fractions ϕ. Comparison is made to the scaling-state derivation of Section 3, where x_c is from eqn (24), as a test of the derivation itself (not of the modelling assumptions). The initial bubble-size distributions are lognormal with σ_R = 1/10, as in Fig. 6. The simulation curves start after ten time steps, to ensure they are sufficiently resolved.

Fig. 9 Relative difference between the simulated critical radius r_c and the scaling-state prediction x_cτ^1/2, versus time. Eqn (24) gives x_c, and the liquid fraction is ϕ = 10%. Runs are plotted for lognormal initial bubble-size distributions with various standard deviations σ_R; five independent runs are plotted per value of σ_R, which overlap until late times. The curves start after ten time steps, as in Fig. 8.

Finally, Fig. 10 compares the polydispersity⁴⁵ [scr P, script letter P] (see Section 3) with its predicted scaling-state value for several liquid fractions. The critical radius r_c is used as the independent variable here, to avoid an offset between the curves. We see that the simulations tend towards the theoretical predictions, which are obtained via numerical integration of eqn (27). The approach to these scaling-state values depends on the initial conditions (data not shown).

Fig. 10 Polydispersity⁴⁵ [scr P, script letter P] versus dimensionless critical radius, for various liquid fractions ϕ. As in Fig. 6 and 8, the initial bubble-size distribution is lognormal with σ_R = 1/10. Five runs per ϕ are plotted, and the curves start after ten time steps. Comparison is made with the prediction of [scr P, script letter P] in the scaling state for each plotted ϕ, which is obtained from the definition [scr P, script letter P] = R₃₂/〈R³〉^1/3 − 1 and the scaling-state distribution (27) by numerical integration.

In Fig. 9 and 10, we have plotted five simulation runs per liquid fraction (differing in the initial bubble-size samples), since this conveniently indicates, by the degree to which the curves ‘fray,’ the importance of statistical fluctuations at a given time.

The results of this section show that the scaling state derived in Section 3 is approached from a range of initial conditions (noting the exceptions^29,42,43 discussed in Appendix B) in simulations which apply the mean-field growth law (10) directly. This validates the derivation of Section 3 and indicates that the derived scaling state is stable.

5 Discussion

The scaling-state bubble-size distribution (27) derived from the mean-field growth law (10) is qualitatively different for ϕ > 0 from the distribution (29) predicted by Lemlich's law. The distribution ρ(R) described by eqn (27) is nonzero at R = 0 for ϕ > 0, whereas ρ(0) = 0 for eqn (29). To show this directly, we use the following result (which we have expressed in terms of R): the relative rate at which bubbles disappear is^13,65

Since dN/dt < 0 during coarsening, and dR/dt takes a finite negative value at R = 0 for [capital Pi, Greek, circumflex] > 0 (and ϕ > 0) by eqn (10), we must have ρ(0) > 0 for our growth law. For Lemlich's law (1), dR/dt diverges as R → 0, and so ρ(0) = 0. The film area vanishes in this limit for ϕ > 0 by eqn (9), thus suppressing the singularity in our model due to the assumption of border blocking.

Simulations²¹ and experiments^22,23 have found a large population of small bubbles in coarsening wet foams, in 2D and 3D respectively. However, the forms of the observed bubble-size distributions (which are similar to each other) differ from eqn (27): they have large peaks at small bubble radius, which decrease in height as ϕ increases,^21–23,54 †† whereas the value of the distribution at R = 0 increases with ϕ in our prediction (see Fig. 5). The differing form of our distribution also means that our predicted polydispersities [scr P, script letter P] (see Fig. 4 and 10) differ from those observed in experiments (see the supplementary information of Galvani et al.²²): the large peak and long tail in their distributions at moderate ϕ result in a larger [scr P, script letter P] than we find, while the reduction in the size of the peak with increasing ϕ causes [scr P, script letter P] to decrease,²² whereas we find it to increase.

The experiments²² are not wholly comparable with our model, since inter-bubble adhesion is believed to have been present in the former.¹⁸ However, such adhesion should increase the shrinkage rate of small bubbles due to the increased contact between bubbles,^22,70 and thus act to reduce the population of small bubbles. Also, we recall that our model neglects gas flow through the bulk liquid, and so we can only compare directly with the simulations reported by Khakalo et al.²¹ which also omit this contribution.⁵⁴

Both Khakalo et al.²¹ and Galvani et al.²² attribute the peaks in their distributions to the presence of rattlers, i.e. bubbles that are not compressed by their neighbours (though they may be adhered to them in the presence of a contact angle²²). In the absence of adhesion, these bubbles lose contact with their neighbours, their films thus shrinking to zero size before their radius has done so (resulting in a border-blocking growth rate of zero).^21,69 Our mean-field model does not predict rattlers of finite size, since A/S > 0 for all R > 0 by eqn (9). However, rattlers might be present in the finite-element simulation results^26,40 shown in Fig. 1(b) and 2(b), as small bubbles with (almost) zero growth rate.

We now consider the possible roles of rattlers in a mean-field coarsening model. One possibility, suggested in 2D by the binned mean of Fig. 1(b), is that very small bubbles have nonzero growth rates, described by eqn (10) or another growth law, when the scatter at fixed R is averaged. In other words, there are no rattlers of finite size in an averaged description, and all rattlers that exist in a real foam are a result of scatter away from the average (such scatter is apparent in Fig. 1 and 2)—differing bubble environments mean that some small bubbles become rattlers, and cease to engage in coarsening (at least for a time²¹), leading to an accumulation of rattlers and hence peaked distributions similar to those observed by Khakalo et al.²¹ and Galvani et al.²² It seems unlikely that such distributions could be captured by any mean-field model in this case.

Another possibility is that the averaged growth rate does reach zero at a finite bubble size; i.e. rattlers exist in the mean field. Fig. 2 suggests this in 3D, although possibly as an artefact of our aggregation of data from several small foams (see Section 2). A stronger argument comes from inspecting Fig. 1(a) closely: extrapolating the trend of the smallest cluster of bubbles (those with three neighbours⁹) suggests that the growth rate reaches zero for R > 0 in 2D (indeed, this is predicted by the theory of Roth et al.⁹ and Schimming and Durian⁷). Galvani et al.²² argue for such a cutoff radius in 3D (below which bubbles are rattlers) based on the size of a Plateau-border junction, and their experimental results support this. However, rattlers would cease to evolve in the absence of adhesion or bulk gas transfer, and so no bubbles would vanish as a result of coarsening. Thus, a conventional scaling state would be precluded for such a border-blocking model without bubble adhesion. This does not contradict the experimental results of Galvani et al.,²² due to the likely presence of adhesion (and bulk gas transfer).‡‡

The second possibility, that the averaged bubble growth rate reaches zero at finite radius, appears more likely to us. However, the simulations of Khakalo et al.²¹ (particularly the results included in their preprint⁵⁴) suggest that there is a conventional scaling state when adhesion and bulk gas transfer are omitted, and hence that there is no cutoff size below which all bubbles are rattlers. Recalling Fig. 1(a), it is conceivable that the growth-rate distribution shown varies with polydispersity—perhaps, as rattlers accumulate, the cutoff size decreases, thus ensuring their eventual disappearance. It may be fruitful to explore the relation between coarsening wet foams and the highly bidisperse soft disk systems studied by Petit and Sperl,⁷¹ who find a second jamming transition at ϕ < ϕ_c for the small disks.

It would be very useful for the development of mean-field models of coarsening to distinguish between these (and any other) possibilities, perhaps through further bubble-scale simulations. Three-dimensional finite-element simulations^40,55 of larger and more polydisperse foams (hence containing more small bubbles) could be performed, although improvement to our current numerics would be needed to make larger systems feasible for us. Two-dimensional finite-element simulations^26,49,72 of more polydisperse foams may also be useful, as might simulations tracking the evolution under coarsening of small bubbles. For comparability, all these simulations should be without bubble adhesion.

Bubble-model^21,30,69,73 simulations in 2D or 3D would seem to be a promising approach, due to the feasibility of simulating large wet foams. It would be interesting to see whether scaling-state bubble-size distributions with rattlers omitted are similar to our prediction (27). Further wet-foam coarsening experiments, such as those described by Galvani et al.,²² in which individual bubbles can be tracked, may also be helpful in clarifying the dynamics of rattlers.

6 Conclusion

In the present work, we derived^11,25 the scaling-state bubble-size distribution (27) predicted by our proposed mean-field border-blocking growth law (10), for foams with zero contact angle and liquid fraction ϕ < ϕ_c. Eqn (10) was validated against averaged data from bubble-scale simulations^26,40 for ϕ = 10% and 30%, and its 2D equivalent (7) was validated for ϕ = 2% and 10%. We showed that the scaling-state growth exponent is α = 1/2 at all ϕ < ϕ_c, as expected since the model incorporates only gas transfer through contact films between bubbles.¹⁸ Using mean-field simulations (with methods adapted from De Smet et al.⁴¹) in which bubbles evolve according to the proposed growth law (10), we also checked that the derived scaling state is approached from various initial conditions (we discuss some apparent exceptions in Appendix B, which are judged to be artefacts of the mean-field approach^29,42,43).

We then compared our predicted bubble-size distribution (27) with prior experimental results^22,23 and simulations using the bubble model.²¹ While, consistently with these, eqn (27) reproduces a large population of small bubbles for ϕ > 0, it exhibits qualitative differences from the previous observations, in which there is a large peak in the distribution at small bubble radii.^21–23 This peak has previously been attributed to rattlers^21,22—small bubbles not pressed into contact with their neighbours—whereas our model does not incorporate rattlers of finite size. The addition of rattlers would be an important refinement, which could be developed using the results of future bubble-scale simulations as discussed in Section 5. Nevertheless, we suggest that our predicted distribution (27) serves as an initial approximation incorporating variation with ϕ (noting that any refinements to our model would likely preclude an analytical derivation), which should allow the qualitative influence of additional effects such as gas transfer through bulk liquid^7,19,24 or finite contact angle^18,32 to be studied by adding approximations for these to the growth law (10). The numerical methods⁴¹ of Section 4.1 can be straightforwardly adapted to use other growth laws. A model for bulk gas transfer may also allow estimation of the liquid fraction above which the border-blocking assumption is a poor approximation, though this will vary with^7,18,24h/〈R〉. Furthermore, it may be fruitful to study whether the tail of the bubble-size distribution is sensitive to the form of the growth law at small radii.⁴²

A version of eqn (10) augmented with an approximation for bulk-liquid gas transfer^7,24 could be used to predict the variation of growth exponent α with ϕ, by adapting the approach of White.²⁸ This may provide a prediction without fitting parameters, thus developing the existing model of Durian.²⁴

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data and simulation code described in this work (including the 3D finite-element simulations⁴⁰) are available from the Aberystwyth University Research Portal (https://doi.org/10.20391/59c8e4a2-08ce-40ce-b3e7-0b7befdcd3f4).

Appendices

Appendix A

In this appendix, we provide derivations which form part of the arguments leading to the proposed 3D mean-field growth law (10), along with the first few moments of the scaling-state bubble size distribution (27) that this growth law predicts in turn.

Derivation of bubble pressure and film area approximations. The following comprises derivations of eqn (8) and (9). The arguments are very similar to those in 2D,²⁶ and are given in more detail by ref. 40. We recall that eqn (8) and (9) generalise results of Höhler et al.³⁹ to individual bubbles in polydisperse foams (although we make additional approximations); our derivations use similar methods to Höhler et al.³⁹

We start with eqn (8), which approximates the pressure p of a bubble (measured relative to the liquid pressure). Let the bubble have effective radius R, and liquid–gas interface Σ with unit normal n, and let δ_ij denote the Kronecker delta. Volume-averaging the stress tensor τ over the single bubble (including its gas and liquid–gas interface, but no liquid) gives^39,74,75

Let p̄ = −(Tr[small tau, Greek, macron])/3 be a bubble-level analogue of the osmotic pressure.^39,76Eqn (32) then gives, after approximating the bubble's surface area S by 4πR²,⁴⁵

As a first approximation, let us assume that p̄ is equal for all bubbles (i.e. the environment of every bubble is similar; some evidence supporting this in emulsions is given by Fig. 2(d) of Jorjadze et al.⁷⁶). Substituting eqn (33) into the relation^39,77Π_c = Π_O/(1 − ϕ) + 2γ/R₃₂ between the foam's capillary and osmotic pressures then gives p̄ ≈ Π_O/(1 − ϕ). By substituting this into eqn (33) we obtain eqn (8), which was validated against the finite-element simulations of Fig. 2 in ref. 40.

We now derive eqn (9), for the ratio of a bubble's film area A to its total surface area S (i.e. the proportion of its interface in contact with other bubbles). Let us retain the above notation, and define f as the force per unit area applied to Σ at point r by the bubble's contacting neighbours (we set the liquid pressure to zero without loss of generality). We start the derivation from an alternative expression of the bubble's volume-averaged stress tensor:⁷⁸

On parts of Σ adjoining contact films, we have f = −Π_Dn, where Π_D is the film's disjoining pressure.² Approximating the pressure of each neighbouring bubble by the capillary pressure Π_c gives⁷⁹Π_D = (p + Π_c)/2. Herein, we substitute eqn (8) for p and³⁹Π_c = Π_O/(1 − ϕ) + 2γ/R₃₂. On parts of Σ not adjoining films, f = 0.

Next, we substitute these results for f into eqn (34), and assume that the bubble geometry does not differ too much from a sphere, so that r ≈ Rn. Approximating 4πR² ≈ S, we obtain

using the above definition p̄ = −(Tr[small tau, Greek, macron])/3, and R̂ = R/R₃₂ and [capital Pi, Greek, circumflex] = Π_O(R₃₂/γ)/(1 − ϕ) from Section 2 (prior to the redefinition of R̂ and [capital Pi, Greek, circumflex] discussed at the end of that section). Again approximating p̄ ≈ Π_O/(1 − ϕ) for all bubbles in the foam, and rearranging eqn (35), gives eqn (9). Like eqn (8), this result was validated in ref. 40 using finite-element simulations.

Normalisation and moments of bubble size distribution. We now provide the normalisation constant C of the bubble size distribution ρ(R̂) derived in Section 3 (and given by eqn (27)), along with this distribution's first few moments 〈R̂^j〉 for j = 1, 2, 3. Expressions for these can be obtained directly by evaluating the definition with the substitution s = 3[R̂₀ + 1/(1 + [capital Pi, Greek, circumflex])]/(R̂₀ − R̂). The constant C is obtained by imposing normalisation; i.e. 〈R̂⁰〉 = 1. Let a = s(R̂ = 0) = 3[1 + (1/R̂₀)/(1 + [capital Pi, Greek, circumflex])]. Then

We see that 〈R̂^j〉 can be evaluated analytically for j = 0, 1, 2, whereas the third and higher-order moments are expressible only in terms of the exponential integral, and thus must be calculated numerically. We note that a, and hence the moments, may be expressed wholly in terms of [capital Pi, Greek, circumflex] viaeqn (28).

Appendix B

Here we discuss convergence of our mean-field simulations (see Section 4) with respect to time-step size Δτ, along with additional apparent scaling states which arise from pathological initial conditions.

Simulation convergence. We recall the definition of Δτ from Section 4.1: Δτ is such that m steps of this size would be needed for a hypothetical bubble initially with the mean volume to vanish if its shrinkage rate were fixed at the current root-mean-square growth rate. We set m = 10³ in Section 4.1. Fig. 11 shows the effect on the bubble-size histograms of increasing m by a factor of 10 (i.e. decreasing Δτ by the same factor), from which we see little change in all bins except that for smallest radius R (the plot is for liquid fraction ϕ = 35% since large ϕ corresponds to the largest population of small bubbles by Fig. 5). To interpret this latter variation, we recall from Section 4.1 that Δτ determines the size below which bubbles are deleted (and their gas redistributed): the minimum volume v is that for which the bubble would vanish at its present shrinkage rate after a time of exactly 2Δτ. Hence, the size threshold increases with Δτ—if the time step is larger, bubbles will be deleted earlier in their shrinkage process, and so fewer bubbles will be present in the bin of smallest radius, in agreement with Fig. 11. This also lowers the polydispersity of the distribution slightly (not shown). Excepting this minor effect on the smallest bin of the histogram, we judge our simulations to be converged with respect to Δτ. Another method for deleting small bubbles (by waiting until their volume drops below zero, for example) might reduce the effect on the smallest bin.

Fig. 11 Bubble-size histogram when approximately 10⁵ bubbles remain in simulations of the type described in Section 4. The liquid fraction is 35%, and results are given for two values of the parameter m defining the time step Δτ (see Section 4.1). Larger m corresponds to smaller Δτ. As in Fig. 6, the two histograms give the average over five distinct runs (the sample standard deviations are smaller than the markers). The initial bubble-size distribution was lognormal with σ_R = 1/10 in both cases. The predicted scaling-state distribution (27) is shown for comparison.

Additional apparent scaling states. We now discuss some apparent exceptions we observe to the universality of the scaling state derived in Section 3, which are interpreted to result from pathological initial conditions.^29,42,43 For a triangular initial bubble-size distribution, which is not continuously differentiable, we see from Fig. 12 that the system appears to approach a different scaling-state distribution. Like Brown,⁴² who studied the LSW law (2), we also find, in our simulations, that a discontinuous initial bubble size distribution (we used a uniform distribution) gives an apparent scaling state which likewise has a discontinuity (not shown). Brown⁴² argued that the scaling state of the LSW law is not unique, and derived a family of scaling-state distributions. It is now believed^29,43 that, with the exception of the previously known scaling state of the LSW law,¹¹ the family found by Brown⁴² consists of unstable scaling states. Indeed, when τu is plotted against R̂ (see Section 3; plots not shown) for the apparent scaling state of Fig. 12, two zeros are observed, with the apparent cutoff R̂₀ near the smaller of these. Thus, this scaling state should be unstable to perturbation by the arguments of Section 3.^25,43 It appears that mean-field simulations of the type we use^41,42 do not supply sufficient perturbations to the distribution during coarsening for the system to exit the spurious scaling state, at least over the timescale of our simulations.^29,43 Mean-field laws specify an equal growth rate for all bubbles of equal R, whereas real foams exhibit a large amount of scatter in the growth rates^59,80 (see our Fig. 1 and 2 for simulations thereof).

Fig. 12 Bubble-size histograms for simulations with ϕ = 10% and an initially triangular radius distribution. The latter is symmetric and has σ_R = 1/10. The data is presented in the same format as for Fig. 6.

We therefore neglect these apparent exceptions to the universality of the scaling state as artefacts of the mean-field growth law (10) which would not arise in a real foam. Modelling approaches which add a diffusion term to the continuity eqn (13) may lack these artefacts.^61,81

Acknowledgements

The authors acknowledge helpful discussions with Sylvie Cohen-Addad, Tudur Davies, François Graner, Reinhard Höhler, and Rory McCranor. JM funded by Cronfa Ymchwil Joy Welch Research Fund and ESA REFOAM MAP contract 4000129502. SC partially supported by the Horizon 2020 Framework Programme for Research and Innovation, Grant Agreement number 101008140 EffectFact. JM and SC funded by EPSRC award number UKRI2390.

References

1. D. Weaire and S. Hutzler, The Physics of Foams, Oxford University Press, Oxford (GB), 1999.
2. I. Cantat, S. Cohen-Addad, F. Elias, F. Graner, R. Höhler, O. Pitois, F. Rouyer and A. Saint-Jalmes, Foams: Structure and Dynamics, Oxford University Press, Oxford (GB), 2013.
3. P. S. Stewart and S. Hilgenfeldt, Annu. Rev. Fluid Mech., 2023, 55, 323–350.
4. T. J. Martin, in Foam Engineering: Fundamentals and Applications, ed. P. Stevenson, John Wiley & Sons Ltd, Chichester (GB), 2012, ch. 17, pp. 411–457.
5. P. Stevenson and N. W. A. Lambert, in Foam Engineering: Fundamentals and Applications, ed. P. Stevenson, John Wiley & Sons Ltd, Chichester (GB), 2012, ch. 11, pp. 229–249.
6. R. Lemlich, Ind. Eng. Chem. Fundam., 1978, 17, 89–93.
7. C. D. Schimming and D. J. Durian, Phys. Rev. E, 2017, 96, 032805.
8. N. Isert, G. Maret and C. M. Aegerter, Eur. Phys. J. E:Soft Matter Biol. Phys., 2013, 36, 116.
9. A. E. Roth, C. D. Jones and D. J. Durian, Phys. Rev. E, 2013, 87, 042304.
10. P. Born, M. Braibanti, L. Cristofolini, S. Cohen-Addad, D. J. Durian, S. U. Egelhaaf, M. A. Escobedo-Sánchez, R. Höhler, T. D. Karapantsios, D. Langevin, L. Liggieri, M. Pasquet, E. Rio, A. Salonen, M. Schröter, M. Sperl, R. Sütterlin and A. B. Zuccolotto-Bernez, Rev. Sci. Instrum., 2021, 92, 124503.
11. I. M. Lifshitz and V. V. Slezov, J. Exp. Theor. Phys., 1959, 35, 331–339.
12. C. Wagner, Z. Elektrochem., 1961, 65, 581–591.
13. W. W. Mullins, J. Appl. Phys., 1986, 59, 1341–1349.
14. J. A. Glazier and D. Weaire, J. Phys.: Condens. Matter, 1992, 4, 1867–1894.
15. J. Stavans, Rep. Prog. Phys., 1993, 56, 733–789.
16. G. L. Thomas, R. M. C. de Almeida and F. Graner, Phys. Rev. E, 2006, 74, 021407.
17. J. Lambert, R. Mokso, I. Cantat, P. Cloetens, J. A. Glazier, F. Graner and R. Delannay, Phys. Rev. Lett., 2010, 104, 248304.
18. M. Pasquet, N. Galvani, A. Requier, S. Cohen-Addad, R. Höhler, O. Pitois, E. Rio, A. Salonen and D. Langevin, Soft Matter, 2023, 19, 6267–6279.
19. I. Fortuna, G. L. Thomas, R. M. C. de Almeida and F. Graner, Phys. Rev. Lett., 2012, 108, 248301.
20. G. L. Thomas, J. M. Belmonte, F. Graner, J. A. Glazier and R. M. C. de Almeida, Colloids Surf., A, 2015, 473, 109–114.
21. K. Khakalo, K. Baumgarten, B. P. Tighe and A. Puisto, Phys. Rev. E, 2018, 98, 012607.
22. N. Galvani, M. Pasquet, A. Mukherjee, A. Requier, S. Cohen-Addad, O. Pitois, R. Höhler, E. Rio, A. Salonen, D. J. Durian and D. Langevin, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2306551120.
23. N. Galvani, O. Pitois and S. Cohen-Addad, J. Colloid Interface Sci., 2025, 696, 137661.
24. D. J. Durian, arXiv, 2023, preprint, arXiv:2304.00415 DOI:10.48550/arXiv.2304.00415.
25. J. A. Marqusee and J. Ross, J. Chem. Phys., 1983, 79, 373–378.
26. J. Morgan and S. Cox, J. Fluid Mech., 2024, 999, A10.
27. F. Bolton and D. Weaire, Philos. Mag. B, 1991, 63, 795–809.
28. R. J. White, Mater. Sci. Eng., 1979, 40, 15–20.
29. Y. Li, Z. Wang, X. Gao, Y. Wang, J. Li and J. Wang, Acta Mater., 2022, 237, 118196.
30. B. S. Gardiner, B. Z. Dlugogorski and G. J. Jameson, J. Phys.: Condens. Matter, 1999, 11, 5437–5453.
31. H. M. Princen, J. Phys. Chem., 1968, 72, 3342–3345.
32. S. J. Cox, A. M. Kraynik, D. Weaire and S. Hutzler, Soft Matter, 2018, 14, 5922–5929.
33. H. M. Princen, J. Colloid Interface Sci., 1979, 71, 55–66.
34. H. W. Yarranton and J. H. Masliyah, J. Colloid Interface Sci., 1997, 196, 157–169.
35. A. T. Chieco and D. J. Durian, Phys. Rev. E, 2023, 108, 034606.
36. M. Hillert, Acta Metall., 1965, 13, 227–238.
37. A. J. Ardell, Acta Metall., 1972, 20, 61–71.
38. P. Streitenberger, Acta Mater., 2013, 61, 5026–5035.
39. R. Höhler, J. Seknagi and A. Kraynik, Soft Matter, 2021, 17, 6995–7003.
40. J. Morgan, PhD thesis, Aberystwyth University, 2025.
41. Y. De Smet, L. Deriemaeker and R. Finsy, Langmuir, 1997, 13, 6884–6888.
42. L. C. Brown, Acta Metall., 1989, 37, 71–77.
43. J. Svoboda, G. A. Zickler, E. Kozeschnik and F. D. Fischer, Acta Mater., 2022, 235, 118085.
44. H. M. Princen, Langmuir, 1986, 2, 519–524.
45. A. M. Kraynik, D. A. Reinelt and F. van Swol, Phys. Rev. Lett., 2004, 93, 208301.
46. A. J. Markworth, J. Colloid Interface Sci., 1985, 107, 569–571.
47. K. A. Brakke, Exp. Math., 1992, 1, 141–165.
48. K. A. Brakke, Surface Evolver v2.70a, 2013, https://kenbrakke.com/evolver/evolver.html, (accessed August 2023).
49. T. Kähärä, T. Tallinen and J. Timonen, Phys. Rev. E, 2014, 90, 032307.
50. A. Boromand, A. Signoriello, F. Ye, C. S. O'Hern and M. D. Shattuck, Phys. Rev. Lett., 2018, 121, 248003.
51. J. Morgan and S. Cox, Data for: Effects of liquid fraction and contact angle on structure and coarsening in two-dimensional foams, 2023, Aberystwyth University (PURE) DOI:10.20391/b2be2c00-151e-43df-90ed-5c900af55160.
52. SciPy contributors, SciPy v1.10.1, 2023, https://scipy.org, (accessed July 2025).
53. D. Mukesh, A. K. Das and P. K. Ghosh, Langmuir, 1992, 8, 807–814.
54. K. Khakalo, K. Baumgarten, B. P. Tighe and A. Puisto, arXiv, 2017, preprint, arXiv:1706.03932 DOI:10.48550/arXiv.1706.03932.
55. D. Wang, J. D. Treado, A. Boromand, B. Norwick, M. P. Murrell, M. D. Shattuck and C. S. O'Hern, Soft Matter, 2021, 17, 9901–9915.
56. P. Van Liedekerke, J. Neitsch, T. Johann, E. Warmt, I. Gonzàlez-Valverde, S. Hoehme, S. Grosser, J. Kaes and D. Drasdo, Biomech. Model. Mechanobiol., 2020, 19, 189–220.
57. R. Quey, P. Dawson and F. Barbe, Comput. Method. Appl. Mech. Eng., 2011, 200, 1729–1745.
58. R. Quey, Neper v4.5.0, 2022, https://neper.info, (accessed Sep 2025).
59. J. Lambert, I. Cantat, R. Delannay, R. Mokso, P. Cloetens, J. A. Glazier and F. Graner, Phys. Rev. Lett., 2007, 99, 058304.
60. A. Maestro, W. Drenckhan, E. Rio and R. Höhler, Soft Matter, 2013, 9, 2531–2540.
61. K. Binder, Phys. Rev. B, 1977, 15, 4425–4447.
62. P. Streitenberger and D. Zöllner, Scr. Mater., 2006, 55, 461–464.
63. D. Zöllner and P. Streitenberger, Scr. Mater., 2006, 54, 1697–1702.
64. A. Requier, S. Varade, N. Galvani, M. Pasquet, S. Cohen Addad, C. Gehin-Delval, R. Höhler, O. Pitois, E. Rio, A. Salonen and D. Langevin, J. Colloid Interface Sci., 2025, 696, 137825.
65. J. Lambert, F. Graner, R. Delannay and I. Cantat, Europhys. Lett., 2012, 99, 48003.
66. NumPy contributors, NumPy v1.24.3, 2023, https://numpy.org, (accessed July 2025).
67. R. Quey and L. Renversade, Comput. Method. Appl. Mech. Eng., 2018, 330, 308–333.
68. R. A. Adams and C. Essex, Calculus: A Complete Course, Pearson, Toronto (CA), 8th edn, 2013.
69. B. S. Gardiner, B. Z. Dlugogorski and G. J. Jameson, Philos. Mag. A, 2000, 80, 981–1000.
70. C. Feng, Z. Jing, X. Ma and H. Wang, Int. J. Multiphase Flow, 2021, 143, 103782.
71. J. C. Petit and M. Sperl, Phys. Rev. Res., 2025, 7, L022073.
72. A. Boromand, A. Signoriello, J. Lowensohn, C. S. Orellana, E. R. Weeks, F. Ye, M. D. Shattuck and C. S. O'Hern, Soft Matter, 2019, 15, 5854–5865.
73. D. J. Durian, Phys. Rev. Lett., 1995, 75, 4780–4783.
74. G. K. Batchelor, J. Fluid Mech., 1970, 41, 545–570.
75. M. E. Evans, A. M. Kraynik, D. A. Reinelt, K. Mecke and G. E. Schröder-Turk, Phys. Rev. Lett., 2013, 111, 138301.
76. I. Jorjadze, L.-L. Pontani and J. Brujic, Phys. Rev. Lett., 2013, 110, 048302.
77. H. M. Princen, Langmuir, 1988, 4, 164–169.
78. L. D. Landau, E. M. Lifshitz, A. M. Kosevich and L. P. Pitaevskiĭ, Theory of Elasticity, Elsevier, Oxford (GB), 3rd edn, 1986.
79. D. Exerowa and P. M. Kruglyakov, Foam and Foam Films: Theory, Experiment, Application, Elsevier, Amsterdam (NL), 1998, p. 90.
80. J. Dittmann, A. Eggert, M. Lambertus, J. Dombrowski, A. Rack and S. Zabler, J. Colloid Interface Sci., 2016, 467, 148–157.
81. D. Zimnyakov, A. Zemlyanukhin, S. Yuvchenko, A. Bochkarev, I. Slavnetskov, S. Gavrilov and D. Tumachev, Phys. D, 2019, 398, 171–182.

---

Footnotes

† There may be other exceptions. For example, it is conceivable that a population of non-evolving rattlers in a border-blocking model³⁰ could preclude a scaling state. This possibility is discussed further in Section 5.

‡ We define these variables partly for use in the simulations of Section 4.

§ This is the reason for replacing R₃₂ by R_c in eqn (10).³⁷ Otherwise, x_c is obtained in terms of x₃₂ = R₃₂τ^−α/L, which can only be obtained numerically.

¶ Linearity and homogeneity of eqn (18) ensures that the argument would not be invalidated if higher order terms were included in eqn (22).

|| Multiple roots arise, but only one has the required property that both x_c and x₀ are positive.

** However, integrating eqn (10) for small R shows that a small bubble does vanish in finite time. Although the film area goes to zero as R → 0, the bubble pressure diverges, and so dR/dt tends to a finite limit. We also emphasise that it is the probability density ρ at R = 0 that is nonzero, not the probability that a given bubble has zero size, which remains zero.

†† The preprint⁵⁴ of Khakalo et al.²¹ contains further plots of the simulated scaling-state distributions in the border-blocking case we consider.

‡‡ We note that the observed peaked distributions^21,22 can be reproduced qualitatively with the simulations described in Section 4.1 (at least during the coarsening transient), if R̂ in the numerator of eqn (9) is replaced by max{0, R̂ − c} (for a small constant c), and a small multiple of the LSW law (2) is added to the resulting growth law (to roughly approximate bulk gas transfer;²¹cf. also ref. 37). However, these adjustments are clearly ad hoc, and physical arguments are needed to fix their parameters.

---