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.
I. INTRODUCTION
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.
II. METHODS
A. Simulation setup
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 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 ( × 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).
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.
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.
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).
III. RESULTS AND DISCUSSIONS
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.
A. Hole formation
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 R0 ≫ h0. 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 R0 ≫ h0 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.
B. Short-term retraction
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
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., . 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 and the existence of a transition regime . ɛγ,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.
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., . 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 and the existence of a transition regime . ɛγ,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.
C. Taylor–Culick retraction
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, , which is moderate for the present case.
D. Temperature effects on surface tension
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 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 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 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.
E. Curvature effects on surface tension
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., , 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 at the tip of the retracting film (r = 0) and γr → γ far from the rim. While the local temperature rise explains an 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 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 deviation, shown by the dashed line denoted by 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).
F. Transition regime
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).
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.
IV. CONCLUSION
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.
SUPPLEMENTARY MATERIAL
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.
ACKNOWLEDGMENTS
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).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
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).
DATA AVAILABILITY
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.