Discrete self-similarity in interfacial hydrodynamics and the formation of iterated structures

The formation of iterated structures, such as satellite and sub-satellite drops, filaments and bubbles, is a common feature in interfacial hydrodynamics. Here we undertake a computational and theoretical study of their origin in the case of thin films of viscous fluids that are destabilized by long-range molecular or other forces. We demonstrate that iterated structures appear as a consequence of discrete self-similarity, where certain patterns repeat themselves, subject to rescaling, periodically in a logarithmic time scale. The result is an infinite sequence of ridges and filaments with similarity properties. The character of these discretely self-similar solutions as the result of a Hopf bifurcation from ordinarily self-similar solutions is also described.

Free-surface flows can produce a great diversity of patterns such as filaments, drops, bubbles, pearls, etc. [1]. Amongst them, probably the most intriguing and elusive to analyze have been the so called "iterated patterns", i.e. "patterns within patterns" where the same structure repeats itself at different time and length scales. Such structures appear in a wide variety of physical, biological and technological settings: from natural phenomena with fractal features to elasticity and composite materials [2]. In the context of interfacial hydrodynamics, in particular, examples of iterated structures are the formation of several generations of satellite drops in capillary breakup [3], the cascade of structures produced in viscous jets [4] and the iterated stretching of viscoelastic filaments [5]. In this letter, we present for the first time a scenario revealing how such structures may appear via a bifurcation from self-similar solutions to discretely selfsimilar ones where scale invariance occurs only at discrete times, resulting in the infinite repetition of some pattern at a discrete sequence of time and length scales. Discrete self-similarity has proven to be present at some instances of gravitational collapse [6] and has also been proposed as a mechanism for the development of turbulence and formation of singularities in Euler's equation through chaotic self-similarity [7]. Our study reveals the mechanism for discrete self-similarity and ensuing complexity on all scales via a model system consisting of a reduced-order hydrodynamic evolution equation.
The physical situation we consider is the rupture of thin films driven by a destabilizing effect. Liquid films are ubiquitous in a wide spectrum of natural phenomena and technological applications [8]. One well-studied effect is that of long-range intermolecular or van der Waals forces; when the film is sufficiently thin, these forces may cause the film to destabilize and eventually rupture and dewet the substrate. In the long-wave approximation (appropriate for slow flows with strong surface tension), the problem may be formulated in terms of an evolution equation for the film profile h(x, t) in the form h t = −∇ · q, where q = −(h 3 /3µ)∇p is the flow rate, with µ being the liquids's viscosity and p = −σ∇ 2 h − Π(h) being the pressure. The first p component is the Laplace pressure (surface tension times linearized curvature) and the second is the disjoining pressure taken to be Π(h) = −A/h n . In one spatial dimension, and after a suitable nondimensionalization, the evolution equation for h reads In the context of rupture by van der Waals forces, A (strictly, 6πA) is the Hamaker constant, while n is almost always taken to be 3 [9,10]. At n = 3, a remarkable property of solutions of (1) is the development of self-similar film rupture (h → 0 at a single point) in finite time [9,10]. There is, however, good reason to examine different values of n. Yatshishin et al. [11] show that the disjoining pressure with n = 3 is an asymptote to DFT as the distance of the chemical potential from saturation vanishes, assuming a Lennard-Jones potential for pairwise molecular interactions and neglecting screening effects. Thus, the usual form with n = 3 is only approached for thick films. For thin films, there is a deviation from the n = 3 behavior, and, dependent on the system, different exponents might be possible. Furthermore, the behavior on h might be non-local, a consequence of the non-local character of the long-range intermolecular interactions (see also [12,13]). However a widely-adopted algebraic dependence offers ease of access to the corresponding equations facilitating their analytical-numerical scrutiny. At the same time, a great deal of experimental study has shown that the functional form of Π(h) is highly dependent on the nature of the dominant intermolecular force, which is influenced by the substrate and liquid properties. For example, [14] found that the effects of screening can lead to a potential better modeled with n = 4 in hydrocarbon-metal experiments with h > 40 nm, while the contribution due to hydrogen bonding in water-silica-glass experiments for very thin (< 30 nm) films was better approximated by n = 1 (Pashley [15]). For water on quartz, Π(h) is estimated to also have n = 1 for h < 80 nm and n = 2 for h > 120 nm [15]. For 80 nm < h < 120 nm, it would then be appropriate to take 1 < n < 2. In any case, given a liquid and a substrate, we can approximate, when appropriate, Π(h) for relatively thin films with a power law by fitting a value of n. We emphasize that assuming a certain form for the disjoining pressure and fitting appropriate values of its parameters is common in the literature (e.g. [16]). There has also been recent interest in generalizing the standard Lennard-Jones potential with attractive exponent λ a = 6 to other exponents, leading to the so-called Mie potential [17]; the corresponding disjoining pressure has exponent n = λ a − 3 [19]. A good summary of different contributions to the disjoining pressure may be found in [18]. The dynamics of rupture under these different values has not previously been examined.
As well as intermolecular forces, (1) models other thinfilm phenomena at different scales, such as destabilization due to thermocapillarity [20] and density contrast (Rayleigh-Taylor instability) [21]; in such cases we may define an equivalent "disjoining pressure" behaving as ln(h) for the thermocapillary effect (essentially n = 0), or as h for the Rayleigh-Taylor instability (n = −1). Instead of self-similar rupture, these two examples exhibit cascades of satellite droplets, similar to those discussed above, so it is of great interest to understand how the two behaviors are connected through variation in n.
Assuming that rupture occurs at a single point x 0 at time t 0 , it is natural to seek solutions in a coordinate system that focuses on the point and time of rupture: where, from simple dimensional arguments based on (1), one finds α = 1/(2n − 1), β = (n + 1)/(4n − 2). For a rupture solution to exist, we must assume n > 1/2, so that α, β > 0. The scaled profile f (ξ, τ ) then satisfies subject to the condition that the interfacial velocity h t remains finite at a finite distance from x 0 . As t → t 0 , one has ξ → ∞ and, in order to cancel out singular dependence on t 0 − t, we must impose Steady states of (3), (4) represent self-similar solutions of (1). Including τ dependence, allows us to examine the stability and dynamics in the vicinity of these solutions. Above a certain value of n, there are infinitely many steady states of (3), (4). This was established for n = 3 in [9] and recently extended to general n by the authors [22]. These solutions are symmetric and can be arranged (for a given n) as a sequence f 1 , f 2 , . . . according to their values at ξ = 0 such that f 1 (0) > f 2 (0) > · · · . We can thus depict solution branches as f j (0) over n (Fig. 1). These solutions were computed using the open source numerical continuation software AUTO-07p [23] (see also [22,24]). As n is decreased, the solution branches merge, with the first two branches merging n = n c ≈ 1.49915 [22]. Our focus is on the change in dynamics of solutions to (1) close to this value.
First, we analyze the stability of f 1 (ξ) as n varies. We linearize (3) about f 1 (ξ) and seek solutions of the linearized problem in the form e στ Φ(ξ), obtaining an eigenvalue problem for σ = σ R + iσ I . As f 1 is symmetric in ξ, Φ(ξ) may be either symmetric or antisymmetric, which we enforce by applying the appropriate boundary conditions at ξ = 0, in addition to the far-field conditions arising from (4). For n = 3, it has been shown [10] that there are two trivial modes of perturbation, symmetric with σ 1 = 1 and antisymmetric with σ 2 = β, which correspond to time and space translation of the singularity, respectively. Otherwise, all eigenvalues have negative real part, and so f 1 is stable. All other branches f 2 , f 3 , . . . have eigenvalues with positive real part and are unstable. Fig. 2 displays results for general n. As well as the two trivial eigenvalues, we compute the two nontrivial eigenvalues with largest real part corresponding to symmetric and antisymmetric modes. These both have negative real part at n = 3 but increase as n decreases, crossing the imaginary axis at Hopf bifurcations close to n = n c where f 1 and f 2 merge. These points are labelled n s and n a . In general, a Hopf bifurcation leads to the existence of a branch of periodic orbits (in scaled time τ , in this case) emanating from the bifurcation.
We now explore the implications of this loss of linear stability on the nonlinear dynamics by computation of the time-dependent equation both in the unscaled (1) and scaled (3) coordinates. To compute solutions to (1) that can capture details close to rupture, we implement an adaptive finite difference scheme that increases local mesh refinement near the minimum of h whenever h min is less than half of its value at the previous mesh refinement. Fig. 3 shows the results of the computations for (a) n = 1.7 and (b) n = 1.5, which are on either side of the Hopf bifurcation structure shown in Fig. 1. Results for other n values are included in the supplementary material. The transition from classical (continuous) self similarity to the onset of cascading oscillations of geometrically decreasing size, is apparent. The inset in Fig. 3a shows that the profiles approach a classical selfsimilar profile [i.e. a steady state of (3)] for n = 1.7. In Fig. 3b, we observe the repetition of the same pattern on geometrically smaller scales, which asymptotically approaches the scaled-time periodic solution to (3), which we describe next.
The computation of solutions to (3) is complicated by the trivial eigenvalues corresponding to shifts in space and time. These instabilities may be thought of as arising from incorrect choices of x 0 and t 0 in scaling the initial condition. We remove these instabilities by letting x 0 and t 0 be time-dependent estimates of the true rupture In each case, the scaled behaviour asymptotes to the stable steady state of (3) for n = 1.7 and the periodic orbit for n = 1.5 (see Fig. 4).
location, which leads to a new equation of the form where P and Q are extra degrees of freedom that may be fixed by applying nonlocal constraints that ensure the rupture remains, or at least asymptotically approaches, ξ = 0. We determine P and Q by approximately fixinĝ f (0,τ ) = 1 in addition to an integral 'pinning' condition [27]. Solutions of (5) are scaled back to those of (3) by In Fig. 4, we plot the results of computations of (5) for n = 1.5, transformed to solutions of (3) via (6), starting with a generic (asymmetric) initial condition, and run until it is clear that a periodic orbit has developed. This provides numerical confirmation that stable periodic solutions to (3) do exist. We found that development of a periodic orbit is sometimes prevented by rupture occurring away from the origin in (5); this is particularly dependent on initial condition and becomes harder to avoid for either n closer to the Hopf bifurcation structure or small values of n. Of particular note in Fig. 4 is that the periodic orbit is asymmetric, and exhibits oscillations that advect outward from the minimum of f to the far field. These oscillations correspond to the cascade of oscillations, seen in Fig. 3, that are asymptotically fixed in unscaled space as t → t 0 .
A periodic orbit in self-similar coordinates implies that the rupture of the film occurs in a discretely, rather than continuously, self-similar fashion; self-similarity of profiles only holds at discrete times t 1 , t 2 , . . . , approaching the rupture time t 0 geometrically; if T is the period of the orbit, then t N +1 /t N = e −T . Such behavior has been referred to as discrete self-similarity [25] and linked to the existence of periodic orbits in scaled coordinates; the results in this letter comprise the first explicit computation of such a periodic orbit [26]. We may understand the outward propagation of peaks and troughs in the solutions to (3) in the scaled coordinates as the creation of 'drops' and necks between drops of geometrically shrinking scale in solutions to the unscaled problem (1), thus leading to fractal-like profiles at rupture (as seen in Fig. 3b).
The geometric factor in question depends both on α, β and T . Suppose the maxima h 1 , h 2 , . . . are located at distances d 1 , d 2 , . . . from x 0 (with d N → 0 as N → ∞). Successive maxima correspond to the same maximum in (ξ, f ) at scaled times τ and τ + T . Using h = e −ατ f and x − x 0 = e −βτ ξ, we deduce The period observed for n = 1.5 is T ≈ 6.1, while the periods at the symmetric/asymmetric Hopf bifurcations are 2π/0.912 ≈ 6.9 and 2π/0.885 ≈ 7.1, respectively.
Recently the transition from continuous to discrete and then chaotic self-similar dynamics has been observed in the context of slip instabilities in elasticity-in particular, simulations of the frictional sliding dynamics of two elastic bodies in contact [28], in which there is a series of Hopf bifurcations on the branch of (self-similar) steady states. As far as we are aware, Eq. (1) is the first model in hydrodynamics in which a periodic orbit in self-similar variables has been observed. Our future aim is to systematically compute solution branches from the Hopf bifurcations. While speculative at this point, one possibility is that on such branches there may be additional bifurcations to quasi-periodic solutions and the system may become chaotic (via e.g. Ruelle-Takens-Newhouse route) as n decreases. For n = 1 (see Fig. 5), a cascade of satellite drops of successively smaller sizes is observed. However, at the space between satellites, new cascades of subsatellites develop. In between subsatellites at the same cascade, subsubsatellites develop etc. The result resembles a fractal-like structure. The minimum height h min does not follow (or oscillate around) a predicatable power law, since the position where it is reached keeps jumping from one cascade of (sub)satellites to another. We do not study this process in detail here, but present it to demonstrate the complexity that develops as n is decreased.
As supplemental information we include more details on linear stability and numerical methods.
We are grateful to Andrew Archer and Amparo Galindo for valuable comments and suggestions on disjoining pressure and intermolecular interactions, respec-tively. We acknowledge financial support from the Engineering and Physical Sciences Research Council (EP-SRC) of the UK through Grants No. EP/K041134/1, EP/K008595/1 and EP/L020564/1 and from the Spanish government through Grant No. MTM2014-57158-R. MCD was employed by Imperial College London while undertaking the work reported on in this paper.