The retraction of thin films, as described by the Taylor–Culick (TC) theory, is subject to widespread debate, particularly for films at the nanoscale. We use non-equilibrium molecular dynamics simulations to explore the validity of the assumptions used in continuum models by tracking the evolution of holes in a film. By deriving a new mathematical form for the surface shape and considering a locally varying surface tension at the front of the retracting film, we reconcile the original theory with our simulation to recover a corrected TC speed valid at the nanoscale.

Thin film rupture is of premier interest in numerous physical and chemical transport processes,1–9 ranging from disease transmission10 to planetary scale environmental and oceanic science.11,12 As such, the film rupture process has been investigated for over a century13–15 and remains an active area of debate. In their pioneering works, Taylor16 and Culick17 independently derived an equation for film retraction by considering that the growth of a nucleated hole of radius, R, is driven purely by the balance between the surface force, Fγ, and the inertia of the liquid mass, m, collected in a rim surrounding the hole, i.e., Fγ=ddt(mṘ). Applying the initial condition at (time) t = 0, R = 0 for finite Ṙ, one obtains, upon integration, the Taylor–Culick (TC) speed,
(1)
where h0 is the thickness of the unperturbed part of the film with density ρ and surface tension γ. In the work of Taylor16 and Culick,17, ϕ = 2, but a later work18 uses smaller values for thinner films. This surprisingly simple, yet elegant equation captures the overall retraction dynamics and corrects Dupre’s inappropriate energy balance assumption.19 However, its accuracy is debatable, and it markedly deviates for thinner films.18,20–25 Among these studies, the seminal work of McEntee and Mysels18 conclusively demonstrates that for films thinner than 100 nm, the growth rate of an expanding hole is distinctly slower than the theoretical prediction—the reasons for such slow retraction are investigated in this work.

Taylor’s resolution is based on a few key assumptions: (i) when rupture propagates, h0 remains uniform except at the rim and (ii) all liquid-mass of the expanding hole is collected inside the rim, i.e., m = ρπR2h0. In a more formal derivation, Culick made the assumptions of uniform thickness and constant surface tension everywhere. While the assumptions leading to Eq. (1) are not incorrect at the global scale, their local deviations become increasingly important as the thickness decreases down to nanometers.

A range of continuum-based models attempted to explain the film rupture process by consideration of the non-uniformity of film thickness,26 surface elasticity,7 surface tension gradients,27 Marangoni stresses,28 and viscous effects.24,29–31 With continuum models, any additional physics must be included by adding another conservation or balance law. In the present work, non-equilibrium molecular dynamics (MD) simulation is applied to understand the retraction process from the atomic scale.

Any rupture event must start at the smallest scale; therefore, molecular dynamics is an ideal technique to realize the origin of this process. It provides the full picture of a fluid down to atomic interactions; the interface dynamics, surface tension, and viscosity emerge as outputs from simulations. Hence, many fluid properties have been successfully studied by molecular dynamics. Koplik and Banavar32 used this technique to show the origins of dewetting, alongside the studies of film breakup,33,34 boiling and nucleation,35 and a range of cases.36–38 Kono et al.39 studied a two-dimensional film breakup with MD and observed a temperature rise due to the rapid rupture, which would cause a decrease in surface tension. In the present work, we model the full three-dimensional film breakdown and systematically outline all possible contributions from the atomic scale. Using local pressure measurement techniques, interface tracking, and dynamic time-evolving mappings, local surface tension is calculated. This uncovers the complete picture of surface force and interface shape, putting the thermal effects in context and exposing other key factors contributing to the dynamics. The TC model does not take into account any surface-species and their contributions to the surface properties of the film.40,41 By construction, the film investigated in this work is pure in nature and, hence, closest to the TC model. The present work thereby tests this model by removing external factors and shows that the classical theory fails for molecularly thin films even in the absence of surface active agents, and only when certain atomic-scale corrections are made, the TC equation can successfully predict the retraction process for films with thickness down to a few nanometers.

Films investigated in this study were modeled using a Lennard-Jones (LJ) fluid simulated in the extensively validated and verified Flowmol MD code.42 The initial simulation domain was a cubic box of dimensions Lx = 76.19 and Ly = Lz = 609.56 in LJ units, with periodic boundary conditions applied in all directions. These dimensions correspond to 207 nm wide films with h0 < 5 nm (for argon, σ = 0.34 nm43). The rest of this communication reports all quantities in LJ units. The middle 20% of the simulation box in the x direction was designated as liquid and setup with a density ρ ≈ 0.7, and the remaining was designated as vapor with a lower density ρ ≈ 0.01. The schematic in Fig. 1(a) shows the system modeled in this study where all atoms (3.54 × 106) are identical, forming a central liquid film coexisting with surrounding vapor at equilibrium. After equilibrating the film at T = 0.78 ± 0.03 in the canonical ensemble (NVT), a hole of initial radius R0 was induced in the otherwise stable film. The system was then set to evolve in the micro-canonical ensemble (NVE), and R was measured over time. The inset in Fig. 1(a) (see also the supplementary material, Fig. 1) shows the liquid atoms surrounding an expanding hole. Displaced atoms (in red) from the hole collect around the hole to form a rim. The radial-averaged film extracted from MD data, in Fig. 1(b), shows its temporal evolution—with the simplified schematic in Fig. 1(c).

FIG. 1.

Thin film retraction. (a) Computational domain considered for this study. A hole is nucleated at the center of the film in the yz plane. A cylindrical needle schematically shows the hole-poking process. Zoomed-in view of the hole: atoms in red denote the atoms initially in the hole at t = 0 that are displaced; faint schematic shows an idealization of the rim surrounding the hole, which borders the unperturbed film. (b) Radial average of a retracting film over time: R(t) is measured from the center of the hole to the tip of the rim. The film is unperturbed only beyond the elliptic (red-dotted lines) blob. (c) Simplified schematic of (b), film thickness, h0, is measured at the unperturbed part of the film.

FIG. 1.

Thin film retraction. (a) Computational domain considered for this study. A hole is nucleated at the center of the film in the yz plane. A cylindrical needle schematically shows the hole-poking process. Zoomed-in view of the hole: atoms in red denote the atoms initially in the hole at t = 0 that are displaced; faint schematic shows an idealization of the rim surrounding the hole, which borders the unperturbed film. (b) Radial average of a retracting film over time: R(t) is measured from the center of the hole to the tip of the rim. The film is unperturbed only beyond the elliptic (red-dotted lines) blob. (c) Simplified schematic of (b), film thickness, h0, is measured at the unperturbed part of the film.

Close modal

Surface deformations arising from thermal motion increase the likelihood of spontaneous rupture.44–46 Therefore, to find a stable and computationally tractable case, films with a range of thicknesses were tested without applying any external perturbation. A thickness of h0 ≈ 14 was selected as films thinner than this were observed to break spontaneously before reaching a constant retraction rate. Other computationally affordable thicknesses were also examined (see the supplementary material, Sec. I A).

The growth of spontaneously nucleated holes (by allowing holes to form naturally due to thermal fluctuations) was found to give similar dynamics, albeit spontaneous holes take considerably longer time to nucleate. Each reported case was averaged over at least three independent simulations. The short-time growth, being more susceptible to thermal noise, was averaged over ten independent simulations (for further details, see the supplementary material, Sec. I).

We chart the hole-growth as predicted by MD and compare to experiments, computational fluid dynamics, and TC theory. Starting from the initial hole formation, through a short-time exponential regime, film retraction eventually reaches a constant speed, which can be shown to match the TC prediction when we account for the local changes in shape and in surface tension.

For systems at macroscopic length scale, growth of holes with an initial radius R0 < h0 is not energetically favored, i.e., an open hole closes.47 Hence, earlier experiments or numerical investigations18,29,48,49 on sheet retraction considered an initial hole dimension satisfying R0h0. Additionally, continuum assumptions require long-wave approximations to describe the rupture process, which is valid for h0 much larger than the atomic scale.50 However, in the limit of a few nanometers thickness (essentially Newtonian black films), the disjoining force acts as a destabilizer opposing the stabilizing effect of surface tension and, therefore, becomes responsible for the puncture of the film.51,52 At this scale, MD is the only appropriate tool to scrutinize the hole opening process. It was observed that R0h0 no longer remains a requirement for holes to expand (see the supplementary material, Sec. I B). This agrees with the seminal study of Koplik and Banavar,32 where the authors reported that for any dry patch larger than a few molecular diameters, the separated molecules suffer from mutual interactions, and the hole always grows.

The hole creation technique may invoke slight differences to the initial conditions, but any difference rapidly ceases.53 We confirm this by applying two separate techniques to create holes by (i) deleting atoms from a circular region of radius, R0, or (ii) applying a small radially decaying force for a short time that models poking a film54,55 (for details of these procedures, see the supplementary material, Sec. I). Henceforth, these two types will be referred to as cut- and poked-holes, respectively.

At the onset of retraction, an exponential growth for a short period of time is observed, which converges to a linear terminal speed at later times. Figure 2(a) portrays this short-term hole expansion for cut- and poked-holes. In agreement with the observations in Refs. 2 and 56, it is evident that the early-time hole growth shows an exponential behavior, i.e., R(t) − R0 = exp  (ct + d), where c and d are fitting parameters. While in Ref. 56, the short-term growth was ascribed to the viscoelastic effects, it was confirmed to be a generic feature of circular hole retraction 29 and was also accredited to the initial shape of the puncture.29,57

FIG. 2.

Retraction regimes and correction to Taylor–Culick speed.  (a) Short-term nonlinear retraction of the film for cut- and poked-holes with initial thickness h0 ≈ 14; time is non-dimensionalized by the inertial time scale, i.e., t*=t/ρh03/2γ. Faint and solid symbols denote, respectively, data from individual simulations and data averaged over ten simulations. The solid line is an exponential fit to the data of the form R(t) − R0 = ect+d. For the poked cases, a brief transitory behavior is seen before the exponential regime (left to the vertical dotted line). (b) Retraction velocity (blue: poked-holes, red: cut-holes) normalized by the TC speed U* = U/UTC with 2% standard deviation. Faint symbols denote cases where classical UTC is used, whereas solid symbols denote cases where UTC is corrected, considering the local surface tension (Uγ*) and the existence of a transition regime (Uγ,α*). ɛγ,T, ɛγ,κ, and ɛα denote, respectively, various corrections arising from the local variation of surface tension due to temperature and curvature and because of the existence of the transition regime. In (a) and (b), the non-linear growth region is shaded. (c) The schematic in the top panel shows an idealized retracting film. The mid-panel sheds molecular insights where a transition region is present between the rim and the unperturbed film. Green arrows show the directions normal (n) to the local film surface. The color map shows a radial-averaged instantaneous temperature field with velocity vectors (black arrows) overlaid. The bottom panel shows the local to global (time-averaged and radial-averaged) surface tension ratio (γ* = γr/γ) against distance r from the tip of the retracting film.

FIG. 2.

Retraction regimes and correction to Taylor–Culick speed.  (a) Short-term nonlinear retraction of the film for cut- and poked-holes with initial thickness h0 ≈ 14; time is non-dimensionalized by the inertial time scale, i.e., t*=t/ρh03/2γ. Faint and solid symbols denote, respectively, data from individual simulations and data averaged over ten simulations. The solid line is an exponential fit to the data of the form R(t) − R0 = ect+d. For the poked cases, a brief transitory behavior is seen before the exponential regime (left to the vertical dotted line). (b) Retraction velocity (blue: poked-holes, red: cut-holes) normalized by the TC speed U* = U/UTC with 2% standard deviation. Faint symbols denote cases where classical UTC is used, whereas solid symbols denote cases where UTC is corrected, considering the local surface tension (Uγ*) and the existence of a transition regime (Uγ,α*). ɛγ,T, ɛγ,κ, and ɛα denote, respectively, various corrections arising from the local variation of surface tension due to temperature and curvature and because of the existence of the transition regime. In (a) and (b), the non-linear growth region is shaded. (c) The schematic in the top panel shows an idealized retracting film. The mid-panel sheds molecular insights where a transition region is present between the rim and the unperturbed film. Green arrows show the directions normal (n) to the local film surface. The color map shows a radial-averaged instantaneous temperature field with velocity vectors (black arrows) overlaid. The bottom panel shows the local to global (time-averaged and radial-averaged) surface tension ratio (γ* = γr/γ) against distance r from the tip of the retracting film.

Close modal

Beyond the short-term exponential regime, the retraction velocity approaches a plateau—the constant TC speed.29 This is shown by the faint symbols in Fig. 2(b) for cut- (red) and poked-holes (blue). Retraction speed U is measured over 50 consecutive time steps and is non-dimensionalized by UTC from Eq. (1) with ϕ = 2, i.e., U* = U/UTC. Evident from Fig. 2(b), UTC is appreciably faster than the measured speed (U* ≈ 0.8), that is, MD predicts a speed that is only 80% of the Taylor–Culick prediction. This is not unexpected and documented in the literature.7,18,40,58,59 Notably, McEntee and Mysels18 studied the retraction of soap-films with a wide range of thicknesses. In their experiments, the retraction speeds for films of thickness smaller than 100 nm were found to be substantially slower than UTC, requiring ϕ ≈ 0.7. UTC is rather closer to our measured speed, requiring ϕ ≈ 1.35 ± 0.05 for h0 ≈ 4.75 nm. However, having the effects of surfactants and visco-elasticity45 completely eliminated in the present study (by considering a pure film), why the TC theory still fails to estimate the retraction rate is explored now.

A number of factors, i.e., temperature variation due to internal heat generation, curvature effects, and viscous contribution,24 might play a critical role in slowing down the retraction. Although viscosity does not appear in Eq. (1) and, therefore, does not affect the terminal speed, its relative importance to the surface tension decides the length over which the film is disturbed and how the momentum is distributed through the film.29 This is often characterized by the Ohnesorge number, Oh=μ/ργh0, which is moderate (0.8) for the present case.

Adjusting the value of ϕ marries the theoretical estimate with the measured speed, but it does not necessarily explain the underlying physics responsible for the slower retraction. This disagreement was attributed to the disturbance caused by the presence of a precursor ahead of the rim arising from a surface tension gradient due to surface shrinkage and surfactant concentration.18,60 It was inferred by Frankel and Mysels53 that the slow growth-rate should be due to a large variation of the local surface tension at the film-front. Although such local fluctuation is not unusual,61,62 capturing these with continuum models is difficult, and the only MD study is the quasi-2D case of Kono et al.39 who reported a considerable temperature rise (>5K) in the liquid surrounding the hole for h0 ≈ 5 nm.

In contrast to 2D systems, fully 3D simulations are computationally more expensive, but avoid the well-known deficiencies of 2D fluids, such as divergent transport coefficients.63 One expects the behavior of an interface to be dominated by the number of interacting neighboring molecules—which is much smaller in 2D. Therefore, a 2D interface would be overly susceptible to rupture.33 Even with periodic boundary conditions, the strong system-size dependency of the diffusion coefficient also cannot be ignored.64 The present study therefore minimizes these effects by simulating a fully three-dimensional (3D) system, which is 8 times larger (in terms of particle number) than the 2D system of Ref. 39. A warmer region surrounding the expanding hole is observed, and the temperature increase is rather less pronounced, 2.9 K vs >5 K; see the supplementary material, Sec. II B. Such an increase in temperature constitutes around a 5%–7% drop in γ. This is, of course, non-negligible, but is not sufficient to explain the slower retraction. Note that the local increase in temperature will slightly lower the viscosity of the liquid surrounding the hole, but should not result in a perceptible variation in Oh.

The retraction process is mainly driven by the tangential component of the pressure tensor,65,66 while the normal component is responsible for the radial growth of the rim. Resolving the pressure tensor spatially along the film profile in local tangential and normal directions (see the supplementary material, Sec. II B and Fig. 3) reveals a noticeable depression of the tangential pressure. From these constituent components, γ can be locally calculated using the Kirkwood–Buff definition,67 i.e., γr=12Pxx(r,x)Prr(r,x)dx, where r denotes the radial distance from the film tip and subscripts xx and rr denote the direct pressure along the x and r axis, respectively. A full rotation and mapping of the pressure tensor using a radially varying normal vector do not give a substantially different surface tension (see the supplementary material, Sec. II B). The tension acts as a “pulling force” in Taylor–Culick formulation, purely in the direction of retraction, seen from the work of Kirkwood and Buff with Prr being the main contribution to γ. The bottom-panel in Fig. 2(c) shows the ratio of the local to global surface tension (time-averaged and radial-averaged) with increasing r. Remarkably, this gives evidence of a substantial drop in γr (15%) at the tip of the retracting film (r = 0) and γrγ far from the rim. While the local temperature rise explains an 5%7% drop (ɛγ,T) in γ, the remaining can be attributed to surface curvature and shrinkage, ɛγ,κ, so ɛγ = ɛγ,T + ɛγ,κ. This is plausible that the latter reduction is linked to the Tolman length effect. However, given the complexity and ambiguity in quantifying the Tolman length,68,69 this is not pursued in the present work. Interested readers are referred to a previous work by two of the co-authors62 on the curvature dependence of surface tension for simpler systems.

Accounting for the local variation in surface tension, an improvement εγ=1γr=0/γ is obtained, and for the present case, using the surface tension averaged over the film-front in Eq. (1) results in agreement with the measured speed within only an 8% deviation, shown by the dashed line denoted by Uγ* in Fig. 2(b). These corrections, however, will diminish for thicker films as the variations of temperature39 and curvature become less pronounced (see the supplementary material, Sec. II B).

The classical formulation discounts any possibility of a transition regime between the rim and the undisturbed film. Our results confirm, for molecularly thin films, the existence of an intermediate (or transition) bridge between the rim and the undisturbed film [Fig. 1(b)]. The top panel in Fig. 2(c) shows an idealized picture of film retraction upon which the TC formalism is based on. This simplified description requires the liquid from the hole to be collected in a spherical rim bordering on the unperturbed film. Contrasting Taylor’s assumption of a sudden transformation, the present study rather agrees with the experimentally7,70 and numerically30,32 observed transition, shown in the mid-panel of Fig. 2(c). Previous studies reported an “aureole-shaped” regime for films of liquid with surfactants and attributed this to a surface tension gradient and surface elasticity.7,27,71 Existence of a similar “extended rim” in the present study—where the film is made of pure liquid—indicates a local variation of the surface force, but caused by curvature effects (geometry).

The transition can be explored through a mass balance argument by assuming an elliptic blob, of aspect ratio α = b/a, connected by a smooth transition region, of length L, to the neck of width h0. The length L is given by
(2)
where β = 2b/h0 (see the supplementary material, Sec. II C, for a full derivation). Transition must take place over the finite length, L,53which we confirm by comparing L to our simulations [see the supplementary material, Fig. 6(b)], thus disproving, for thinner films, the assumption of an unperturbed film bordering on the rim. As in Eq. (2), the transition length becomes more prominent with thinner films and disappears when bh0/2 corresponding to viscous film retraction.29,72
The MD data indicate that the normal pressure difference of the surface is minimal at the point where the rim meets the transition region. In order to satisfy the Young–Laplace equation, it follows that the absolute value of the curvature must also be minimized. Under these assumptions, the atomistic form of Eq. (1) can be expressed as
(3)
For simplicity, assuming a spherical rim, i.e., α = 1, yields a surface-shape correction factor, ɛα = 2−1/4 ≈ 0.84. When ɛα is accounted for, in addition to ɛγ, the corrected TC model, (Uγ,α*) shows close agreement with the MD results. This is shown by the solid symbols in Fig. 2(b), and as evidenced in Fig. 2(b), Uγ,α*1.

There remain a few other aspects that plausibly affect the dynamics of film retraction. Any mutual dependency of the two correction factors developed in this study cannot be ruled out. The growth of the blob is non-negligible as compared to Ṙ. In situations like this, it has been argued24,27,29 that the velocity within the blob is non-uniform; the velocity vectors overlaid on the temperature field in Fig. 2(c) confirm this hypothesis. The effect of this velocity would contribute to the departure of the blob shape from being circular and, therefore, is assumed to be included in the surface shape correction, ɛα. The effect of evaporation, however, in agreement with the previous study,39 was found to be negligible.

To summarize, this study examines the Taylor–Culick theory for film retraction at atomic scale using large-scale MD simulations. A similar growth rate of holes made by two separate methods, either poking or cutting, confirms that the dynamics is independent of hole creation mechanism. This then allows for a systematic study, where we repeatedly show a short term exponential growth period followed by a constant velocity regime over an ensemble of cases—both behaviors observed in the experimental literature. A Lennard-Jones model is used to remove viscoelastic effects, so various driving contributions to the growth dynamics can be exposed. The theoretical conjecture of a constant terminal retraction rate is accurate at this scale; however, we show that estimating the “true” retraction speed requires corrections for (i) the local variation of surface tension arising from surface curvature and temperature rise and (ii) the existence of a transition regime between the rim and the unperturbed film—both of these factors were ignored in the classical formulation. Surface tension and interface shapes are outputs of MD simulations, and we observe how they contribute to the dynamics and quantify their impact. While their contributions are marginal for macroscopic films, these are the key reasons for the failure of the classical theory as the film thickness is reduced to nano-scale.

This work provides compelling evidence that the pioneering theory of Taylor and Culick, which in its original form fails at the atomic scale, can be corrected to accurately predict the film retraction speed.

See the supplementary material for further details on the simulation methodology, the derivation of the classical theory, the measurement of the local surface tension, and the derivation for the correction factor arising from a transition regime.

M.R.R. acknowledges Shell and the Beit Trust for Ph.D. funding through a Beit Fellowship for Scientific Research. L.S. acknowledges the Engineering and Physical Sciences Research Council (EPSRC) for a Postdoctoral Fellowship (Grant No. EP/V005073/1). J.P.E. was supported by the Royal Academy of Engineering (RAEng) through their Research Fellowships scheme. B.C. was supported by Shell and the EPSRC via an iCASE Ph.D. studentship (Grant No. EP/T517690/1). D.D. acknowledges a Shell/RAEng Research Chair in Complex Engineering Interfaces and the EPSRC for an Established Career Fellowship (Grant No. EP/N025954/1).

The authors have no conflicts to disclose.

Muhammad Rizwanur Rahman: Formal analysis (lead); Investigation (lead); Methodology (equal); Writing – original draft (lead); Writing – review & editing (equal). Li Shen: Formal analysis (supporting); Funding acquisition (supporting); Supervision (supporting); Writing – review & editing (equal). James P. Ewen: Methodology (supporting); Supervision (supporting); Writing – review & editing (equal). Benjamin Collard: Formal analysis (supporting). D. M. Heyes: Supervision (supporting); Writing – review & editing (equal). Daniele Dini: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Supervision (equal); Writing – review & editing (equal). E. R. Smith: Conceptualization (equal); Methodology (equal); Software (lead); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Codes to reproduce and analyze the data reported in this work can be found at https://github.com/MuhammadRRahman/Thin-Film-Rupture-NEMD.git.

1.
U.
Thiele
,
M. G.
Velarde
, and
K.
Neuffer
, “
Dewetting: Film rupture by nucleation in the spinodal regime
,”
Phys. Rev. Lett.
87
,
016104
(
2001
).
2.
J.-L.
Masson
and
P. F.
Green
, “
Hole formation in thin polymer films: A two-stage process
,”
Phys. Rev. Lett.
88
,
205504
(
2002
).
3.
P.-G.
De Gennes
,
F.
Brochard-Wyart
,
D.
Quéré
et al,
Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves
(
Springer
,
2004
), Vol. 315.
4.
P. G.
de Gennes
,
Bubbles, Foams and Other Fragile Objects
(
Springer Netherlands
,
Dordrecht
,
2008
), pp.
77
85
.
5.
J. C.
Bird
,
R.
De Ruiter
,
L.
Courbin
, and
H. A.
Stone
, “
Daughter bubble cascades produced by folding of ruptured thin films
,”
Nature
465
,
759
762
(
2010
).
6.
J.
Bico
, “
Cracks in bursting soap films
,”
J. Fluid Mech.
778
,
1
4
(
2015
).
7.
P. C.
Petit
,
M.
Le Merrer
, and
A.-L.
Biance
, “
Holes and cracks in rigid foam films
,”
J. Fluid Mech.
774
,
R3
(
2015
).
8.
S.
Poulain
,
E.
Villermaux
, and
L.
Bourouiba
, “
Ageing and burst of surface bubbles
,”
J. Fluid Mech.
851
,
636
671
(
2018
).
9.
D.
Langevin
, “
On the rupture of thin films made from aqueous surfactant solutions
,”
Adv. Colloid Interface Sci.
275
,
102075
(
2020
).
10.
L.
Bourouiba
, “
The fluid dynamics of disease transmission
,”
Annu. Rev. Fluid Mech.
53
,
473
508
(
2021
).
11.
J.
Wu
, “
Evidence of sea spray produced by bursting bubbles
,”
Science
212
,
324
326
(
1981
).
12.
F.
Veron
, “
Ocean spray
,”
Annu. Rev. Fluid Mech.
47
,
507
538
(
2015
).
13.
C.
Marangoni
,
P.
Stefanelli
, and
R.
Liceo
, “
Monografia sulle bolle liquide
,”
Il Nuovo Cimento (1869-1876)
7
,
301
356
(
1872
).
14.
W. E.
Ranz
, “
Some experiments on the dynamics of liquid films
,”
J. Appl. Phys.
30
,
1950
1955
(
1959
).
15.
E.
Chatzigiannakis
and
J.
Vermant
, “
Breakup of thin liquid films: From stochastic to deterministic
,”
Phys. Rev. Lett.
125
,
158001
(
2020
).
16.
G. I.
Taylor
, “
The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets
,”
Proc. R. Soc. London, Ser. A
253
,
313
321
(
1959
).
17.
F. E. C.
Culick
, “
Comments on a ruptured soap film
,”
J. Appl. Phys.
31
,
1128
1129
(
1960
).
18.
W. R.
McEntee
and
K. J.
Mysels
, “
Bursting of soap films. I. An experimental study
,”
J. Phys. Chem.
73
,
3018
3028
(
1969
).
19.
A.
Dupré
,
Théorie Mécanique de la Chaleur
(
Gauthier-Villars
,
1869
).
20.
L. J.
Evers
,
S. Y.
Shulepov
, and
G.
Frens
, “
Rupture of thin liquid films from Newtonian and viscoelastic liquids bursting behaviour of Newton-black films
,”
Faraday Discuss.
104
,
335
344
(
1996
).
21.
A. T.
Florence
and
K. J.
Mysels
, “
Bursting of soap films. VI. Effect of surfactant purity
,”
J. Phys. Chem.
78
,
234
235
(
1974
).
22.
L. J.
Evers
,
S. Y.
Shulepov
, and
G.
Frens
, “
Bursting dynamics of thin free liquid films from Newtonian and viscoelastic solutions
,”
Phys. Rev. Lett.
79
,
4850
(
1997
).
23.
C.
Planchette
,
F.
Marangon
,
W.-K.
Hsiao
, and
G.
Brenn
, “
Breakup of asymmetric liquid ligaments
,”
Phys. Rev. Fluids
4
,
124004
(
2019
).
24.
J.-L.
Pierson
,
J.
Magnaudet
,
E. J.
Soares
, and
S.
Popinet
, “
Revisiting the Taylor-Culick approximation: Retraction of an axisymmetric filament
,”
Phys. Rev. Fluids
5
,
073602
(
2020
).
25.
D.
Lohse
and
E.
Villermaux
, “
Double threshold behavior for breakup of liquid sheets
,”
Proc. Natl. Acad. Sci. U. S. A.
117
,
18912
18914
(
2020
).
26.
J. B.
Keller
, “
Breaking of liquid films and threads
,”
Phys. Fluids
26
,
3451
3453
(
1983
).
27.
M.
De Corato
,
D.
Tammaro
,
P. L.
Maffettone
, and
N.
Fueyo
, “
Retraction of thin films coated by insoluble surfactants
,”
J. Fluid Mech.
942
,
A55
(
2022
).
28.
C. R.
Constante-Amores
et al, “
Role of surfactant-induced Marangoni stresses in retracting liquid sheets
,”
J. Fluid Mech.
949
,
A32
(
2022
).
29.
N.
Savva
and
J. W. M.
Bush
, “
Viscous sheet retraction
,”
J. Fluid Mech.
626
,
211
240
(
2009
).
30.
L.
Gordillo
,
G.
Agbaglah
,
L.
Duchemin
, and
C.
Josserand
, “
Asymptotic behavior of a retracting two-dimensional fluid sheet
,”
Phys. Fluids
23
,
122101
(
2011
).
31.
H.
Deka
and
J. L.
Pierson
, “
Revisiting the Taylor-Culick approximation. II. Retraction of a viscous sheet
,”
Phys. Rev. Fluids
5
,
093603
(
2020
).
32.
J.
Koplik
and
J. R.
Banavar
, “
Molecular simulations of dewetting
,”
Phys. Rev. Lett.
84
,
4401
(
2000
).
33.
J.
Koplik
and
J. R.
Banavar
, “
Molecular dynamics of interface rupture
,”
Phys. Fluids A
5
,
521
536
(
1993
).
34.
C.
Zhao
,
J. E.
Sprittles
, and
D. A.
Lockerby
, “
Revisiting the Rayleigh–Plateau instability for the nanoscale
,”
J. Fluid Mech.
861
,
R3
(
2019
).
35.
J.
Diemand
,
R.
Angélil
,
K. K.
Tanaka
, and
H.
Tanaka
, “
Direct simulations of homogeneous bubble nucleation: Agreement with classical nucleation theory and no local hot spots
,”
Phys. Rev. E
90
,
052407
(
2014
).
36.
K.
Kadau
,
J. L.
Barber
,
T. C.
Germann
,
B. L.
Holian
, and
B. J.
Alder
, “
Atomistic methods in fluid simulation
,”
Proc. R. Soc. London, Ser. A
368
,
1547
1560
(
2010
).
37.
B. D.
Todd
and
P. J.
Daivis
,
Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications
(
Cambridge University Press
,
2017
).
38.
Y.
Zhang
,
J. E.
Sprittles
, and
D. A.
Lockerby
, “
Molecular simulation of thin liquid films: Thermal fluctuations and instability
,”
Phys. Rev. E
100
,
023108
(
2019
).
39.
S.
Kono
,
T.
Kaneko
, and
I.
Ueno
, “
Elevation of the temperature of liquid films caused by rapid rupturing
,”
Phys. Rev. E
90
,
051004
(
2014
).
40.
D.
Tammaro
et al, “
Flowering in bursting bubbles with viscoelastic interfaces
,”
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2105058118
(
2021
).
41.
H.
Manikantan
and
T. M.
Squires
, “
Surfactant dynamics: Hidden variables controlling fluid flows
,”
J. Fluid Mech.
892
,
P1
(
2020
).
42.
E.
Smith
, “
On the coupling of molecular dynamics to continuum computational fluid dynamics
,” Ph.D. thesis,
Imperial College London
,
2013
.
43.
D.
Frenkel
,
B.
Smit
, and
M. A.
Ratner
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Academic Press San Diego
,
1996
), Vol. 2.
44.
A.
Vrij
, “
Possible mechanism for the spontaneous rupture of thin, free liquid films
,”
Discuss. Faraday Soc.
42
,
23
33
(
1966
).
45.
A.
Vrij
and
J. T. G.
Overbeek
, “
Rupture of thin liquid films due to spontaneous fluctuations in thickness
,”
J. Am. Chem. Soc.
90
,
3074
3078
(
1968
).
46.
I. B.
Ivanov
,
B.
Radoev
,
E.
Manev
, and
A.
Scheludko
, “
Theory of the critical thickness of rupture of thin liquid films
,”
Trans. Faraday Soc.
66
,
1262
1273
(
1970
).
47.
G. I.
Taylor
and
D. H.
Michael
, “
On making holes in a sheet of fluid
,”
J. Fluid Mech.
58
,
625
639
(
1973
).
48.
A. B.
Pandit
and
J. F.
Davidson
, “
Hydrodynamics of the rupture of thin liquid films
,”
J. Fluid Mech.
212
,
11
24
(
1990
).
49.
V.
Sanjay
,
U.
Sen
,
P.
Kant
, and
D.
Lohse
, “
Taylor–Culick retractions and the influence of the surroundings
,”
J. Fluid Mech.
948
,
A14
(
2022
).
50.
D.
Vaynblat
,
J. R.
Lister
, and
T. P.
Witelski
, “
Rupture of thin viscous films by van der Waals forces: Evolution and self-similarity
,”
Phys. Fluids
13
,
1130
1140
(
2001
).
51.
B.
Néel
and
E.
Villermaux
, “
The spontaneous puncture of thick liquid films
,”
J. Fluid Mech.
838
,
192
221
(
2018
).
52.
E.
Villermaux
, “
Fragmentation versus cohesion
,”
J. Fluid Mech.
898
,
P1
(
2020
).
53.
S.
Frankel
and
K. J.
Mysels
, “
Bursting of soap films. II. Theoretical considerations
,”
J. Phys. Chem.
73
,
3028
3038
(
1969
).
54.
B.
Mandracchia
et al, “
Quantitative imaging of the complexity in liquid bubbles’ evolution reveals the dynamics of film retraction
,”
Light: Sci. Appl.
8
,
20
(
2019
).
55.
A. T.
Oratis
,
J. W. M.
Bush
,
H. A.
Stone
, and
J. C.
Bird
, “
A new wrinkle on liquid sheets: Turning the mechanism of viscous bubble collapse upside down
,”
Science
369
,
685
688
(
2020
).
56.
G.
Debrégeas
,
P.
Martin
, and
F.
Brochard-Wyart
, “
Viscous bursting of suspended films
,”
Phys. Rev. Lett.
75
,
3886
(
1995
).
57.
C.
Roth
,
B.
Deh
,
B.
Nickel
, and
J.
Dutcher
, “
Evidence of convective constraint release during hole growth in freely standing polystyrene films at low temperatures
,”
Phys. Rev. E
72
,
021802
(
2005
).
58.
F.
Müller
,
U.
Kornek
, and
R.
Stannarius
, “
Experimental study of the bursting of inviscid bubbles
,”
Phys. Rev. E
75
,
065302
(
2007
).
59.
T.
Trittel
,
T.
John
,
K.
Tsuji
, and
R.
Stannarius
, “
Rim instability of bursting thin smectic films
,”
Phys. Fluids
25
,
052106
(
2013
).
60.
K. J.
Mysels
and
J. A.
Stikeleather
, “
The bursting of soap films. III. Spontaneous and induced bursting
,”
J. Colloid Interface Sci.
35
,
159
162
(
1971
).
61.
M. P.
Moody
and
P.
Attard
, “
Curvature dependent surface tension from a simulation of a cavity in a Lennard-Jones liquid close to coexistence
,”
J. Chem. Phys.
115
,
8967
8977
(
2001
).
62.
J.
Wen
,
D.
Dini
,
H.
Hu
, and
E. R.
Smith
, “
Molecular droplets vs bubbles: Effect of curvature on surface tension and Tolman length
,”
Phys. Fluids
33
,
072012
(
2021
).
63.
J. H.
Irving
and
J. G.
Kirkwood
, “
The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics
,”
J. Chem. Phys.
18
,
817
(
1950
).
64.
I.-C.
Yeh
and
G.
Hummer
, “
System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions
,”
J. Phys. Chem. B
108
,
15873
15879
(
2004
).
65.
C.
Braga
,
E. R.
Smith
,
A.
Nold
,
D. N.
Sibley
, and
S.
Kalliadasis
, “
The pressure tensor across a liquid-vapour interface
,”
J. Chem. Phys.
149
,
044705
(
2018
).
66.
K.
Shi
,
E. R.
Smith
,
E. E.
Santiso
, and
K. E.
Gubbins
, “
A perspective on the microscopic pressure (stress) tensor: History, current understanding, and future challenges
,”
J. Chem. Phys.
158
,
040901
(
2023
).
67.
J. G.
Kirkwood
and
F. P.
Buff
, “
The statistical mechanical theory of surface tension
,”
J. Chem. Phys.
17
,
338
343
(
1949
).
68.
Y. A.
Lei
,
T.
Bykov
,
S.
Yoo
, and
X. C.
Zeng
, “
The Tolman length: Is it positive or negative?
,”
J. Am. Chem. Soc.
127
,
15346
15347
(
2005
).
69.
A. E.
van Giessen
and
E. M.
Blokhuis
, “
Direct determination of the Tolman length from the bulk pressures of liquid drops via molecular dynamics simulations
,”
J. Chem. Phys.
131
,
164705
(
2009
).
70.
H.
Lhuissier
and
E.
Villermaux
, “
Soap films burst like flapping flags
,”
Phys. Rev. Lett.
103
,
054501
(
2009
).
71.
H.
Lhuissier
and
E.
Villermaux
, “
Destabilization of flapping sheets: The surprising analogue of soap films
,”
C. R. Mec.
337
,
469
480
(
2009
).
72.
M. P.
Brenner
and
D.
Gueyffier
, “
On the bursting of viscous films
,”
Phys. Fluids
11
,
737
739
(
1999
).

Supplementary Material