Snap evaporation of droplets on smooth topographies

Droplet evaporation on solid surfaces is important in many applications including printing, micro-patterning and cooling. While seemingly simple, the configuration of evaporating droplets on solids is difficult to predict and control. This is because evaporation typically proceeds as a “stick-slip” sequence—a combination of pinning and de-pinning events dominated by static friction or “pinning”, caused by microscopic surface roughness. Here we show how smooth, pinning-free, solid surfaces of non-planar topography promote a different process called snap evaporation. During snap evaporation a droplet follows a reproducible sequence of configurations, consisting of a quasi-static phase-change controlled by mass diffusion interrupted by out-of-equilibrium snaps. Snaps are triggered by bifurcations of the equilibrium droplet shape mediated by the underlying non-planar solid. Because the evolution of droplets during snap evaporation is controlled by a smooth topography, and not by surface roughness, our ideas can inspire programmable surfaces that manage liquids in heat- and mass-transfer applications.


Supplementary Note 1: Preparation of liquid-impregnated rough surfaces (LIRS)
Surface topographies were designed using standard 3D computer modelling software (Solid-Works). Samples of dimensions 5 cm × 5 cm × 0.5 cm were produced using an Objet30 3D printer in Vero White Plus RGD835 resin.
Planar-wave surfaces were modelled using a local elevation where x is a lateral coordinate (see Supplementary Figure 1a). Surfaces were produced for three different values of the amplitude, = {50, 100, 200} µm, and a fixed wavelength, λ = 2 mm.
Egg-box surfaces were modelled using a local elevation η(x, y) = cos 2πx Λ cos 2πy Λ ,  the surface using a micro-syringe. The droplet was allowed to evaporate within a closed chamber.
The humidity of the chamber was controlled by using hygroscopic silica beads. The temperature and relative humidity within the chamber were tracked using a digital sensor (Adafruit HTU21DF), and ranged between 20-25 • C, and 20-30%, respectively. Under such conditions, the typical evap- VIEW program at 1 frame every 300 seconds. The height and width of the droplet were measured by using the best fit of an ellipse to the apparent edge of the droplet (see Supplementary Figure 5b).
The experiment was repeated five times.

Supplementary Note 4: Lattice-Boltzmann simulations
Simulations were carried out using a binary-fluid lattice-Boltzmann algorithm. The geometry of the lattice is a cubic grid where lattice nodes are connected to their zeroth, first and second nearest neighbours (so-called D2Q9 and D3Q15 models for 2D and 3D simulations, respectively) [1]. Lattice nodes can be either "solid" or "fluid". At any given fluid node, indicated by a position vector r, we define two probability distribution functions, f i and g i , where the index i refers to the advection lattice propagation directions, c i . The time evolution of the distribution functions is given by the single-relaxation-time lattice-Boltzmann equations, These include a collision step where the distribution functions relax towards equilibrium values, indicated by the superscript "eq" over the respective relaxation timescales τ f and τ g (set to unity in the simulations), followed by a propagation step A suitable choice of the equilibrium distribution functions recovers the macroscopic equations of motion of the fluid in the limit of small Mach numbers. This is done such that the first moments read, i f eq i = ρ, i c i f eq i = ρv, i c i c i f eq i = P + ρvv, i g eq i = φ, i c i g eq i = φv, and i c i c i g eq i = M µ/(τ g − 1/2) + φvv, where P is the pressure tensor, and µ is the chemical potential [2]. After the collision step is performed, external forces can be applied to include gravitational effects.
We consider two types of boundary conditions: no-slip boundary conditions, applied at the solidfluid interface, and open boundaries that allow evaporation fluxes exit the simulation domain. The no-slip boundary conditions are enforced using a bounce-back algorithm [3] enhanced with a linear interpolation scheme to model the curved boundaries [4]. This ensures that no flows permeate through the solid. To model the wetting properties of the fluids, the derivatives of the phase-field are subject to a Neumann boundary condition [5]. Explicitly, the boundary conditions at any point on the solid surface (denoted by the position vector r w ) read ρv(r w ) = 0, n · ∇µ(r w ) = 0, where n is a vector normal to the solid, h is a parameter that controls the equilibrium contact angle, γ is the fluid-fluid surface tension and ξ is the fluid-fluid interface thickness.
Supplementary Equations (3)-(5) allow the pinning-free motion of the contact lines by means of diffusive fluxes of the chemical potential along the solid surface. Such a mechanism produces an apparent slip-length and a dynamic contact angle that deviates from its equilibrium value [6][7][8], in a manner consistent with contact-line hydrodynamics [9,10].
We set the simulation parameters, in lattice-Boltzmann units, to the following values: density ρ = 1, viscosity ν = 1/6, surface tension γ = 10 −3 , mean contact angle θ e = 105 • , contact angle variance σ(θ e ) = 4 • for the chemical noise in the solid surface, and the value of the chemical potential at the open boundaries µ b = −10 −5 , which is sufficient to drive the evaporation of the droplets in the diffusion-limit regime [11]. The 3D simulation domain is contained in a box of size 120 × 120 × 80, the mobility M = 1, interface thickness ξ = 4, and gravitational acceleration g = −2 × 10 −6ẑ . The 2D simulations were carried out in a box of size 240 × 120, and parameters: M = 4, ξ = 6, and g = 0. As an initial condition for all simulations, a spherical droplet of radius 1.5λ was centred at either a peak or a valley. The simulations were run for ∼ 10 7 iterations, which was a sufficiently long simulation time for the droplet to evaporate completely from the solid surface. We consider a two-dimensional (2D) sessile droplet on a smooth solid substrate which has a sinusoidal topography given by: where is the amplitude, and k = 2π/λ with λ the wavelength. In equilibrium, the shape of the droplet, h(x), is given by the Young-Laplace equation: where y(x) = h(x) + η(x), ∆p is the pressure difference across the liquid-gas interface and γ is its surface tension. For a droplet located either on a valley (with midpointx = x V = 0, ±λ, ±2λ, . . . ) or a peak (x = x P = ± 1 2 λ, ± 3 2 λ, . . . ), whose contact points are at x =x ± R, where R is the base radius, the apparent contact angle is [12]: where θ is the contact angle, which at equilibrium corresponds to the equilibrium contact angle θ e .
The area A of the droplet is then given by 2θ a − sin 2θ a sin 2 θ a + 2 R cos kx sin kR kR − cos kR .
For a given area, we numerically solve the above equation and find the different equilibrium configurations with θ = θ e . Supplementary Figures 6a , b and c show examples of the different droplet shapes found at a valley (x = 0) for A/λ 2 = 2, 1.5 and 1.35. In all cases, the equilibrium contact angle is θ e = 105 • .

Local stability analysis
We study the local stability of the equilibrium solutions by analysing the interfacial energy which is given by: By locally perturbing the equilibrium points with axisymmetric perturbations (R → R + δR), and asymmetric perturbations (x →x + δx), we can construct, for the fixed area given by Supplementary Equation (12), the energy landscape F (R,x). We find that the equilibrium configurations correspond to extrema of the energy function, for which ∂xF = ∂ R F = 0, and can be a local minimum (stable to both δR and δx), a local maximum (unstable to both δR and δx), or a saddle node (stable to δR perturbations but unstable to δx perturbations). Supplementary Figures 6d, e, and f show plots of the energy as function of R/λ for the three cases of A considered in Supplementary   Figures 6a, b, and c, where we can see how equilibrium solutions emerge as A is varied.
Therefore, by continuously changing the droplet area we can construct the bifurcation diagrams for a droplet on a valley or peak of the solid surface. The results for different surface aspect ratios are shown in Supplementary Figure 7a-f. We can see that for small aspect ratios and small droplet areas, the stability of the equilibrium solutions alternates between saddle nodes (black solid lines) and stable points (blue/red solid lines for peak/valley solutions). The transition between these sates is a consequence of a pitchfork bifurcation where two additional saddle node branches emerge