Self-Similar Cusp Formation in Thin Liquid Films By Runaway Thermocapillary Forces

Many physical systems give rise to dynamical behavior leading to cuspidal shapes which represent a singularity of the governing equation. The cusp tip often exhibits self-similarity as well, indicative of scaling symmetry invariant in time up to a change of scale. Cusp formation can even occur in liquid systems when the driving force for fluid elongation is sufficiently strong to overcome leveling by capillarity. In almost all cases reported in the literature, however, the moving interface is assumed to be \textit{shear-free} and the operable forces orient exclusively in the direction normal to the advancing boundary. Here we focus on a system in which a slender liquid film is exposed to large thermocapillary stresses, a system previously shown to undergo a linear instability resembling microlens arrays. We demonstrate by analytic and numerical means how in the nonlinear regime these surface forces undergo self-similar runaway behavior leading to cusp formation with a conical tip whose slope can be prescribed from the analytic relation derived. On a fundamental level, this finding broadens our understanding of known categories of flows capable of cusp formation. More practically, the system geometry proposed offers a potentially novel lithographic method for one-step non-contact fabrication of cuspidal microarrays.

-Liquid (a) line type and (b) point type cusp formation in a thin film subject to a negative disjoining pressure from van der Waals forces that promotes dewetting of the film from the bottom solid substrate [12]. Liquid (c) line type and (d) point type conical cusp formation caused by thermocapillary forces which draw fluid away from the lower warm substrate toward the colder top substrate, described in more detail in the text.

Cusp formation in physical systems
Despite that capillary forces always act to repress regions of high curvature, nature nonetheless finds clever ways of forming and sustaining cusps in many physical systems. In fact, cusps are rather ubiquitous and occur in such diverse phenomena as thermal grooving at grain boundaries [1], surface diffusion and pinchoff in annealed or sintered systems [2], complex plasma formations [3], wavefront propagation in systems described by the linear [4] or nonlinear Schrödinger equation [5], critically charged droplets [6], microbranching instabilities in fast moving cracks [7], line attractor states in neural computation models [8] and many more. A recent delightful book by J. Eggers [9] describes as well the complex dynamics governing cusp formation in many liquid systems including thread and droplet breakup, Hele-Shaw sink flow, and thin film rupture caused by a negative disjoining pressure which induces a dewetting process [2,10,11]. The latter system is sketched in Fig. 1 (a) and (b) where the receding air/liquid interface is observed to form a cuspidal curve.
In these and other systems [13,14,15,16,17], the apical region of the evolving cusp exhibits self-similar behavior characterized by universal exponents, some of which have been confirmed experimentally [18,19,20,21,22]. The resulting power laws stem from scaling symmetries that are invariant in time up to a change of scale. In almost all cases reported in the literature, however, the moving interface is assumed to be shear-free where the operable forces orient exclusively in the direction normal to the advancing boundary. The interface therefore experiences no shear forces and plays no active role in corralling fluid into a sharpened tip. And while there have been observations of cusp formation leading to tip streaming in droplet systems subject to interfacial shear from surfactant concentration gradients [23,24], the dynamics of cusp formation there  [26]. (b) AFM image of cicada wing [27].
remains an unsolved problem.
To explore cuspidal formation driven by shear forces at a free interface, we here focus on a thin film system designed to elicit self-reinforcing thermocapillary stresses at the air/liquid interface. We analyze the dynamics by which the ensuant self-similar process gives rise to fluid elongations shaped like liquid cusps whose conical tips further promote self-focusing. Shown in Fig. 1 is an example of a thermocapillary (c) line and (d) point cusp caused by runaway thermocapillary forces. While Figs. 1 (a) and (b) depict cusp formation arising from forces exclusively oriented normal to the free interface (disjoining pressure counterbalanced by capillary pressure), Figs. 1 (c) and (d) depict formation of a cusp from thermocapillary (shear) forces which orient parallel to the moving interface. An additional distinction between these two thin film systems is that the thermocapillary problem exhibits multiscale dynamics in the apical region, which although quite interesting, considerably complicates any stability analysis.
Aside from such fundamental considerations, there is practical motivation for this study as well. We are interested in exploring thermocapillary based techniques for patterning thin liquid films which can be rapidly solidified in situ. The system geometry examined in this work offers a potentially novel lithographic method for onestep non-contact fabrication of cuspidal microarrays. This development can facilitate design and manufacture of specialty microarrays such as biomimetic cuspidal substrates. Two recent important examples include infrared (IR) antireflective moth eye surfaces patterned with quintic cusps for eliminating Fresnel reflections in the mid-IR [25,26], as shown in Fig. 2 (a), and superhydrophobic, self-cleaning antimicrobial surfaces mimicking the surface of a cicada wing [27,28], as shown in Fig. 2(b). Such surfaces can likely be architected using thermocapillary forces to create substrates for which form follows function i.e. imprinted cuspidal shapes relate directly to their intended function.
Our group has previously demonstrated experimentally [29] how patterned thermocapillary forces can be used to sculpt nanofilms into liquid microlens arrays, which are then solidified rapidly in situ. The resulting ultrasmooth surfaces are ideally suited to micro-optical applications such as beam shaping. The analysis presented in this work now suggests that were the microlens configuration allowed to evolve further in time, the system would transition to a microcuspidal array. The local analysis presented in this work indicates how initial protrusions of any sort, whether triggered by the linear instability [30,31] or triggered by large amplitude perturbations [32,29], are expected to evolve into individual or array-like cuspidal patterns.
The outline of this work is as follows. In Section 2 we present the thin film evolution equation for a molten Newtonian which are imposed by thermal conduction across a very slender confined system. This system gives rise to an initial linear instability whose wavelength characterizing the fastest growing mode is subsequently used to rescale the original equation to dimensionless form. Further rescaling to parameter-free form yields an equation belonging to the general class of so-called gradient flows, which in this case also contains a virtual singularity where thermocapillary stresses diverge to infinity. In Section 3, it is shown that this evolution equation equation does not support any stable stationary states since for small excursions about a stationary state there exist states of even lower energy. This demonstration proves that the dynamics incurred by the confined geometry imposed on the film leads to a runaway thermocapillary process in which the liquid can reduce its free energy by advancing ever closer to the top colder substrate. In Section 4, 2D and 3D numerical solutions of the nonlinear evolution equation reveal stable formation of a cusp with a conical tip that undergoes continuous sharpening by a self-similar process exhibiting power law growth in the tip speed and tip curvature. In Section 5 we present an asymptotic analysis of the apical region about the virtual singular point which reveals the presence of a stable fundamental mode which appears to act an attractor state. Various measures characterizing this fundamental mode are shown to be in excellent quantitative agreement with numerical simulations of the late time asymptotic behavior of the apical region. The asymptotic analysis also reveals a simple relation for the conical tip slope that can be used to prescribe the shape for experimental applications. In Section 6, we conclude with some final thoughts on how these findings may held advance a novel lithographic method for fabrication of specialty cuspidal microarrays.

Long wavelength thermocapillary model for growth of protrusions by increasing interfacial shear forces
A theoretical model has previously been derived [30,31] to describe the evolution and stability of a gas/liquid interface for the system sketched in Fig. 3. A molten nanofilm of initial uniform or average thickness h o overlay by a slender gas film is confined within a very narrow gap d o (typically less than a micron) by two opposing substrates maintained at a uniform temperature difference ∆T = T hot − T cold > 0. The model assumes that the film thickness is much smaller than any characteristic lateral scale, that inertial forces are negligible, and that the viscosity of the film µ = µ(T hot ) is relatively constant given the very small gap width dimension d o . The film thickness is restricted to the range 0 < h(x, t) < d o where x = (x, y). Since for single component fluids the variation in surface tension γ with temperature T given by dγ/dT is a negative quantity, those portions of the liquid film which are closer to the cold substrate experience a colder temperature and consequently a higher local value of surface tension. Temperature variations along the liquid interface therefore give rise to spontaneous interfacial thermocapillary stresses given by ∇ γ = (dγ/dT )∇ T , which act to pull liquid from warmer to cooler regions of the film. Within the long wavelength approximation, the operator ∇ denoting the surface gradient simply reduces to (∂/∂x, ∂/∂y). In this limit, the corresponding energy equation describing heat transfer across the gas/liquid bilayer reduces simply to the 1D Laplace equation d 2 T /dz 2 = 0 from which can be derived the temperature distribution along the liquid interface z = h(x, t): The material parameter κ denotes the ratio of gas to liquid thermal conductivity evaluated at the temperatures of the respective adjacent substrates. Since the gas layer is always more thermally insulating than the liquid layer, the ratio κ is restricted to the range 0 < κ < 1. Depending on the materials of choice, however, the magnitude of κ can range anywhere from about 1/4 or higher for molten polymer films like polystyrene overlay by an air film [31] to 10 −4 or smaller for liquid metal films such as indium [33] overlay by a xenon gas layer [34]. The confined geometry leads to self-reinforcing thermocapillary stresses, which promote growth of elongations toward the colder substrate. This process is mitigated only by capillary forces which try to repress formation of regions of high interfacial curvature. Stabilizing gravitational forces, which are orders of magnitude smaller than the thermocapillary forces, are completely negligible. The system described has also been shown to be susceptible to a linear instability [30,31] which establishes irrespective of the size of the applied thermal gradient. At early times, infinitesimal disturbances generate periodic undulations in film thickness which undergo exponential growth. The fastest growing undulations are characterized by the wavelength where γ o = γ(T hot ) and γ T = |dγ/dT | T hot . All else equal, a larger difference in temperature ∆ leads to growing undulations of smaller wavelength. Recent [35,36,37] and ongoing experiments to confirm the mechanism leading to instability so far indicate good agreement with analytic predictions for the fastest growing wavelength and its growth rate. The dimensionless evolution equation describing the long wavelength thermocapillary model is given by where u c is chosen to be a characteristic fluid speed based on in-plane thermocapillary flow. Details of the analysis and derivations leading to this form have been previously presented elsewhere [31]. The thin film behavior is therefore controlled by two dimensionless numbers, namely a modified Capillary number Ca = µu c / 3 γ o and a modified Marangoni number M a = γ T ∆T /µu c . These numbers differ from their usual definitions by factors of the small parameter = h o /λ max intrinsic to the long wavelength approximation. (This parameter should not be confused with the small parameter ε pertaining to temporal behavior introduced in Section 5.) As evident, Eq. (3) exhibits a virtual singularity at H s = D/(1−κ) (equivalently in dimensional variables h = d o /(1 − κ)). This singularity lies outside the physical domain and beyond the top cold substrate since κ is always less than one. In general, the system described by Eq. (3) is not limited to initially flat liquid configurations and describes equally well the response of any initial non-uniform liquid state to thermocapillary forces, in which case h o denotes the average initial film thickness.
For purposes of this current study, it proves convenient to recast Eq. (3) into parameter-free form such that where H = H/H s , X = X/X c , ∇ = X c ∇ and τ =τ /τ c . The scalings for this reduction are given by X c = (2 DH s /3κM a Ca) 1/2 and τ c = 4 D 2 /(3κ 2 H s M a 2 Ca). In this final form, the top cold substrate is located at H = 1 − κ while the virtual singularity occurs at H = 1. Since in this work we wish to investigate the long time behavior of film protrusions which evolve into cuspidal shapes, we restrict attention to small values of κ 2 × 10 −4 (characteristic of a liquid metal film overlay by a highly insulating gas layer). This allows for a longer evolution interval not prematurely terminated by contact with the top substrate. Such contact, of course, would modify the fluid behavior in ways not described by Eq. (4).

Stability considerations by analogy to general gradient flows
In previous work [30,31], we presented the linear stability analysis of Eq. (3) which exclusively focused on early time behavior of infinitesimal fluctuations in interfacial temperature or film thickness. That analysis showed that the instability is of Type II [38] where all modal fluctuations of wavelength λ > λ max / √ 2 are linearly unstable irrespective of the value ∆T . Eliciting the stability characteristics of stationary states of the full nonlinear equation given by Eq. (4) requires a different approach based on the system free energy F[H]. By exploiting an analogy to gradient flows in general, we show next that Eq. (4) does not admit any stable stationary states on a periodic or infinite domain so long as H > 0.
Mitlin [39] has previously shown that the interface equation describing thin film dewetting by van der Waals forces, the process depicted in Figs. 1 (a) and (b), can be rewritten in Cahn-Hilliard form described by known more generally as gradient flow form [40]. The equation governing the thermocapillary model can also be written in this form for a free energy functional given by which reflects the surface pressure required for maintaining stationary states of constant volume V . It has been shown that for a general class of thin film equations [41], which includes the form of Eq. (4), small perturbations to periodic stationary states (i.e. δH ∝ ∂ 2 H/∂X 2 ) lead to negative values of the second variation whenever the potential function satisfies the relation d 4 U/dH 4 | H∈H < 0. This negative value indicates that there are always nearby states with same periodicity as H of lower free energy. The proof is provided in Appendix B. Since for the thermocapillary model the curve d 4 U/dH 4 shown in Fig. 4 is always negative, this therefore proves that Eq. (4) cannot support stable stationary periodic states. This analysis is quite general and can be applied to many other thin film systems (even volume non-conserving systems) so long as the governing interface equation can be cast into the form of Eq. (5). We next focus on the nonlinear evolution of liquid shapes and show how the thermocapillary model promotes formation of self-similar cusps. , which also exhibits self-similar collapse. The virtual singularity H apex = 1 appears therefore to act as an attractor state for formation of the conical tip.

Numerical solution of nonlinear thermocapillary model equation
The top panel shown in Fig. 6 represents 3D views of an evolving cusp with a conical tip at for the four times designated, as obtained from finite element simulation of the full nonlinear Eq. (4). Additional information about the simulation and accompanying video clip can be found in Appendix C. The bottom panel displays the value of the curvature of the gas/liquid interface at every point within the computational domain. The orange curves delineate concave from convex regions. The last image in the bottom panel clearly reveals that the interface evolves into a true cusp capped by a conical tip of decreasing radius flanked by a broader convex surface.

Asymptotic analysis of self-similar cusp formation
The exponents extracted from the numerical simulations described in the previous section are also confirmed by analysis of Eq. (4) by considering a Taylor expansion  about the virtual singular point H = 1, which yields the asymptotic evolution equation Balancing the first and second term with the second and third term yields the same asymptotic relation obtained previously, namely X ∼ 1 − H apex ∼ (τ s − τ ) 1/4 . Based on these scalings, we introduce the stretched variables Note that if Eq. (4) were truly scale invariant, and not just asymptotically so as H apex → 1, the expansion in Eq. (10) would terminate at n = 1. The appearance of the 1 − H term in the denominator of Eq. (4), however, precludes such global scaling and instead leads to multiscale expansions of the form: ε n S n (w 1 , . . . , w n ) (12) where the operator symbols ∇ , ∇ · and ∇ 2 are understood to reduce to the appropriate rectilinear (X) or cylindrical (R) form for the gradient, divergence and Laplacian operations. To leading order n = 1, Eq. (9) then reduces to the nonlinear, fourth order equation given by where the operators T 1 , S 1 and M 1 are defined as Here and in what follows, operator subscripts denote differentiation with respect to the self similar variable η. Required symmetry about the axis of origin yields two boundary conditions, namely (dw 1 /dη) η=0 = 0 and (d 3 w 1 /dη 3 ) η=0 = 0. An additional boundary condition is obtained from the requirement that Eq. (11) remain bounded as ε → 0, or equivalently as η → ∞, which requires that the leading term T 1 vanish. This then leads to the Robin condition T 1 (w 1 )| η→∞ = 0. To leading order then, the asymptotic solution to Eq. (14) is satisfied by the Laurent series Convergence to w ∞ 1 can be obtained by linearizing Eq. (14) about the solution w 1 (η → ∞) = w ∞ 1 (η) + f (η) which leads to the non-homogeneous linear equation In the limit |a 1 | 1, this equation leads to a singular perturbation problem whose inner region is influenced by the fourth order capillary term (not shown). Here we only focus on the global outer region solutions of the linearized equation obtained by WKBJ analysis where f (ση) = exp σ −4/3 ∞ n=0 σ 4n/3 S n (ση) for σ 1. Matching terms of order σ −4/3 and σ 0 , then solving for the resulting two ordinary equations yields the general solution where α = 1 for rectilinear and α = 5/3 for axisymmetric geometry. To preclude the first two terms in the summation from undergoing diverging oscillatory behavior, it is required that β 1 = β 2 = 0. The two remaining non-vanishing terms proportional to β 0 and β 3 simply reflect an infinitesimal shift in the far field slope and a rapidly decaying function, respectively. Were the analytic solution to Eq. (14) known within the apical region, then the coefficients β 0 and β 3 could be obtained by asymptotic matching. Absent that information, the solutions to Eq. (14) are still constrained by the symmetry requirement about at the origin. This constraint imposes that the solutions correspond only to discrete values of the far field slope, as discussed next.
The numerical solutions to Eq. (14) were computed on a finite domain sufficiently long to preclude finite size effects. Simulations with increasing mesh refinement were conducted to assure convergent solutions. Shown in Fig. 7 are the first six similarity solutions with selected numerical values listed in Table 1. The asymptotic interface slopes in the conical region for axisymmetric geometry are always smaller than the slopes for rectilinear geometry, as expected. The axisymmetric solutions also display weaker oscillatory behavior, likely due to suppression by the capillary pressure associated with the additional term in the interface curvature. The fundamental mode p = 1 exhibits no oscillatory behavior unlike the higher order solutions p ≥ 2. Table 1 -Asymptotic values of the interface slope, apex height and apex curvature for the leading order solution w 1 to Eq. (14). Numbers in blue and red denote values for rectilinear and axisymmetric geometry, respectively.   Table 1 into these expressions yields intercept values for (∂H/∂τ ) apex equal to −1.632 (rectilinear) and −1.681 (axisymmetric). Likewise, the intercept values for (∇ 2 H) apex equal −0.175 (rectilinear) and −0.078 (axisymmetric). These predicted values are in excellent agreement with the numerical intercept values (shown in parentheses) in Fig. 5 (c). Additionally, the asymptotic values of the interface slope lim η→∞ dw (1) 1 /dη given in Table 1 also show excellent agreement when superposed on the profiles in Fig. 5 (d). The asymptotic values are predicted to be 1.0437 (rectilinear) and 0.7639 (axisymmetric), while the numerical results yield 1.044 and 0.764. Converting back to dimensional form, the value of the interface slope in the region of the conical tip is given by the relation With the value of the asymptotic slope known, the remaining parameters in Eq. (21) are set by the material constants of the gas/liquid system of choice and the temperature drop applied to the confining substrates. In Section 4, it was shown that the numerical solution to the full nonlinear equation given by Eq. (4) asymptotes to a fluid shape resembling a cusp capped by a conical tip. The asymptotic analysis in this Section reveals that the numerical solution obtained corresponds identically to the fundamental solution w is beyond the scope of this paper. Consideration of this issue by implementing a conventional linear stability analysis of Eq. (9) is non-trivial due to the multiscale nature of the asymptotic, self-similar base state solutions, which evolve on multiple time scales {ε n } ∞ n=1 . Since both the numerical and analytic solutions suggest that the late stage dynamics of Eq. (10) is dominated by the term w (p) 1 , it suffices then to consider infinitesimal perturbations described by where |φ (p) m (η)| 1 denotes an infinitesimal modal perturbation to w (p) 1 (η), θ is the polar angle in cylindrical coordinates, and ε is defined in Eq. (10). The resulting eigenvalue problem is given by where 1 ) 2 and where differential operators in Eq. (14) have been expanded to include the appropriate θ-dependence. In order for localized perturbations in the far field to preserve constant slope, invariant, there also exist for each value of p two eigenvalues reflecting these symmetries, namely the eigenfunction cos θ × dw 1 is therefore the only solution that is linearly stable to perturbations. The remaining positive eigenvalues increase in magnitude with increasing p, indicating more rapid growth and instability associated with the coefficient ε 1−4λ multiplying the last term in Eq. (22). The numerical simulations described in Section 4 and plotted in Figs. 5 and 6 were always found to asymptote to the bounded fundamental solution w (1) 1 . Similar strong convergence to the stable fundamental solution has also been reported for the thin film equation describing van der Waals rupture [43] (shown in Fig. 1). In that example, initialization of the thin film equation by the corresponding solution w (p≥2) 1 for that problem leads to a different global liquid film configuration -however, the local behavior in the vicinity of a line or point of rupture nonetheless converges to the fundamental mode w (1) 1 . A full investigation of the local scaling behavior leading to self-similar cuspidal formation in the thermocapillary system for initial conditions resembling higher order eigenmodes is left for further study. It is anticipated that irrespective of the initial condition, simulation of the full nonlinear evolution equation given by Eq. (4) will still yield film shapes dominated by w 1 is linearly stable.

Conclusion
The analysis and simulations presented in this work predict how surface shear forces from self-reinforcing thermocapillary stresses at a gas/liquid interface will produce sharp protrusions resembling cusps capped by conical tip. This finding expands the category of hydrodynamic flows known to form stable cusps to include thin film systems subject to interfacial shear, where the driving force is oriented parallel to the moving interface. The asymptotic analysis reveals how the liquid tip undergoes self-focusing toward a virtual attractor state characterized by a line (rectilinear case) or point (axisymmetric case) singularity via a persistent self-similar process. The asymptotic derivation also yields an simple analytic relation for the slope of the conical tip which should prove useful to experimentalists who may require microarrays with specified tip slope, for beam shaping purposes, design of super antireflective coatings [26], or other applications.
The original system described, based on a thin uniform molten film confined by parallel solid boundaries maintained at different uniform temperature, is known to support a linear instability that forms arrays of rounded protrusions resembling microlenses. These protrusions are expected to evolve into arrays of conical cusps by the nonlinear dynamical process described since the thermal gradient in the region above the fluid tip becomes increasingly large with time, leading to a runaway process. We anticipate that any initial film configuration that contains local maxima in film thickness, whether or not periodically arranged and however initially seeded, will also trigger cusp formation at those locations given the local, self-similar nature of the underlying growth dynamics.
We have previously shown [35,29,36,37] that the evolution process leading to rounded lenslet microarrays can be terminated on demand and the liquid shapes affixed in place by dropping the temperature of both substrates below the solidification point. Rapid solidification of these liquid structures is made possible by two advantageous features: the large surface to volume ratios intrinsic to microscale or nanoscale films which facilitates rapid cooling and digital control over the temperature of the confining substrates. It is fully expected then that similar rapid solidification can be achieved once the desired conical protrusions have formed in order to solidify and affix their shape on demand. Perhaps alternative methods of flow control by laser manipulation, previously applied to thin film thermocapillary spreading along a solid substrate, can also be used [44]. In summary, the findings and implications outlined in this work offer a novel lithographic method for direct, non-contact fabrication of cuspidal microarrays, whose shapes are more difficult, costly or even impossible to fabricate by other means. Figure C1 -Image of progressively refined mesh used to resolve details in the apical region.