Taylor-Culick retractions and the influence of the surroundings

When a freely suspended liquid film ruptures, it retracts spontaneously under the action of surface tension. If the film is surrounded by air, the retraction velocity is known to approach the constant Taylor-Culick velocity. However, when surrounded by an external viscous medium, the dissipation within that medium dictates the magnitude of the retraction velocity. In the present work, we study the retraction of a liquid (water) film in a viscous oil ambient (\emph{two-phase} Taylor-Culick retractions), and that sandwiched between air and a viscous oil (\emph{three-phase} Taylor-Culick retractions). In the latter case, the experimentally-measured retraction velocity is observed to have a weaker dependence on the viscosity of the oil phase as compared to the configuration where the water film is surrounded completely by oil. Numerical simulations indicate that this weaker dependence arises from the localization of viscous dissipation near the three-phase contact line. The speed of retraction only depends on the viscosity of the surrounding medium and not on that of the film. From the experiments and the numerical simulations, we reveal unprecedented regimes for the scaling of the Weber number $We_f$ of the film (based on its retraction velocity) or the capillary number $Ca_s$ of the surroundings vs. the Ohnesorge number $Oh_s$ of the surroundings in the regime of large viscosity of the surroundings ($Oh_s \gg 1$), namely $We_f \sim Oh_s^{-2}$ and $Ca_s \sim Oh_s^{0}$ for the two-phase Taylor-Culick configuration, and $We_f \sim Oh_s^{-1}$ and $Ca_s \sim Oh_s^{1/2}$ for the three-phase Taylor-Culick configuration.

When a freely suspended liquid film ruptures, it retracts spontaneously under the action of surface tension. If the film is surrounded by air, the retraction velocity is known to approach the constant Taylor-Culick velocity. However, when surrounded by an external viscous medium, the dissipation within that medium dictates the magnitude of the retraction velocity. In the present work, we study the retraction of a liquid (water) film in a viscous oil ambient (two-phase Taylor-Culick retractions), and that sandwiched between air and a viscous oil (three-phase Taylor-Culick retractions). In the latter case, the experimentally-measured retraction velocity is observed to have a weaker dependence on the viscosity of the oil phase as compared to the configuration where the water film is surrounded completely by oil. Numerical simulations indicate that this weaker dependence arises from the localization of viscous dissipation near the three-phase contact line. The speed of retraction only depends on the viscosity of the surrounding medium and not on that of the film. From the experiments and the numerical simulations, we reveal unprecedented regimes for the scaling of the Weber number We f of the film (based on its retraction velocity) or the capillary number Ca s of the surroundings vs. the Ohnesorge number Oh s of the surroundings in the regime of large viscosity of the surroundings (Oh s 1), namely We f ∼ Oh −2 s and Ca s ∼ Oh 0 s for the two-phase Taylor-Culick configuration, and We f ∼ Oh −1

Introduction
Liquid films, sheets, and shells have piqued the interest of fluid dynamicists for nearly two centuries (Savart 1833a,b,c;Taylor 1959a,b;Clanet 2007;Villermaux 2020). A freely suspended liquid film warrants further attention among these configurations as it is inherently metastable owing to its high surface area. Indeed, if a large enough † Both authors contributed equally to this work.
‡ Email address for correspondence: vatsalsanjay@gmail.com ¶ Email address for correspondence: u.sen@utwente.nl Email address for correspondence: d.lohse@utwente.nl of water film (f ) in an oil (s) environment (two-phase configuration), (c) retraction of a water film sandwiched between air and oil (three-phase configuration). The dot-dashed line represents the axis of rotational symmetry, and R(t) is the radius of the growing hole centered at this axis. In all the schematics, the water film is retracting from left to right with velocity v f , as indicated by the arrow, and γij denotes the surface tension coefficient between fluids i and j. (Taylor & Michael 1973;Villermaux 2020) hole nucleates on the film, the sheet will spontaneously retract to reduce its surface area. Such interfacial destabilization leading to film rupture and bursting can also result in waterborne disease transmission (Bourouiba 2021). The bursting of liquid films at an oil-air interface is important for various industrial applications in the chemical and petrochemical engineering sectors as well. One area of particular interest is underwater oil spills in oceans, such as the Deepwater Horizon spill in 2010 in the Gulf of Mexico (Summerhayes 2011). For these spills, droplets (or slugs) of oil may rise to the free surface of water via buoyancy, and then rupture the free surface of water directly above it. The water film will retract upon rupture, and the oil will spread on the water surface, thus perpetuating an environmental hazard. Perhaps the most widely studied example of sheet destabilization and retraction is during the bursting of liquid (e.g., soap) films in air -an area of active research since the pioneering works (Dupré 1867(Dupré , 1869Rayleigh 1891;Taylor 1959b;Culick 1960;McEntee & Mysels 1969) in the late nineteenth and mid-twentieth century to the more recent investigations (Bremond & Villermaux 2005;Müller et al. 2007;Lhuissier & Villermaux 2012;Munro et al. 2015;Deka & Pierson 2020). In these studies, the outer medium is assumed passive (inviscid and zero-inertia). The origin of the nucleation of the initial hole in the film can be manifold (Lohse & Villermaux 2020). After film rupture, the internal viscous stresses in the film do not contribute to the momentum balance, but dictate the distribution of momentum within the film (Savva & Bush 2009), as long as the Ohnesorge number of the film (ratio of its visco-capillary to inertio-capillary time scales, see § 4) is less than its aspect ratio (see Deka & Pierson 2020). Nonetheless, half of the surface energy released goes into internal viscous dissipation (see appendix A and Culick (1960);de Gennes (1996); Sünderhauf et al. (2002); Villermaux (2020)).
A representative schematic of the situation mentioned above is shown in figure 1a (henceforth referred to as the classical Taylor-Culick configuration), where the water film (f ) of thickness h 0 is retracting in air (a) under the action of surface tension. The retraction velocity, v f , in such a scenario is constant (after a period of initial transience) and approaches the Taylor-Culick velocity given by where 2γ af is the net surface tension driving the retraction (γ af being the interfacial tension coefficient between the film and air) and ρ f is the density of the liquid film (figure 1a). Furthermore, it was observed that during the retraction, the liquid collects in a thicker rim at the retracting edge of the film, particularly for low viscosity liquids (Rayleigh 1891;Ranz 1959;Pandit & Davidson 1990, not depicted in figure 1a). The seminal work of Keller (1983) further explores the retraction of these films of non-uniform thickness.
The effect of viscosity of the film (η f ) during its retraction process has also been studied (Debrégeas et al. 1995(Debrégeas et al. , 1998. Brenner & Gueyffier (1999) showed that although viscosity does not have any effect on the constant retraction velocity, it can have a significant effect on the shape of the retracting edge of a planar film. They report that if the radial extent of the film is greater than its Stokes length (= η f /(ρ f v TC )), a growing rim is formed, whereas the rim is absent for the converse situation. Savva & Bush (2009) extended the work by Brenner & Gueyffier (1999) for highly viscous films, and also developed a lubrication model for the retraction dynamics of a circular hole. They concluded that although viscosity does not determine the magnitude of the constant retraction velocity, it does dictate the time required (post rupture) to attain that constant velocity, which increases with increasing viscosity. Recently, Pierson et al. (2020); Deka & Pierson (2020) revisited the viscous retraction dynamics by exploring self-similar solutions for slender filaments and sheets of finite length.
In all the aforementioned studies, the surrounding medium is assumed to play no role in the rupture dynamics. A question naturally arises: what happens when the outer medium also interacts with the retracting film? In particular, how do the viscosity and the inertia of the outer medium influence the rupture dynamics (Mysels & Vijayendran 1973;Joanny & de Gennes 1987;Reyssat & Quéré 2006;Jian et al. 2020b)? A representative schematic for such a scenario is shown in figure 1b, where a water film is retracting in a viscous oil ambient. This geometry will henceforth be referred to as the two-phase configuration where the net surface tension force responsible for retraction is 2γ sf (γ sf being the interfacial tension coefficient between the film and the surrounding medium, see figure 1b). In such a situation, viscous dissipation is not limited only to the retracting film, but is also present in the ambient. If the ambient happens to be significantly more viscous than the film, then the dissipation in the ambient dominates. In such a situation, the retraction velocity is still a constant. However, unlike the classical case, the velocity depends on the viscosity η s of the ambient medium (Martin et al. 1994;Reyssat & Quéré 2006). Common realizations of this configuration include relaxation of filaments and droplets in a viscous medium (Stone & Leal 1989), or that of air-films during drop impact (Jian et al. 2020a,b). Additionally, in this context, Anthony et al. (2020) showed that the so-called "inertially limited viscous regime" in the early-times of drop coalescence (Paulsen et al. 2012;Paulsen 2013) stems from a Taylor-Culick type retraction of the air film between the deformable drops.
In the present work, we study the influence of the surrounding medium on the retraction velocity of a ruptured liquid film using both force balance and energy conservation arguments. To accomplish this goal, along with the two canonical configurations shown in figures 1a and 1b, we also study the retraction dynamics of a liquid film sandwiched between air and a viscous oil bath. A representative schematic is shown in figure 1c. This geometry will henceforth be referred to as the three-phase configuration. This paper elucidates this case experimentally by inflating an oil droplet at the water-air interface and letting the water film rupture. Such a configuration can also be found in the early stages of water film retraction when an air bubble approaches a water-oil interface if the oil-layer is thick enough (Feng et al. 2014(Feng et al. , 2016. Furthermore, we also use direct numerical simulations (DNS) to demystify the retraction dynamics by using a precursor film-based three-fluid volume of fluid (VoF) method. We show that the film in this three-phase configuration still retracts with a constant velocity, and similar to the twophase case, the retraction velocity depends on the viscosity η s of the oil bath. However, this dependence is weaker in the three-phase configuration. Furthermore, we reveal an unprecedented scaling relationship for the retraction velocity of the film, which arises from the localization of the viscous dissipation near the three-phase contact line. The paper is organized as follows: § 2 describes the problem statement for the threephase Taylor-Culick retractions along with the experimental method employed to probe this configuration. The results from these experiments are discussed in § 3. § 4 presents the numerical framework, and § 5 describes the simulation results for both the twophase and three-phase configurations. § 6 demonstrates the balance of forces in Taylor-Culick retractions, followed by the corresponding scaling relationships in § 7. Further, § 8 analyzes the overall energy balance, highlighting the differences in the viscous dissipation mechanisms between the two-phase and three-phase configurations. The work culminates with conclusions in § 9. Throughout the manuscript, we refer to Appendix A for discussions on the classical Taylor-Culick retractions, and use the experimental datapoints from Reyssat & Quéré (2006) for the two-phase configuration.

Film bursting at an air-liquid interface: experimental method
We study the three-phase configuration experimentally by inflating an oil drop ('s' for 'surroundings' in figure 1c) at a water-air free interface and capturing the retraction of the water film (f in figure 1c). The schematic of the experimental setup is shown in figure 2a. A plastic box of dimensions 25 mm × 25 mm × 15 mm (length × width × height, Bodemschat) filled with purified water (Milli-Q) was used as the liquid bath for most of the experiments. To study the effect of the viscosity of the retracting film, the water in the bath was replaced by glycerol (Sigma-Aldrich)-water mixtures (concentrations in the range 50% -70% by wt.) for some experiments. A dispensing needle (inner diameter = 0.41 mm, HSW Fine-Ject) was submerged within the bath such that its dispensing end was at a depth of 2.4 mm from the free surface (depth kept constant during all experiments). A silicone oil (Wacker) droplet was created at the tip of the needle by connecting it to an oil-filled plastic syringe (5 mL, Braun Injekt) via a flexible plastic PEEK tubing (Upchurch Scientific). The oil flow rate was maintained at 0.05 mL/min with the help of a syringe infusion pump (Harvard Apparatus). In the present experiments, silicone oils of different viscosities were used, and their densities (ρ s ), kinematic (ν s ), and dynamic viscosities (η s ) are listed in Table 1. It is to be noted that for only the AK 0.65 oil (η s = 4.94 × 10 −4 Pa.s), the drop is less viscous than the water film (η f = 8.9 × 10 −4 Pa.s), while for all the other oils, the film is less viscous. The oil-water interfacial tension (γ sf ) was considered to be 0.040 N/m (Peters & Arabali 2013). The drop volume was increased by a slow infusion using the syringe pump. The needle depth below the free surface was chosen such that the drop remained anchored to the needle during inflation. As a result, the water film right above the oil droplet progressively thins with increasing volume of the drop. Below a certain thickness, the film ruptures due to van der Waals forces (Vaynblat et al. 2001), and subsequently retracts into the bath. This situation is analogous to the rupture and retraction of a water film sandwiched between air and a viscous oil droplet. This is also a configuration that is flipped vertically as compared to the early time scenario studied by Feng et al. (2016). Another key difference is that we ensure negligible vertical velocity at the point of rupture of the water film, making this scenario ideal for studying three-phase Taylor Culick retractions. High-speed imaging of this rupture and retraction phenomena was performed at 50000 fps (frames-per-second) for the lower viscosity oils and at 10000 fps for the higher viscosity ones, with a 2.5 µs exposure time, by a high-speed camera (Fastcam Nova S12, Photron) connected to a macro lens (DG Macro 105 mm, Sigma) with 64 mm of lens extender (Kenko). The camera was pointed at a plane mirror (Thorlabs) inclined at 45 • to the horizontal to capture the top view of the retraction phenomenon (figure 2a), while the experiments were illuminated from the top by a LED light source (KL 2500 LED, Schott). A typical bursting event is shown in figure 2b, where time t = t 0 indicates the instant when rupture is optically discernible. With increasing time, the size of the hole formed due to rupture increases as the film retracts. The phenomenological observations shown in figure 2b will be discussed in detail in § 3. The captured images were then further analyzed using the open-source software FIJI (Schindelin et al. 2012) and an in-house OpenCV-based Python script to obtain quantitative information presented in the following sections.

Film bursting at an air-liquid interface: experimental results
The rupture and retraction of a water film on the surface of an oil drop of η s = 4.94 × 10 −4 Pa.s is shown in figure 3a (and supplementary movie SM1). The timestamps indicate (t − t 0 ), where t 0 is the time instant when rupture is optically discernible, and t is the current time. As mentioned earlier, in this particular case, the water film is more viscous than the oil. It is to be noted that in the present experiments, we could not precisely control the location of rupture as it was sensitive to experimental noise (see § 4.2 of Villermaux 2020). Hence, the rupture in the present experiments did not always occur at the apex of the thinning film. Such behavior was also observed in other similar experiments of film rupture (Oldenziel et al. 2012;de Maleprade et al. 2016). The rupture location may also be determined by a 'prehole' formation, also observed by Vernay et al. (2015) for the bursting of emulsion-based liquid sheets. In their work, Vernay et al. (2015) show that the presence of emulsion oil droplets at the air-water interface results in lowering of the local interfacial tension, leading to Marangoni flows away from that location. This flow leads to a local thinning of the film, which ultimately ruptures at that location. They also report that the prehole formation always precedes rupture in their experiments. In the present experiments, the water surface is never pristine and always contains small impurities (which are practically unavoidable). It is possible that these impurities might have reduced the local surface tension, resulting in a similar Marangoni flow leading to a prehole. For the discussion on the origin of the hole nucleation, we also refer to Lohse & Villermaux (2020). In any case, upon rupture, a circular hole is formed in the film, which grows radially in time. Therefore, the oil bounded by the periphery of the hole gets into contact with air and not with the water film. It is also noticeable that the edge of the retracting film forms a thick rim -an observation also made for the retraction of liquid films in air (Pandit & Davidson 1990;Brenner & Gueyffier 1999;Sünderhauf et al. 2002). The presence of the rim can be qualitatively surmised from the experimental snapshots (figures 2b, 3a, and 3b), where the change in the curvature of the film downstream of the rim introduces a difference in the color intensity. As the hole increases in size (or as the film retracts further), this rim also becomes thicker. Finally, since silicone oil prefers to spread on water (Li et al. 2020), the retraction process ceases when the oil droplet has completely spread on water, thus creating a macroscopic film whose thickness is controlled by volume conservation and thermodynamics (de Gennes et al. 2004). When the viscosity of the oil phase is increased to η s = 9.30 × 10 −3 Pa.s, the hole opening (or film retraction) dynamics (as seen in figure 3b and supplementary movie SM1) are qualitatively similar to that for the lower viscosity described before (figure 3a). Here also, the film forms a thick rim at its retracting edge. Nonetheless, an increase in the oil's viscosity decreases the retraction speed, as indicated by the timestamps (corresponding to t − t 0 ). This behavior is expected since the physical situation is analogous to a retracting water film shearing the free surface of viscous oil: increasing η s increases the resistance to shearing, which in turn makes the retraction process slower.
Furthermore, the experimental snapshots show that the oil-air-water contact line exhibits corrugations during the retraction process, and fine streams of droplets are released from these corrugations. Similar observations were also made for film retraction in the two-phase configuration (Reyssat & Quéré 2006;Oldenziel et al. 2012) and during the rupture of the intermediate film when a drop coalesces with a pool of the same liquid in the presence of an external medium (Aryafar & Kavehpour 2008;Kavehpour 2015). The nature of the corrugations is reminiscent of the sharp tips observed during selective withdrawal (Cohen & Nagel 2002;Courrech du Pont & Eggers 2006 or tip streaming (Montanero & Gañán Calvo 2020). In a frame of reference co-moving with the rim, the film sees highly viscous oil being aspirated away from it, resulting in the formation of the sharp tips. Indeed, such a mechanism was also hinted at by Reyssat & Quéré (2006) for the instabilities observed in their experiments for film retraction in the two-phase configuration. Tseng & Prosperetti (2015) showed that such instabilities are formed due to the local convergence of streamlines in the neighborhood of a zero-vorticity point or line on the interface. However, a detailed and quantitative investigation of the formation and subsequent breakup of these liquid tips is beyond the scope of the present work.
For an even higher viscosity of the oil phase (η s = 9.60 × 10 −2 Pa.s, see figure 3c and supplementary movie SM1), the retraction of the ruptured water film is further slowed down (as evident from the timestamps in figure 3c). Furthermore, the retracting edge also does not possess a thick rim. This observation is similar to the case of Brenner & Gueyffier (1999) for the retraction of viscous films in air, where films of higher viscosity do not form a rim. Moreover, although the expanding holes for the lower η s cases (as shown in figures 3a and 3b) are almost circular, the one for the high viscosity case shown in figure 3c is highly asymmetric. This asymmetry can be attributed to the location of the rupture not being at the film's apex. Since the rupture is happening at an off-apex location, the film thickness at the location of rupture is not spatially uniform due to the curvature of the oil droplet. Hence, the retraction velocity is faster on the part of the film which has a lower thickness. Presumably, this effect is more pronounced when the overall film retraction dynamics are slower, as is the case for the experiments shown in figure 3c. To confirm this hypothesis, one requires high-resolution measurements of the spatial variation of the film thickness, which is challenging in the present experiments (further discussed in § 7). The corrugations at the oil-air-water contact line are also observed in this case. However, since the retraction velocity itself is considerably smaller than for the case shown in figure 3b (see figure 4b for specific values), the tips are not as sharp, and no droplet streams are observed.
To quantify the retraction dynamics, we measure the hole opening radius from each snapshot captured using the high-speed camera. For each experimental snapshot, the area of the hole A(t) is measured, and subsequently an equivalent hole opening radius R(t) is calculated as A(t) = π(R(t)) 2 . A typical measurement from the optical images is depicted in the inset of figure 4a. The temporal variation of the measured hole radius, R, is shown in figure 4a. The time instant corresponding to the first frame in which rupture is optically discernible is denoted by t 0 . Each datapoint in figure 4a denotes the mean of measurements from five independent experiments, and the error bars correspond to ± one standard deviation. In the present work, we focus on the early moments following rupture, as indicated by the red rectangle in figure 4a. Zooming into this early time regime, as shown in figure 4b, it is observed that R(t) varies linearly with time (as  The size of the domain is much larger than the hole radius (Lmax where the oil is more viscous than water (η f = 8.9 × 10 −4 Pa.s), the retraction velocity varies as as evident from the line in figure 4c. This is a weaker dependence as compared to the expected 1/η s variation observed for retraction in the two-phase configuration (Martin et al. 1994;Eri & Okumura 2010). We will attempt to explain the scalings for the twophase and three-phase configurations in § 7. Furthermore, the reason for not fitting the datapoint for the case where the oil is less viscous than the water film in figure 4c will also be addressed therein.

Governing equations
In this section, we discuss the governing equations that describe the retraction of a ruptured liquid film in the three configurations we study in this paper, namely, the classical, two-phase, and three-phase Taylor-Culick retractions. We perform axisymmetric direct numerical simulations using the free-software volume of fluid (VoF) program, Basilisk C (Popinet & collaborators 2013(Popinet & collaborators -2022Sanjay 2021b), which uses the one-fluid approximation (Tryggvason et al. 2011) to solve the continuity and the Navier-Stokes equations: where, v and p are the velocity vector and pressure fields, respectively, η is the viscosity of the fluid, and t denotes time. Furthermore, D is the symmetric part of the velocity gradient tensor D = ∇v + (∇v) T /2 , and f γ is the singular surface tension force needed in the one-fluid approximation to comply with the dynamic boundary condition at the interfaces (Brackbill et al. 1992).

Non-dimensionalization of the governing equations
We non-dimensionalize the governing equations by using the inertio-capillary velocity scale v γ , the thickness of the film h 0 , and the capillary pressure p γ . These scales also define the characteristic inertio-capillary time as τ γ : Here, γ sf is the surface tension coefficient between the film (f ) and the surrounding (s) medium, ρ f the film density, and h 0 its thickness. The dimensionless form of the Navier-Stokes equation (4.2) is where the expressions for the Ohnesorge number (Oh, ratio of visco-capillary to inertiocapillary time scales), the dimensionless density (ρ), and the singular surface tension force (f ) depend on the specific configurations that we discuss below.

Two-phase Taylor-Culick configuration
In this configuration, a liquid film (f ) retracts in a viscous surrounding (s) medium (figure 5a). We use the volume of fluid (VoF) tracer Ψ to differentiate between the film (Ψ = 1) and the surroundings (Ψ = 0), which follows the VoF scalar advection equation, Furthermore, the singular surface tension force is given by (Brackbill et al. 1992): where the curvature κ is calculated using the height-function approach (Popinet 2009). We follow the same sign convention as Tryggvason et al. (2011, see page 33): the curvature is positive if the interface folds towards it normaln, i.e., κ = −∇ · n. Note that the surface tension scheme in Basilisk C is explicit in time. So, we restrict the maximum time step as the characteristic inertio-capillary time based on the wavelength of the smallest capillary wave. Additionally, the density of the film is the same as that of the surroundings, givingρ = 1. Lastly, the Ohnesorge number (Oh) is given by where are the Ohnesorge numbers based on the film and surroundings viscosities, respectively. For this configuration, we keep Oh f constant at 0.05 (based on the experiments of Reyssat & Quéré 2006), and vary the control parameter Oh s in § 5. Note that the computational domain in figure 5a along with (4.5) -(4.8) can be used to simulate classical Taylor-Culick retractions as well by replacing the surroundings (s) with air (a). We discuss the details of the classical configuration in appendix A.

Three-phase Taylor-Culick configuration
In this configuration, we model the bursting of a water film at an oil drop-air interface by simulating the retraction of a fluid film (f ) on an initially flat oil bath (s), while ignoring the effects of the oil drop's curvature (as the retraction length in the early time regime of figure 4b is much smaller than the oil drop radius, see figure 5b). We extend the traditional volume of fluid (VoF) method described in § 4.2.1 to tackle three fluids by using two VoF tracers: Ψ 1 , which is tagged as 1 for the liquids (water film, f , and oil surroundings, s) and 0 for air (a), and Ψ 2 which is 1 for the water film (f ) and 0 everywhere else (figure 5b). Note that this implementation requires an implicit declaration of the surrounding phase (s), given by Ψ 2 (1 − Ψ 1 ) (Sanjay et al. 2019;Sanjay 2021a,b;Mou et al. 2021). Additionally, both Ψ 1 and Ψ 2 follow the VoF tracer advection equation, (4.9) and the dimensionless density ratio is (with ρ f = ρ s ) The Ohnesorge number (Oh) is now given by where Oh f and Oh s follow (4.8), and Oh a = η a / ρ f (2γ sf ) h 0 is the Ohnesorge number based on the viscosity of air. Both Oh f and Oh a are fixed at 10 −1 and 10 −3 , respectively, for all the three-phase simulation data presented in this paper (see § 7), and we vary the control parameter Oh s in § 5. Lastly, the surface tension body force takes the form with γ sa and γ sf being the surface tension coefficients for the surroundings-air and surroundings-film interfaces, respectively.
Physically, such a configuration (figure 5b) and (4.9) -(4.11) ideally imply the presence of a zero thickness precursor film of the surrounding liquid (s, represented by (1−Ψ 2 )Ψ 1 = 1, (4.11)) over the liquid film (f , Ψ 1 Ψ 2 = 1, (4.11)). Note that this numerical assumption is applicable only when it is thermodynamically favorable for one of the fluids (here s) to spread over all the other fluids, i.e., it has a positive spreading coefficient (de Gennes et al. 2004;Berthier & Brakke 2012), S ≡ γ af − γ sf − γ sa > 0, and the Neumann triangle collapses at the three-phase contact line. In reality, this precursor film will have a finite thickness controlled by microscopic forces (like van der Waals forces, Vaynblat et al. 2001), and is much smaller than the length scales that we can resolve numerically in the continuum framework. Indeed, for our numerical simulations, this precursor film has an effective thickness of ∆/2, where ∆ is the size of the finest grid employed in this work. We further assume that, on the time scale of film retraction, the effective spreading coefficient of the surrounding liquid (s) is 0 (Bonn et al. 2009). Consequently, the effective surface tension coefficient between the film and air is γ af = γ sf + γ sa . This precursor film (Thoraval & Thoroddsen 2013) is analogous to the mathematical model for spreading of a perfectly wetting liquid on a solid substrate (de Gennes et al. 2004;Bonn et al. 2009), which regularizes the contact line singularity owing to the numerical slip (with an effective slip length of ∆/2, Afkhami et al. 2018) due to the discretization of the interface.

Note on non-dimensionalization in the viscous regime
For highly viscous surroundings (Oh s > 1), it is convenient to scale the velocities with the visco-capillary velocity scale v η , owing to the dominant interplay between viscous and capillary stresses (Stone & Leal 1989). Further, we can use the visco-capillary time τ η , film thickness h 0 , and capillary pressure p γ to normalize the time, length, and pressure dimensions, respectively: where γ sf is the surface tension coefficient between the film (f ) and the surroundings (s), h 0 the film thickness, and η s the viscosity of the surrounding medium. These viscocapillary scales modify the momentum equation as ρ Oh 2 s ∂ṽ ∂t +∇· (ṽṽ) = −∇p +∇· 2ηD +f γ . (4.14) Here, Oh s is the surroundings Ohnesorge number (4.8),ρ followsρ = 1 and (4.10) for the two-phase and the three-phase configurations, respectively, andf γ equals the corresponding expressions for the two configurations. Additionally, the dimensionless viscosities are given bỹ (4.15) 4.4. Domain size and boundary conditions Figure 5 depicts the computational domains. The left boundary represents the axis of symmetry with origin marked at (0, 0). We set no-penetration and free-slip bound-ary conditions to all other domain boundaries along with zero gradient conditions for pressure. These boundaries are far away from the expanding hole and do not affect its growth. Furthermore, the size of the domain is chosen such that L max max (Oh f , Oh s ), with a minimum L max of 200 for Oh s 1. We have varied this domain size to ensure that the simulations are independent of its value. Note that, if this condition is not met, the assumption of infinite film, which is essential for the theoretical scaling relations developed in this work, will fail (Deka & Pierson 2020).
We employ Adaptive Mesh Refinement (AMR) to correctly resolve the different interfaces as well as regions of high velocity gradients (and hence, high viscous dissipation, see appendix B). To ensure that the velocity field is captured accurately, these refinement criteria (see Sanjay 2021b) effectively maintain a minimum of 40 cells across the thickness of the film (i.e., h 0 /∆ 40). As the apparent three-phase contact line and the viscous boundary layer are critical in the present work, the refinement criteria maintain a minimum of 40 cells in the wedge region near the apparent three-phase contact line. Furthermore, the viscous boundary layer is almost 10 times larger than the film thickness (see § 8.2, figure 11). Consequently, a minimum of 400 cells in the viscous boundary layer in the surrounding medium is needed to properly resolve the velocity gradients. We have conducted extensive grid independence studies so that the final results (energy transfers and the retraction velocity) are independent of the number of grid cells. timescale τ γ ), and the insets contain the growth rate of this hole:Ṙ γ (t) = τ γ dR(t)/dt. After the initial transients, the hole grows (i.e., the film retracts) linearly in time with a constant velocity (v f ). We can use this retraction velocity to calculate the film Weber number,

Taylor-Culick retractions: numerics
which is represented by the black dotted lines in figures 6a and 7a. We f is an output parameter of the retraction process.Note that, for very low Oh s , as the rim grows with time, the inertial drag on the moving rim due to the surrounding medium overcomes the driving capillary forces, resulting in a decrease of the tip velocity (see insets of figure 6a and Jian et al. 2020b). However, we can still calculate a velocity scale (and hence We f ) associated with the Taylor-Culick like retraction immediately after the initial transients (as marked by the black dotted lines in the insets of figure 6a for the lowest Oh s ). Furthermore, when the surroundings is highly viscous (Oh s 1), we plot the growing hole radiusR(t) as a function of time, which is normalized by the visco-capillary timescale τ η , see § 4.3, and figures 6b and 7b. The insets of these panels contain the growth rate of the hole, calculated asṘ η = τ η dR/dt. Once again, we observe that the growth of the hole (and the film retraction) depend linearly on time with a constant velocity, which can be used to calculate the surroundings capillary number,  . Insets of these panels show the variation of the dimensionless growth rate of the hole radius at different Ohs , and mark the definitions of We f and Cas . Lastly, panel (c) illustrates the morphology of the flow at different Ohs atR = 30. In each snapshot, the left hand side contour shows the velocity magnitude normalized with the (terminal) film velocity v f and the right hand side shows the dimensionless rate of viscous dissipation per unit volume normalized using the inertio-capillary scales, represented on a log 10 scale to differentiate the regions of maximum dissipation. Here, the film Ohnesorge number is Oh f = 0.10 and that of air is Oha = 10 −3 . Also see supplementary movie SM3. marked with the black dotted lines in figures 6b, 7b, and the corresponding insets. Ca s is another output parameter of the retraction process. Note that the velocity of the retracting film (v f ) is the same as the velocity scale in the surrounding medium (v s ), following the kinematic boundary condition at the circumference of the growing hole. Consequently, the two output parameters, We f (5.1) and Ca s (5.2) are related as Ca s = Oh s We f (see § 7).
Lastly, figures 6c and 7c illustrate the flow morphologies for the two-phase and threephase configurations, respectively, when the hole has grown toR = 30. Readers can refer to supplementary movies SM2 and SM3 for temporal dynamics of the two-phase and three-phase configurations, respectively. Similar to the classical Taylor-Culick retraction case (appendix A and supplementary movie SM4), both the film and the surroundings move. However, unlike the classical case, even for low Oh s , the surrounding medium takes away momentum from the film owing to inertia (added mass-like effect), thus reducing the retraction velocity (see insets of figures 6a and 7a). Furthermore, contrary to the classical case where the dissipation is highest at the neck connecting the rim to the rest of the film (see appendix A), the dissipation in the other two configurations is spread out, and also occurs in the surrounding medium.
As the hole grows, the retracting film collects fluid parcels from upstream of the moving front and forms a rim (Culick 1960;de Gennes 1996;Villermaux 2020). Essentially, the moving fluid parcels of the retracting tip collide with the fluid parcels upstream of the tip, which were initially at rest. The collisions are inelastic as a fraction of the available energy is dissipated by internal viscous fluid friction. For the classical and two-phase configurations, this rim entails a top-bottom symmetry, which is lost in the three-phase configuration. This is due to the air medium having significantly less inertia (added masslike effect from the properties of the film) than the oil bath, causing the film to dig into the bath and forming a hook-shaped rim (see figures 6c: i-iii and 7c: i-iii). Furthermore, the surrounding bath (s) engulfs the retracting film in order to feed the precursor film. This is a result of the high capillary pressure (high curvature) in the wedge region near the apparent three-phase contact line (Sanjay et al. 2019;Cuttle et al. 2021), which also aids in the formation of the hook-shaped rim (Peschka et al. 2018). Moreover, for the cases where there is a density contrast between the film and the surroundings, the topbottom symmetry can break down even for the two-phase configuration due to a flapping instability at very low Oh s , as discussed by Lhuissier & Villermaux (2009);Jian et al. (2020b). Furthermore, as Oh s increases, the bulbous rim disappears, leading to slender, more elongated retracting films. In the two-phase case, the retraction film maintains (topbottom) symmetry (see figures 6c: iii-v), and the dissipation is highest in the viscous boundary layer in the surrounding medium (see figure 6c; further elaborated upon in § 8.2). However, the three-phase case features (top-bottom) asymmetric films owing to the accumulation of fluid towards the low-resistance air medium (see figures 7c: iii-v), and the dissipation is highest near the apparent three-phase contact line (see figure 7c; further elaborated upon in § 8.2). The disappearance of bulbous rims matches with the experimental observations (see § 3 and Reyssat & Quéré 2006).
Note that the numerical results presented in this section is complementary to the experiments on a film retracting in a submerged oil bath (two-phase case, Reyssat & Quéré 2006) and a film bursting at an air-liquid interface ( § 3). The numerical simulations give us access to the cross-sectional view to elucidate the shape of the retracting films (figures 6 and 7), which is difficult to resolve experimentally. On the other hand, our axisymmetric (by definition) simulations do not show the azimuthal instabilities resulting in the corrugations at the oil-air-water contact line. Furthermore, as we focus only on the early time dynamics, we also neglect the curvature of the oil drop in the case of the three-phase retractions. Nonetheless, we can still sufficiently compare the dependences of the retraction velocity on the Ohnesorge number Oh s of the surroundings (see § 7), along with the scaling relations that we develop in the next section for both the experimental and numerical datapoints.

Taylor-Culick retractions: a force perspective
The capillary and viscous forces, along with the inertia of the film and the surrounding media, govern the retraction dynamics. For the classical configuration (figure 1a), the viscosity and inertia of the outer medium are negligible. Furthermore, the film viscosity η f plays no role in determining the magnitude of the retraction velocity owing to the internal nature of the associated viscous stresses (Savva & Bush 2009), as long as Oh f is less than the aspect ratio of the film (see Deka & Pierson 2020). Using these features, Taylor (1959b) calculated v f = v TC (1.1), resulting solely from momentum equilibrium while disregarding the fate of the liquid accumulated in the rim (Villermaux 2020). In terms of the dimensionless numbers introduced earlier (see § 5), (1.1) implies that is constant and equal to 1 (see appendix A for details of the retraction dynamics in the classical configuration). In this section, we delve into the different realizations of the dominating forces, and their implications, in the two-phase and three-phase configurations.

Two-phase Taylor-Culick retractions
For the two-phase configuration (figure 1b), if the viscosity of the oil phase is small (i.e., Oh s = η s / ρ f γ sf h 0 1), the Weber number based on the film velocity v f and the driving surface tension coefficient (2γ sf ), We f = ρ f v 2 f h 0 /(2γ sf ) (5.1) has a value smaller than 1 (see inset of figure 6a). Nonetheless, the driving surface tension force F γ (t) ∼ γ sf (2πR(t)) (see figure 1b) still balances the inertial force F ρ (t) ∼ ρ f v 2 f (2πR(t)) h 0 . Note that since the oil (surrounding) and water (film) densities are very similar (ρ f ≈ ρ s ), we can still use ρ f for the density scale despite the added mass-like effect. Consequently, in this regime, the Weber number is still a constant during retraction (We f ∼ O (1), inset of figure 6a).
On the other hand, if the viscosity of the oil phase (η s ) is significantly higher (i.e., Oh s 1), the resistive viscous force F η (t) dominates over the inertial effects, as the surroundings Reynolds number Re s ≡ ρ s v s h 0 /η s ∼ O(10 −2 ). In such a scenario, the retraction dynamics will be governed by the balance between the capillary (F γ (t)) and viscous (F η (t)) forces (Fraaije & Cazabat 1989;Reddy et al. 2020), given by For F η (t) in (6.1), one can consider the retracting rim to be a cylinder translating in a viscous flow (Reyssat & Quéré 2006;Eri & Okumura 2010). Thus, the viscous drag can be described by the Oseen approximation to the Stokes flow (Lamb 1975;Happel & Brenner 1983), which to the leading order is expressed as where the factor 2πR(t) is due to the axisymmetric geometry. On equating (6.2) and (6.3), we get Moreover, v f = v s (where v s is the velocity scale in the surrounding medium, see § 5). As a result, (6.4) implies that the capillary number Ca s (5.2) is constant, i.e., Further, upon dividing both sides of (6.4) by the inertio-capillary velocity scale v γ = 2γ sf /(ρ f h 0 ) and squaring, we obtain The aforementioned equations (6.5) -(6.6) denote the scaling laws for viscous two-phase Taylor-Culick retractions.

Three-phase Taylor-Culick retractions
For the three-phase configuration (figure 1c), in the viscous limit (Oh s 1), the force balance is still given by (6.1), especially for the oils that are significantly more viscous than water. Here, the driving surface tension force can be expressed as (from figures 1c and 5b) F γ (t) = (γ sf + γ sa + γ sf − γ sa ) (2πR(t)) = 2γ sf (2πR(t)) , (6.7) assuming the presence of a precursor film of oil on top of the water film (see § 4.2.2 and de Gennes et al. 2004;Bonn et al. 2009;Thoraval & Thoroddsen 2013). However, writing an expression for F η (t) is not as straightforward as the two-phase configuration. As can be observed from figure 7c, during the retraction of the film, the oil climbs on top of the water, resulting in a strong flow in the wedge-like region close to the oil-air-water contact line. The rate of local viscous dissipation in this region is also very high (right panels of figure 7c). Similar wedge flows have also been observed for moving contact lines on solid substrates (de Gennes 1985;Marchand et al. 2012;Snoeijer & Andreotti 2013). It has been reported that the wedge flow results in a viscosity-dependence of velocity that is weaker than 1/η s (Marchand et al. 2012), but the exact nature of the dependence has hitherto not been quantified. The presence of a deformable liquid substrate on which the wedge flow occurs (the retracting water film in this case) complicates the situation even further -making it extremely difficult to arrive at the experimentally-observed v f (η s ) ∼ η −1/2 s dependence (3.1) from a simple force balance. In § 8.2, we attempt to explain this scaling from an energetics point of view. Nevertheless, from the experiments, we know that the v f (η s ) scaling is given by (3.1), which can be rewritten as Dividing both sides of (6.8) by v γ from (4.3) and squaring, we obtain Therefore, from (6.4), (6.6), (6.8), and (6.9), we hypothesize that the presence of the oil-air-water apparent contact line in the three-phase configuration dramatically alters the scaling relationships as compared to the two-phase configuration for Oh s > 1 (see figures 6b and 7b). This will be further elaborated upon in § 8.2. Contrary to this scenario, for low Oh s numbers, the retraction velocities in both these configurations have the same scaling behavior. Despite the presence of a hook-shaped rim in the three-phase case (figure 7c: i-ii), we can still treat the moving rim and the surroundings as lumped elements. As a result, the driving surface tension force γ sf (2πR(t)) still balances the inertial force that scales with ρ f v 2 f (2πR(t)) h 0 , thus giving We f ∼ O (1) (see figure 7a). In the next section, we demonstrate the validity of the scaling relations developed in this section. Figure 8 illustrates the dependence of We f and Ca s on the Ohnesorge number Oh s of the surroundings for the retraction configurations described in figure 1. Note that the same datapoints are presented in both panels 8a and 8b, following the relation Ca s = Oh s We f (as v f = v s , see § 5). In figure 8a, We f = 1 marks the classical Taylor-Culick retraction limit, whereas for the two-phase and three-phase configurations, we identify two regimes: inertial (Oh s < 1) and viscous (Oh s > 1).

Demonstration of the scaling relationships
The inertial scaling is identical for both the two-phase and three-phase configurations: We f ∼ O (1), which also implies Ca s ∼ Oh s (see § 6.1 and 6.2). The brown lines in figure 8 represent these two scaling relations.
The datapoints corresponding to the two-phase numerical simulations (from figure 6) are shown by the dark blue triangles. As Oh s increases, the retraction transitions from the inertial scaling (brown lines), to the viscous two-phase Taylor-Culick scaling: Ca s ∼ Oh 0 s (6.5) or We f ∼ Oh −2 s (6.6). We also plot the experimental datapoints from Reyssat & Quéré (2006) for their silicone oil (surroundings, s) -soap water (film, f ) dataset, shown in figure 8 by the light blue pentagrams. In order to make these datapoints dimensionless, we use h 0 = 100 µm and γ sf = 7 mN/m, denoting the thickness of the soap film and the surroundings-film interfacial tension coefficient, respectively, in their experiments. We also neglect any Marangoni flow, or dynamic surface tension effects. Our simulations and scaling relationships are in reasonable agreement with the experimental datapoints of Reyssat & Quéré (2006). Note that Reyssat & Quéré (2006) tried to fit a trend line of (ln η s ) /η s (higher order Oseen correction) through all of their experimental datapoints to obtain a good fit. When the same datapoints are plotted in figure 8a and b, it is observed that some of their datapoints (corresponding to the low Oh s numbers) are, in fact, in the transition between the inertial and the viscous regimes, while the rest of the datapoints show reasonable agreement with the viscous two-phase retraction dynamics given by (6.5) or (6.6).
We also plot the datapoints corresponding to our experiments (figure 4) and simulations (figure 7) for the three-phase configuration (figure 1c) in figure 8. We observe that, at low Oh s , these datapoints follow the inertial dynamics (brown lines), while at higher Oh s , the datapoints follow the scaling relationships given by (6.8) and (6.9): Ca s ∼ Oh 1/2 s and We f ∼ Oh −1 s , respectively (represented by the black lines in figure 8). Note that in order to non-dimensionalize the experimental datapoints shown in figure 4c (so that they can be plotted in figure 8), one needs to know the film thickness h 0 . In the present experiments, the optical resolution was insufficient for accurate measurement of the film thickness prior to rupture. Moreover, as mentioned earlier, the breakup process itself is highly sensitive to experimental noise (see § 4.2 of Villermaux 2020). Similar difficulties were also presumably experienced by Eri & Okumura (2010) in their experiments of two-phase retraction, and they used a fitting parameter in their v f (η s ) relation, which was a function of h 0 . We also know from bubble bursting experiments (Doubliez 1991;Lhuissier & Villermaux 2012) that the film thickness prior to breakup varies in the range O (100 nm) -O (10 µm). Moreover, in similar studies (Lhuissier & Villermaux 2012;Thoroddsen et al. 2012), the film thickness is retroactively calculated from the retraction velocity measurements. We can see from figure 8 that for low Oh s , the dynamics are independent of the specific nature of the configuration (classical, two-phase, or threephase). Knowing v f , η s , and γ sf , we can calculate the Ca s for the datapoint in figure 4c corresponding to η s = 4.94 × 10 −4 Pa.s. Fitting that Ca s value to the Ca s ∼ Oh s trend line (brown line) in figure 8b, a value of h 0 = 1.5 µm can be calculated, which is within the range observed for previous experiments in a similar system (Lhuissier & Villermaux 2012). Using h 0 = 1.5 µm for the remaining experimental datapoints in figure 4c (for η s > 4 × 10 −3 Pa.s) and calculating Ca s , Oh s , and We f , we find that those datapoints (red circles) also collapse on the trend lines (black lines) along with the numerical simulations (yellow triangles) in figure 8. Note that a water film thickness of h 0 = 1.5 µm sets the Oh f at 0.1 for the three-phase case, which is different from the Oh f that we use for the two-phase case (Oh f = 0.05 based on their experiments of Reyssat & Quéré 2006). Therefore, to justify comparison between the two cases, we varied Oh f in simulations from 0.01 to 0.1 and found that the dimensionless retraction velocities (We f and Ca s ) are Oh f -independent for both the two-phase and three-phase configurations (for Oh f < 1). We also verify the Oh f -independence experimentally by replacing the water in our bath by glycerol-water mixtures, and the measurements thus obtained (green circles) also follow the We f ∼ Oh −1 s and Ca s ∼ Oh 1/2 s trendlines (black lines) in figures 8a and 8b, respectively.
In summary, in § 6 - § 7, we discussed the forces involved during the retraction of liquid films owing to the unbalanced capillary traction, followed by identification of the inertial (Oh s < 1) and viscous (Oh s > 1) regimes in the We f vs. Oh s and Ca s vs. Oh s dependences. We also checked the validity of the corresponding scaling behaviors in this section. To further understand the retraction dynamics due to the capillary traction, we focus on the different thermodynamically consistent energy transfer modes in the next section. Particularly, we try to understand the scaling relationship for the viscous threephase Taylor-Culick retraction that still eludes understanding from a momentum balance point of view (see § 8.2).

Taylor-Culick retractions: an energetics perspective
A retracting liquid sheet loses surface area and consequently releases energy (Dupré 1867(Dupré , 1869Rayleigh 1891;Culick 1960), which further increases the kinetic energy of the system (i.e., film and surroundings). A part of this energy is lost in the process due to viscous dissipation. So, the overall energy budget entails where E γ is the surface energy, E k the kinetic energy, and E d the viscous dissipation. The superscripts account for the film (f ), the surroundings (s), and air (a). Of course, for the two-phase case, the terms associated with air (a) do not exist as there is no air phase. Figure 9 depicts (8.1) for both the two-phase and three-phase configurations. Coincidentally, even for the three-phase configuration, the energies associated with the air phase are negligible (see figure 9, the dashed and dot-dashed lines overlap), even though the velocity field in air is not negligible (figure 7c). We keep E a k (R) and E a d (R) in the energy budget for the sake of completeness. In general, the sum of all these energies at any hole radius R(t) equals the total surface energy at R = 0, i.e., the total energy available to the system. As the film retracts, it continuously releases energy, as its surface energy decreases. Therefore, to calculate (8.1), one can choose a reference for surface energy arbitrarily. In figure 9, the surface energy at a hole radius of R = R max is used as this arbitrary instance. This datum is chosen such that by the time the hole expands to R max , the film would have reached a constant velocity. Furthermore, we can normalize the energies in (8.1) with the total surface energy released as the film retracts to a hole of radius R max . The energy budget now reads Here,Ē(R) = E(R)/(E γ (0) − E γ (R max )), ∆E γ (R) = E γ (R) − E γ (R max ), and R =R(t) = R(t)/h 0 is the dimensionless hole radius. The reader is referred to appendix B for details of the energy budget calculations. Removing the arbitrary datum described above and noting that there is a continuous injection of surface energy (−Ė γ , minus sign because the surface energy is decreasing The energies E are normalized by the total surface energy released as the film retracts creating a hole of radiusRmax = 100 for Ohs 1, andRmax = 1000 (two-phase case) andRmax = 1200 (three-phase case) for Ohs = 100. Note that thisRmax, and hence the surface energy datum, are arbitrarily chosen. We use hole radii that are large enough such that the sheets approach a constant velocity. The superscripts account for the film (f ), the surroundings (s), and air (a).
with the growing hole) into the system, we can also write the energy budgets in terms of rates:  Figure 10. Variation of the rate of change of kinetic energy (Ė k ) and viscous dissipation (Ė d ) as proportions of the rate of energy injection (−Ėγ) with Ohs at steady state for the (a) two-phase and (b) three-phase configurations. For both cases, in the inertial limit (Ohs 1), the fraction of energy that goes into kinetic energy and viscous dissipation are comparable. However, in the viscous limit (Ohs 1), viscous dissipation in the surroundings dominates. Insets show the representative temporal variations of the ratio of the rate of change of energy (Ė) to the rate of energy injection (−Ėγ) with dimensionless hole radiusR at three different Ohs . The superscripts account for the film (f ), the surroundings (s), and air (a). .
We devote the rest of this paper to understanding the distribution of the energy injection rate into the rates of increase of kinetic energy and viscous dissipation for both the inertial and viscous regimes.

Energy transfers in the inertial regime
We first focus on the energy balance in the classical Taylor-Culick retraction and the famous Dupré-Rayleigh paradox (Villermaux 2020). Dupré (1867Dupré ( , 1869 hypothesized that the total surface energy released during retraction manifests as the kinetic energy of the film (Rayleigh 1891). As a result, the predicted retraction velocity was off by a factor of √ 2 (see appendix A), leading to discrepancies with experiments (Ranz 1959;Culick 1960). Nonetheless, it is noteworthy that Dupré (1867Dupré ( , 1869 reached the correct scaling relationship by identifying the essential governing parameters of classical sheet retractions. Culick (1960) identified that the rate of surface energy released (8.9) should be distributed into an increase in kinetic energy of the rim and the viscous dissipation inside the film: The viscous dissipation can be attributed to the inelastic acceleration of the undisturbed film up to the velocity of the edge of the rim. Note that the dissipation is independent of the fluid viscosity and is given by (Culick 1960 where m(t) is the mass of the retracting film. Coincidentally, this rate of viscous dissipation in the film is the same as the rate of increase in its kinetic energy (constant rim velocity, Culick 1960; Villermaux 2020). We confirm this hypothesis in appendix A (see figures 12c, d, and Sünderhauf et al. 2002) Next, we delve into the energy transfers in the two-phase and three-phase configurations. In the inertial limit, in a manner akin to the classical case, the fraction of the rate of energy injection that goes into increasing the kinetic energy is similar to that of viscous dissipation. However, unlike the classical case, the kinetic energy as well as viscous dissipation are distributed among the film and the surrounding medium (figures 9 and 10, Oh s 1). We observe that In a manner reminiscent of Dupré (1867Dupré ( , 1869, we can write where v f = v s (kinematic boundary condition at the tip of the film) and ρ s = ρ f . Additionally, following Bohr & Scheichl (2021) and appendix B, the rate of change of surface energy is given by Using (8.8) -(8.9), and rearranging, we get which is the same as the inertial scaling derived using the force balance (insets of figures 6a and 7a).

Demystifying dissipation in the viscous regime
In the viscous limit (Oh s 1), for both the two-phase and three-phase configurations, the surface energy released is entirely dissipated in the surrounding medium (figures 9 and 10), i.e., −Ė γ (t) ∼Ė s d (t). (8.11) In fact, this interplay between the surface energy and the viscous dissipation sets the velocity scale (v s ) in the surrounding medium, which is equal to the retraction velocity (v f , kinematic boundary condition at the hole). Therefore, to estimate this velocity, we first calculate the rate of viscous dissipationĖ s d (t), which depends on both the viscosity η s of the surrounding medium and the velocity gradients D, following the relation (see appendix B)Ė Here, ε s is the rate of viscous dissipation per unit volume, and the integrals are evaluated over the volume Ω s of the surrounding medium. Note that ε s is highest at the expanding hole, i.e., the tip of the retracting film in the case of two-phase retractions (figures 6c), and the macroscopic contact line in the case of three-phase retractions (figures 7c). The latter is analogous to wetting and dewetting of rigid surfaces ( Figures 11a-i and 11b-i show that the local viscous dissipation increases as we move away from the hole (increasingr). Furthermore, the energy dissipated increases in time as the region of flow expands, owing to the increasing hole radius and the dominant radial flow. To rationalize this increase, we plot the rate of local viscous dissipation per unit circumference of the hole in figures 11a-ii and 11b-ii.
For the two-phase case, the viscous dissipation occurs in the viscous boundary layer δ ν ∼ Oh s √t and saturates atr ≈δ ν (figure 11a-ii). However, for the three-phase case, we can identify two distinct regions of viscous dissipation, the wedge region close to the macroscopic contact line, where the viscous dissipation per unit circumference of the expanding hole increases steeply r < 0.01δ ν , and the viscous boundary layer r < 0.1δ ν , beyond which it saturates (figure 11b-ii). Furthermore, this saturation value gives the total viscous dissipation per unit circumference of the hole, V. Sanjay, U. Sen, P. Kant, and D. Lohsė ( 8.15) which is shown in figure 11c as a function of Oh s . We observe that for the two-phase case, the total dissipation is independent of Oh s , whereas in the three-phase case, it scales with Oh 1/2 Oh 1/2 s 2πR(t) three-phase case. (8.16) Moreover, upon non-dimensionalizing (8.9) using the same scales as used in (8.14), and noting that v f = v s and Ca s = η s v s /(2γ sf ), we get In summary, in this section, we confirmed our hypothesis that the presence of the oil-air-water contact line in the three-phase configuration dramatically alters the scaling relationships and dynamics as compared to the two-phase configuration (see § 6.2). We also relate the dimensionless retraction velocity Ca s with the control parameter Oh s in the viscous limit by following the location and magnitude of the local rate of viscous dissipation during Taylor-Culick retractions in viscous surroundings.

Conclusion and outlook
In this paper, we have studied the effects of the surrounding media on the retraction dynamics of liquid sheets in three canonical configurations. In the classical Taylor-Culick configuration, the interplay between capillarity and inertia of the film results in a constant retraction velocity. We can further neglect the surrounding medium as it does not influence the retraction process. However, for a film retracting in a dense and viscous oil (two-phase configuration), and that at an oil-air interface (three-phase), both inertia and viscosity of the oil phase influence the retraction process. The former presents itself as an added mass-like effect. Even though capillarity still governs the constant retraction velocity, the surrounding medium's inertia reduces the magnitude of the film's momentum as it retracts.
Moreover, when the viscosity of the oil is significantly higher than that of the film, the viscous stresses dictate the retraction process and set the velocity scale. To further demystify the energy balance in this process, we used thermodynamically consistent energy transfer mechanisms to understand the fate of the released surface energy owing to the loss of surface area of the retracting film. This energy is injected into the system and manifests itself as kinetic energy and viscous dissipation. In the inertial regime, the proportions of kinetic energy and viscous dissipation are the same, conforming to the analyses of Culick (1960). However, in the viscous regime, the total surface energy released goes into viscous dissipation in the surroundings.
Following the lumped elements analysis, motivated by Taylor (1959b); Culick (1960), we also developed scaling relations to relate the non-dimensionalized retraction velocity (We f and Ca s ) with the control parameter Oh s . In the inertial limit, the Weber number We f based on the retraction velocity is a constant for all three configurations. On the other hand, in the viscous limit, the retraction velocity in the two-phase configuration scales with the visco-capillary velocity scale (constant capillary number, Ca s ∼ O (1)); while for the three-phase configuration, the capillary number Ca s increases with increasing Oh s , owing to the localization of viscous dissipation near the three-phase contact line.
A natural extension of the present work would be to understand the retraction of non-Newtonian sheets and filaments (Sen et al. 2021) in similar surroundings. In such scenarios, the retraction dynamics will depend not only on capillarity and viscosity as described in this work, but also on the rheological properties of both the film and the surroundings. Furthermore, in a broader perspective, the precursor film-based three-fluid volume of fluid method can be used to elucidate several spreading phenomena, both at small and large scales, e.g., drop-film interactions in the inkjet printing process (Lohse 2022) and late time spreading during oil spillage (Hoult 1972), respectively. The left hand side contour shows the velocity magnitude normalized with the inertio-capillary velocity scale ( v /vγ), while the right hand side shows the dimensionless rate of viscous dissipation per unit volume normalized using the inertio-capillary scales 2Oh (D : D) τ 2 γ , represented on a log 10 scale to differentiate the regions of maximum dissipation. (b) Temporal evolution ofR(t). Time is normalized using the inertio-capillary time scale, τγ = ρ f h 3 0 /γ sf . Inset of panel (b) shows the variation of dimensionless growth rate of the hole radius. Notice that We f = lim R→∞Ṙ γ = 1. (c) Energy budget where the energies (E) are normalized using the total surface energy released as the film retracts, creating a hole of radiusRmax = 150. (d) Variations of the rate of change of energyĖ(t) as a fraction of the rate of energy injection into the system (−Ėγ(t)) with dimensionless hole radiusR(t). The superscripts account for the film (f ) and air (a). The Ohnesorge number of the film for this simulation is Oh f = 0.05, and that of air is Oha = 10 −5 to respect the assumption that the surrounding medium has negligible effect on the retraction process (Taylor 1959b;Culick 1960). Additionally, the air-to-film density ratio is ρa/ρ f = 10 −3 . Also see supplementary movie SM4. area dominates (Ȧ sf Ȧ sa ). Therefore, even for the three-phase configuration, in the steady state,Ė γ ≈ γ sfȦsf .
(B 9) As a result, we only need to evaluateȦ sf for developing a scaling for the rate of change of surface energy. For doing this, we use the analysis presented in Bohr & Scheichl (2021), written in our notations aṡ for a control volume bounded by the control surface A sf (free surface of the film without the rim, figure 13b). Here, U is the velocity of differential control volume bounded by dA sf , κ the curvature at this location, and n is a unit vector normal to dA sf . Lastly, the control surface A sf is bounded by the contour C, and m is a unit vector perpendicular to this contour. Note that capillary traction acts perpendicular to C away from the axis of symmetry. The first term on the right hand side of (B 10) accounts for the change in surface area due to inflation normal to A sf , and the second term is a consequence of the distortion of A sf in the tangential direction, i.e., stretching, or in this case, compression (the growing hole). With this choice of the control surface, the dilation normal to A sf is zero (area inflation only occurs at the rim which we ignore), anḋ where the factor 2 comes in because of the two surfaces (top and bottom). Therefore, for both the two-phase as well as three-phase Taylor-Culick retractions, the rate of injection of energy in the system is −Ė γ (t) ≈ 2γ sf v f (2πR(t)) .

(B 12)
Note that while calculating the rate of change of surface energy, we did not account for the growth of the rim because it is much slower than the growth of the hole, and the flow is predominantly in the radial direction (see figures 6 and 7, and Gordillo et al. 2011).