# Small Leidenfrost droplet dynamics

Benjamin Sobac, Alexey Rednikov, Pierre Colinet

► To cite this version:

Benjamin Sobac, Alexey Rednikov, Pierre Colinet. Small Leidenfrost droplet dynamics. Journal of Fluid Mechanics, 2025, 1010, pp.A47. 10.1017/jfm.2025.208. hal-04797907

HAL Id: hal-04797907

https://hal.science/hal-04797907v1

Submitted on 19 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Small Leidenfrost droplet dynamics

Benjamin Sobac¹ ID, Alexey Rednikov² ID and Pierre Colinet² ID

¹CNRS, Université de Pau et des Pays de l'Adour, E2S UPPA, LFCR, 64600 Anglet, France

²TIPs Lab, Université libre de Bruxelles, 1050 Brussels, Belgium

**Corresponding author:** Benjamin Sobac, benjamin.sobac@univ-pau.fr

(Received 25 September 2024; revised 12 January 2025; accepted 16 February 2025)

An isolated Leidenfrost droplet levitating over its own vapour above a superheated flat substrate is considered theoretically, the superheating for water being up to several hundred degrees above the boiling temperature. The focus is on the limit of small, practically spherical droplets of several tens of micrometres or less. This may occur when the liquid is sprayed over a hot substrate, or just be a late life stage of an initially large Leidenfrost droplet. A rigorous numerically assisted analysis is carried out within verifiable assumptions such as quasi-stationarities and small Reynolds/Péclet numbers. It is considered that the droplet is surrounded by its pure vapour. Simple formulae approximating our numerical data for the forces and evaporation rates are preliminarily obtained, all respecting the asymptotic behaviours (also investigated) in the limits of small and large levitation heights. They are subsequently used within a system of ordinary differential equations to study the droplet dynamics and take-off (drastic height increase as the droplet vapourises). A previously known quasi-stationary inverse-square-root law for the droplet height as a function of its radius (at the root of the take-off) is recovered, although we point out different prefactors in the two limits. Deviations of a dynamic nature therefrom are uncovered as the droplet radius further decreases due to evaporation, improving the agreement with experiment. Furthermore, we reveal that, if initially large enough, the droplets vanish at a universal finite height (just dependent on the superheat and fluid properties). Scalings in various distinguished cases are obtained along the way.

**Key words:** boiling, drops, condensation/evaporation

## 1. Introduction

When a volatile liquid droplet is placed on a hot solid surface, superheated well above the boiling temperature, it neither touches the substrate nor boils, but rather floats on a thin
film of its own vapour. This fascinating phenomenon, known as the Leidenfrost effect, has not ceased to attract attention since its first descriptions approximately 300 years ago (Boerhaave 1732; Leidenfrost 1756). This is due to not only the myriad of intriguing and unexpected behaviours a droplet can exhibit in this state, but also its relevance across a wide range of industrial and technological processes, spanning from the traditional heat transfer applications to the emerging field of multiphase milli-/micro-fluidics. See, for example, the review articles of Quéré (2013), Ajaev & Kabov (2021) and Stewart (2022), dedicated book chapters in Brutin (2015) and Marengo & De Coninck (2022) and the many references therein.

The vapour film, a key feature of the Leidenfrost state, ensures the droplet levitation while acting as a thermal insulator, resulting in relatively low evaporation rates and hence long lifetimes of the droplet. In this state, the droplet’s weight is balanced by the pressure within such a vapour cushion squeezed by the slowly and steadily evaporating droplet. With no contact with the substrate, the observable shapes of the droplets are governed by a balance between capillarity and gravity similarly to a perfectly non-wetting (superhydrophobic) situation. Denoting the capillary length by $ℓ_c$, droplets with radii $R$ smaller than $ℓ_c$ remain quasi-spherical while puddles larger than $ℓ_c$ are flattened by gravity, whose height is limited by $≈ 2ℓ_c$ (Biance et al. 2003). The profile of the underlying vapour film is non-trivial. For a large droplet with $R \gtrsim ℓ_c$, the vapour film exhibits a pocket-like structure composed of an internal vapour ‘pocket’ surrounded by a thin neck. As the droplet gets smaller, the vapour film slims down, the droplet getting closer to the substrate. When the droplet radius is small enough as compared with $ℓ_c$, the vapour pocket disappears completely, and the droplet becomes quasi-spherical, with a small circular area slightly flattened at the bottom. Accurate interferometric measurements of the vapour film thickness profile (Burton et al. 2012) turn out to be in a good agreement with a refined theoretical model (Sobac et al. 2014, 2021) coupling lubrication vapour flow, capillarity and hydrostatic pressure effects, which was recently confirmed numerically by Chakraborty et al. (2022). Note that the main scaling laws featuring the shapes of a Leidenfrost droplet and its evaporation dynamics can be found in Biance et al. (2003), Pomeau et al. (2012) and Sobac et al. (2014, 2021).

In practice, this absence of contact between the Leidenfrost droplet and the substrate leads to very rich dynamics. For large puddle-like droplets, the vapour pocket grows until it eventually pops up as a central ‘chimney’ due to a Rayleigh-Taylor mechanism (Biance et al. 2003; Snoeijer et al. 2009). Instability of large droplets can also occur (either spontaneously or forced) in the form of ‘star-faceted’ shapes when azimuthal surface oscillations develop along the periphery of the droplets (Brunet & Snoeijer 2011; Ma et al. 2017; Ma & Burton 2018; Bergen et al. 2019; Bouillant et al. 2021a). Self-induced spontaneous oscillations can also occur in the vertical plane, yielding the recently reported bobing, bouncing or trampolining dynamics when Leidenfrost droplets reach moderate and small sizes with $R \le ℓ_c$ (Liu & Tran 2020; Graeber et al. 2021). Other spectacular behaviours related to their high mobility have been observed. These include Leidenfrost wheels, when a droplet initially at rest spontaneously rolls and moves over a flat surface like a wheel due to symmetry breaking in the internal flow of the liquid (Bouillant et al. 2018), and self-propelling of Leidenfrost droplets when interacting with substrates breaking the axisymmetry, either due to surface topography such as ratchets or herringbones (Linke et al. 2006; Dupeux et al. 2011; Marín et al. 2012; Baier et al. 2013; Soto et al. 2016; Dodd et al. 2019), or temperature gradients (Sobac et al. 2017; Dodd et al. 2020; Bouillant et al. 2021b). Thus, droplets move in a direction dictated by the patterns due to symmetry breaking of the vapour layer. Strategies have emerged to control the motion and manipulate these droplets. In addition to geometric and thermal
heterogeneities, chemical patterns of the surface can also be exploited to tailor the vapour film, enabling the stretching, sloshing, spinning, propelling or trapping of a Leidenfrost droplet (Li et al. 2023).

As compared with large and moderate size droplets, the dynamics of small, near-spherical Leidenfrost droplets ($R \ll l_c$) has not received so much attention. In their seminal work, Celestini et al. (2012) first explored the final fate of Leidenfrost droplets as they became very small, just moments before disappearing. By spraying tiny droplets of water or ethanol in the size range of approximately 1–30 µm onto a superheated substrate, they discovered that when small enough (i.e. with $R$ below a characteristic radius corresponding to the breakup of the lubrication approximation), Leidenfrost droplets took off from the heated substrate with an elevation $h \propto R^{-1/2}$, as predicted also by Pomeau et al. (2012). Remarkably, in this regime, droplets become too light to withstand the upward force generated due to evaporation, and they reach higher and higher elevations while vapourising. This behaviour drastically contrasts with what is observed for larger Leidenfrost droplets. More recently, Lyu et al. (2019) observed that a second final fate, other than lift-off, is possible for Leidenfrost droplets. Namely, if the liquid droplet is not pure or contaminant-free enough, small Leidenfrost droplets are unable to take off, but instead disappear by exploding with an audible crack.

Here, we propose to theoretically revisit the dynamics of small spherical Leidenfrost droplets with the aim of comprehensively and thoroughly analysing the mechanisms involved in their final fate. Thanks to a model including a realistic description of the coupling between hydrodynamics, heat transfer and evaporation, this work seems to be the first to provide exact estimates of the droplet elevation as a function of the physical parameters and without any fitting parameter. After numerically computing the entirety of fluxes, evaporation rates and forces, a master curve for droplet elevation as a function of its size is derived by simply balancing the droplet weight with the upward evaporation-induced hydrodynamic force. While the scaling law agrees with Celestini et al. (2012) and Pomeau et al. (2012), there appear subtleties concerning the prefactor. Moreover, the analysis reveals that such a classical quasi-steady description is not fully sufficient to describe the take-off phenomenon. Even at these small scales, further dynamical effects must be taken into account to achieve a good agreement with the original experimental data of Celestini et al. (2012).

## 2. Statement of the problem, premises and outlook

Consider a small evaporating spherical droplet of radius $R$ in a Leidenfrost state levitating at a height $h$ above a superheated substrate at a ‘wall’ temperature $T_w$, as sketched in figure 1. The substrate is flat and horizontal. We shall be interested in the take-off of Leidenfrost droplets (Celestini et al. 2012; Pomeau et al. 2012), which occurs in the realm of small droplets ($R$ of the order of tens of µm) with a negligible deviation from the spherical shape. The (immediate) surroundings of the droplet are assumed to be saturated with vapour (totally displacing the air) and heated through to the substrate temperature. Thus, $T_w$ is here also an effective overall ambient temperature (i.e. $T_{amb} = T_w$). The droplet is assumed isothermal at saturation (boiling) temperature $T_{sat}$. The superheat is given by $\Delta T = T_w - T_{sat} > 0$. The typical values considered here are $T_{sat} = 100$ °C (water at atmospheric pressure), $T_w = 400$ °C and $\Delta T = 300$ °C. Such superheat occurs e.g. in the experiments by Celestini et al. (2012).

Mathematically, the goal of the present consideration is obtaining an interrelation between $h$, $R$ (in particular, as functions of time $t$ due to the droplet evaporation) and the parameters of the problem, such as $\Delta T$, $g$ (gravitational acceleration, 9.81 m s$^{-2}$)
<figure><img src="image_3.png" /><figcaption>Figure 1. Sketch of the problem.</figcaption></figure>

and the liquid and vapour properties. The latter include $\rho_l$ (liquid density, defined at $T_{sat}$) as well as the following vapour properties: $\rho_v$ (density), $\mu_v$ (dynamic viscosity), $\nu_v = \mu_v/\rho_v$ (kinematic viscosity), $\lambda_v$ (thermal conductivity), $c_{p,v}$ (heat capacity at constant pressure), $\alpha_v = \lambda_v/(\rho_v c_{p,v})$ (thermal diffusivity) and $\mathcal{L}$ (latent heat of evaporation). These may vary considerably in the temperature range between $T_{sat}$ and $T_w$. However, for the sake of simplicity, we shall here assume them constant and defined at the mid-temperature $(1/2)(T_w + T_{sat})$, except for $\rho_l$ and $\mathcal{L}$ defined at $T_{sat}$, similarly to the approach used elsewhere (Sobac et al. 2014, 2015). The relevant property values are provided in Appendix A.

Other key assumptions include negligible advective/convective effects (small Péclet and Reynolds numbers), so that the temperature field in the vapour is governed by heat conduction, while the evaporative flow from the droplet can be considered by means of the Stokes approximation. Quasi-steadiness of the temperature and velocity fields, in spite of $R$ and $h$ changing in time due to evaporation, is another key assumption of the analysis. In other words, these fields and the evaporation fluxes and forces they determine, are merely functions of the instantaneous values of $R$ and $h$ and do not depend on the history. It is under this premise that a preliminary calculation of these quantities is carried out in § 3. The validity of this and other assumptions is verified a posteriori in their due course.

The quasi-steadiness assumption is also applied at first when it comes to the force balance on the droplet in § 4, permitting to predict the levitation height $h$ and make a first comparison with experiment. Yet, certain limitations are thereby disclosed, inspiring consideration of a more general droplet dynamics in §§ 5–7. However, even in such a situation, the quasi-steadiness of the quantities such as those calculated in § 3 is still assumed to hold.

An important geometric parameter of the configuration (figure 1), meriting a special notation, is the ratio of the droplet’s height and radius

$$ \delta = \frac{h}{R}. \qquad (2.1) $$

In the present study, we shall be interested in a full range of this relative-height parameter, from very small to very large. The large values are expected for small droplets (small $R$) upon a take-off (Celestini et al. 2012; Pomeau et al. 2012; Sobac et al. 2021; Chakraborty
$$
\begin{array}{ccccccc}
[r, z, s] & [S] & [T - T_{sat}] & [j] & [J] & [\mathbf{v}] & [p] & [F_{ev}] \\
R & R^2 & \Delta T & \frac{\lambda_v \Delta T}{L^2 R} & [j] R^2 = \frac{\lambda_v \Delta T R}{L} & \frac{[j]}{\rho_v} = \frac{\lambda_v \Delta T}{\rho_v L^2 R} & \frac{\mu_v [\mathbf{v}]}{R} & \mu_v R [\mathbf{v}] = \frac{\mu_v \lambda_v \Delta T}{\rho_v L}
\end{array}
$$

<caption>Table 1. Scales used to render the various quantities dimensionless in § 3: cylindrical coordinates r and z, arclength s, surface area S, temperature difference T − Tsat, evaporation flux density j (kg m−2 s−1), evaporation rate J (kg s−1), evaporative velocity v and pressure p fields and evaporation/levitation force Fev. The square brackets denote the scale of the quantity inside.</caption>

et al. 2022). In contrast, small δ are attained for larger droplets. In this way, we arrive
at a transition from spherical Leidenfrost droplets to Leidenfrost droplets for which
deformation (at first at the bottom slightly flattened by gravity) becomes essential. Such a
transition is touched upon in § 8.

3. Basic calculations: fields, fluxes and forces

Dimensionless variables are introduced using the scales given in table 1 (definitions to be given in their due course). For simplicity and expecting no confusion, no notation distinction is made between the original, dimensional variables and their dimensionless versions in the present section (the distinction being clear from the context). We just note that a dimensionless temperature is introduced as

$$
\bar{T} = \frac{T - T_{sat}}{\Delta T}, \qquad (3.1)
$$

where recall that $\Delta T = T_w - T_{sat}$. Hereafter, in the same spirit, the bars are omitted for
the sake of brevity.

3.1. *Temperature field*

As stipulated in § 2, the heat transport in the gas phase is conductive and quasi-steady.
Thus, the thermal problem is decoupled from the evaporative velocity field, and the
dimensionless temperature field *T* is governed by the Laplace equation

$$
\nabla^2 T = 0. \tag{3.2}
$$

It is subject to the boundary conditions

$$
T = 1 \quad \text{on the hot substrate,} \tag{3.3}
$$

$$
T \to 1 \quad \text{far away from the droplet,} \tag{3.4}
$$

$$
T = 0 \quad \text{on the droplet surface.} \tag{3.5}
$$

Although an exact solution in bipolar coordinates can be found e.g. using the methods by Lebedev (1972), it is rather cumbersome so that we eventually opt for a numerical solution using COMSOL Multiphysics.

The results of the simulations are shown in figure 2(a) for three different values of the separating distance δ: a large droplet very close to the substrate with δ = 0.1; a droplet at a distance from the substrate comparable to its radius with δ = 1; and a small droplet beginning to be far away from the substrate with δ = 5. One immediately observes that, at small δ, the temperature difference is squeezed into a thin film between the droplet and the substrate. At large δ, the temperature field approaches a spherically symmetric one, as expected. Other results displayed in figure 2 will be discussed later.
<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Dimensionless temperature (right) and velocity (left) fields for various values of $\delta$. Streamlines are also shown by white lines (left). (b) Corresponding profiles of the dimensionless evaporation flux $j$ along the droplet surface as a function of the dimensionless arclength $s$ ($s=0$ at the droplet apex, $s=\pi$ at its lowest point; the plot is formally continued up to $s=2\pi$ back to the apex for aesthetic purposes). Insets are simply zooms of the main plot.</figcaption></figure>

## 3.2. Evaporation rate

At the droplet surface, the evaporation mass flux $j$ [$kg m^{-2} s^{-1}$] is at the expense of the heat coming from the superheated surroundings through the vapour phase: $j = (\lambda_v / \mathcal{L}) n \cdot \nabla T$, where $n$ is the (external) unit normal vector. In dimensionless terms (cf. table 1), this reads

$$j = n \cdot \nabla T. \qquad (3.6)$$

Using the temperature field computed in § 3.1, the profiles of the evaporation flux $j$ are calculated from (3.6) and shown in figure 2(b). One can appreciate that, due to the presence of the hot substrate, $j$ is maximum at the base and decreases towards the apex, where the minimum is attained. The closer the droplet is to the substrate (the smaller $\delta$ is), the more the profile of $j$ is non-uniform and the values of $j$ are large. At small $\delta$, one obviously obtains $\max(j) \propto 1/\delta$ (heat conduction across a thin vapour layer). When the relative droplet height $\delta$ increases, the non-uniformity of $j$ weakens and the average of $j$ decreases. Eventually, $j$ tends to a uniform value of 1 for $\delta \gg 1$ (as for a droplet in an unbounded medium).

The (global) evaporation rate $J$ can be directly deduced by integrating the evaporation flux all over the droplet surface

$$J = \iint j \, dS. \qquad (3.7)$$

Figure 3 reports the computed values of $J$ as a function of $\delta$. As expected from the knowledge of the $j$ behaviour, $J$ diverges as $\delta \to 0$ and decreases to saturate at $4\pi$ as
<figure><img src="image_5.png" /><figcaption>Figure 3. Evaporation rate $J$ (dimensionless) as a function of the relative droplet height $\delta$. Numerical results (blue open circles) are fairly well approximated by (3.8) (red line). Black dashed lines correspond to the asymptotic behaviours (Appendix C).</figcaption></figure>

$\delta \to +\infty$. Such asymptotic behaviours are investigated in detail in Appendices C1 and C2 and also represented in figure 3. A good simple approximation of the numerical data for $J$, respecting the leading-order asymptotic behaviours, is given by

$$J(\delta) = 4\pi \left[ 1 + \frac{1}{2} \ln \left( 1 + \frac{1}{\delta} \right) \right], \quad (3.8)$$

where a maximum deviation from the data does not exceed 2.7%. However, a more precise fit is also provided for reference in Appendix B1.

### 3.3. Velocity field

In accordance with the approach followed in the present paper (§ 2), the evaporative flow field, generated by droplet evaporation, is considered in the Stokes and quasi-steady approximation. Thus, we proceed from the continuity and Stokes equations

$$\nabla \cdot \mathbf{v} = 0, \quad (3.9)$$

$$\nabla^2 \mathbf{v} - \nabla p = 0. \quad (3.10)$$

The following boundary conditions are used:

$$\mathbf{v} = 0 \quad \text{on the hot substrate,} \quad (3.11)$$

$$\mathbf{v} \to 0 \quad \text{far away from the droplet,} \quad (3.12)$$

$$\mathbf{v} \cdot \boldsymbol{\tau} = 0 \quad \text{and} \quad \mathbf{v} \cdot \mathbf{n} = j \quad \text{on the droplet surface.} \quad (3.13)$$

Here, $\mathbf{v}$ and $p$ are the dimensionless velocity and pressure fields in the vapour (cf. the scales in table 1), and $\boldsymbol{\tau}$ is the unit tangential vector. In the first boundary condition (3.13), a possible internal flow in the droplet is neglected relative to the velocity scale in the vapour, hence no slip. The last boundary condition (3.13) contains the driving factor of the flow field, where the normal velocity at the droplet surface is determined by the evaporation flux (3.6), which is in turn determined by the temperature field obtained from the formulation (3.2)–(3.5) in § 3.1. Similarly to § 3.1, this hydrodynamic part of the problem is also solved numerically using COMSOL Multiphysics.
<figure><img src="image_6.png" /><figcaption>Figure 4. Evaporative force $F_{ev}$ (dimensionless, in terms of $\delta^2 F_{ev}$) as a function of the relative droplet height $\delta$. Numerical results (blue open circles) are well approximated by (3.15) (red line). Black dashed lines correspond to the asymptotic behaviours.</figcaption></figure>

The computation results for the vapour velocity fields are added into figure 2(a), mirroring the temperature fields and also displaying the streamlines. For large relative heights $\delta$, the streamlines remain straight near the droplet’s surface, indicating that the flow is almost spherically symmetric there and only slightly disturbed by the substrate. Farther from the droplet, the streamlines are significantly bent due to the substrate presence. Higher velocity field values are attained for smaller $\delta$. This is not only due to a profound maximum in the evaporation flux due to the substrate proximity at small $\delta$ (as in figure 2b1), but also additionally due to a confinement effect in a thin vapour layer between the droplet and the substrate, when the longitudinal velocity becomes even higher than the evaporation-flux-driven normal one at the droplet surface.

### 3.4. Levitation force

The bending and asymmetry of the evaporative flow due to the substrate give rise to a hydrodynamic force acting on the evaporating droplet in the sense of its repulsion from the substrate. We refer to it as an evaporative force $F_{ev}$. In our configuration (figure 1), this amounts to a force acting on the droplet vertically upwards (along the z axis), which is responsible for droplet levitation against gravity. The force balance on the droplet and its levitation height are considered in § 4 later on. Here we simply calculate $F_{ev}$ in dimensionless terms (cf. the scales in table 1) as a function of the relative height $\delta$. Namely, we evaluate

$$F_{ev} = \left( \iint_S (-p\mathbf{n} + (\nabla\mathbf{v} + \nabla\mathbf{v}^T) \cdot \mathbf{n}) \, dS \right) \cdot \mathbf{e}_z \quad (3.14)$$

using the velocity and pressure fields computed in § 3.3, where $\mathbf{e}_z$ is a unit vector along z.

The result is reported in figure 4 in terms of $F_{ev}\delta^2$. The overall tendency is $F_{ev} \propto \delta^{-2}$, as is already known in the literature (Celestini et al. 2012; Pomeau et al. 2012). However, it is less known that the prefactor is different in the limits $\delta \to 0$ and $\delta \to +\infty$. For instance, Celestini et al. (2012) attempted to fit the experimental data using a single prefactor. We obtain a prefactor $3\pi$ as $\delta \to 0$, which can be deduced from the lubrication approximation (cf. Pomeau et al. 2012; Sobac et al. 2021, see also Appendix C3), although Pomeau et al. (2012) obtained $3\pi/8$ here (erroneously, in our opinion). In contrast, the prefactor is $6\pi$
as $\delta \to +\infty$, which is confirmed by an asymptotic analysis described in Appendix C4, where a number of contributions in terms of the droplet-substrate interaction are followed through. The following simple expression nicely approximates our numerical result while respecting the prefactor values in both limits:

$$F_{ev}(\delta) = \frac{3\pi}{\delta^2} \frac{1 + 2\delta}{1 + \delta}, \qquad (3.15)$$

which is also shown in figure 4. It covers the numerical data with a relative error of 1.4%, whereas a more precise fit is provided for reference in Appendix B2. In what follows, for purely presentational reasons, we shall use the simple-form approximations like (3.8) and (3.15), whose precision is deemed to be already rather satisfactory.

### 3.5. Validity of assumptions

1. Negligible advection. We start with the estimation of an evaporative Péclet number $Pe = [v]R/\alpha_v$, where the evaporative velocity scale $[v]$ from table 1 is used. We obtain

$$Pe = \frac{c_{p,v} \Delta T}{\mathcal{L}}, \qquad (3.16)$$

which incidentally turns out to be a version of the Jakob number often used in the literature. For the typical $\Delta T$ value (cf. § 2) and parameter values (cf. the first row of table 3 in Appendix A), we obtain $Pe \approx 0.19 \ll 1$, which justifies the approximation used in § 3.1.

2. Stokes approximation. Likewise, the Reynolds number is $Re = \rho_v R[v]/\mu_v = Pe Pr^{-1}$, where $Pr = \mu_v/(\rho_v \alpha_v)$ is the Prandtl number. As $Pr = 0.71$ (cf. ibid), we obtain $Re \approx 0.26 \ll 1$, confirming the approximation used in § 3.3.

3. Negligible natural convection. This is related to small values of the Grashof number at the droplet scale

$$Gr = \frac{\rho_v g R^3}{\mu_v \nu_v}; \qquad (3.17)$$

(written in this form given that the variations of $\rho_v$ are here of the order of $\rho_v$ itself). One typically obtains $Gr < 0.01$ for our small droplets ($R \lesssim 50$ µm).

4. Gas phase quasi-steadiness. The results of the present section imply the quasi-steadiness of the temperature and velocity fields in the entire region between the droplet and the substrate, which may be especially questionable for large levitation heights $h$. The appropriate thermal and viscous time scales can be chosen as $\tau_{th} = \max(R, h)^2/\alpha_v$ and $\tau_{vis} = \rho_v \max(R, h)^2/\mu_v$. The quasi-steadiness takes place when $\tau_{th} \ll \tau$ and $\tau_{vis} \ll \tau$, where $\tau$ is the typical time scale of the process. As $Pr = O(1)$ here, we just limit our attention to the first one of these conditions for the sake of brevity, and hence

$$\frac{\max(R, h)^2}{\alpha_v} \ll \tau. \qquad (3.18)$$

An immediately obvious time scale of the process is here the evaporation time scale of the droplet $\tau_{ev} = \rho_l R^3/[J]$ (cf. table 1 for [J]), i.e.

$$\tau_{ev} = \frac{\rho_l \mathcal{L} R^2}{\lambda_v \Delta T} = \frac{\rho_l}{\rho_v} \frac{1}{Pe} \frac{R^2}{\alpha_v}. \qquad (3.19)$$

Using it as $\tau = \tau_{ev}$ in (3.18), we arrive at

$$\max(1, \delta^2) \ll \frac{\rho_l}{\rho_v Pe}. \qquad (3.20)$$
Given that $\rho_l \gg \rho_v$ and $Pe \ll 1$ here, the condition (3.20) leaves a considerable margin for possible large values of the relative height $\delta = h/R$. We shall come back to it later on, after having considered concrete solutions for $h$.

## 4. Levitation and take-off: quasi-steady approach

### 4.1. Formal application of the quasi-steady approach

A vertical balance of the (evaporative) levitation force (3.15) against the droplet weight directly yields the equation for the levitation height of the droplet $h$ as a function of its radius $R$, which we write in dimensional form:

$$3\pi \frac{\mu_v \lambda_v \Delta T}{\rho_v \mathcal{L}} \frac{R^2}{h^2} \frac{R+2h}{R+h} = \frac{4\pi}{3} \rho_l g R^3, \quad (4.1)$$

where we have used the scale from table 1 and the definition (2.1). There exists a single natural length scale $\ell_*$, the non-dimensionalisation with which renders (4.1) parameter-free:

$$\frac{1}{\hat{h}^2} \frac{\hat{R} + 2\hat{h}}{\hat{R} + \hat{h}} = \frac{4}{9} \hat{R}, \quad (4.2)$$

where

$$[\hat{R}] = [\hat{h}] = \left( \frac{\mu_v \lambda_v \Delta T}{\rho_v \rho_l g \mathcal{L}} \right)^{1/3} \equiv \ell_*, \quad \hat{R} = \frac{R}{[\hat{R}]}, \quad \hat{h} = \frac{h}{[\hat{h}]}. \quad (4.3)$$

The solution of (4.2) is shown in figure 5(a) and adheres to the following asymptotic behaviour:

$$\hat{h} = \frac{3}{2} \frac{1}{\hat{R}^{1/2}} \quad \text{i.e.} \quad \delta = \frac{3}{2} \frac{1}{\hat{R}^{3/2}} \quad \text{as} \quad \hat{R} \to +\infty, \quad (4.4)$$

$$\hat{h} = \frac{3}{\sqrt{2}} \frac{1}{\hat{R}^{1/2}} \quad \text{i.e.} \quad \delta = \frac{3}{\sqrt{2}} \frac{1}{\hat{R}^{3/2}} \quad \text{as} \quad \hat{R} \to 0. \quad (4.5)$$

The length scale $\ell_*$ indicates the characteristic size $R \sim \ell_*$ at which the droplet takes off at a height of the order of itself, with $h \sim R$ (i.e. $\delta \sim 1$). At smaller sizes ($R \ll \ell_*$), the droplet soars even higher ($h \gg \ell_* \gg R$, $\delta \gg 1$), whereas at larger sizes ($R \gg \ell_*$), the droplet levitates lower ($h \ll \ell_* \ll R$, $\delta \ll 1$). Typically, $\ell_*$ is in the range of a few tens of micrometres. For our reference case of a water droplet on a superheated substrate with $\Delta T = 300 \,^\circ\text{C}$, we obtain $\ell_* = 28.46 \, \mu\text{m}$, which is much smaller than the capillary length $\ell_c = \sqrt{\gamma/(\rho_l g)}$ ($\gamma$ being the liquid-air surface tension, $\ell_c \sim 2.5 \, \text{mm}$ for water). This ‘take-off scale’ $\ell_*$ has earlier been pointed out by Celestini et al. (2012) and Pomeau et al. (2012) (their notation $R_l$) as the scale at which a drastic take-off takes place (although note that a mere increase of $h$ as $R$ decreases already starts from much larger sizes $R$, cf. Sobac et al. 2021, as well as § 8 here). They also interpret it as the droplet size starting from and below which the lubrication approximation in the vapour film between the droplet and the substrate becomes invalid (since $h \ll R$ ceases indeed to hold). Note that $h = R$ for $R = 3/2 \ell_*$.

Similarly to what was commented for $F_{ev}$ in § 3.4, the overall tendency $h \sim R^{-1/2}$ (or $h/R \sim R^{-3/2}$) has been well known since Celestini et al. (2012) and Pomeau et al. (2012), who first pointed out this exponent. However, we here calculate the prefactor and point out that it is actually not fully constant, as highlighted in the inset of figure 5(a) and further put into evidence by the two different limiting values in (4.4) and (4.5).
<figure><img src="image_7.png" /><figcaption>Figure 5. (a) Relative height as a function of the droplet radius as predicted by the quasi-steady model. Full solution (solid blue line, 'master curve') and the asymptotic behaviours for smaller and larger radii (dashed and dotted purple lines, respectively). The log-log inset highlights the dominant power law and the prefactor change between the two limits. (b) First comparison with the experimental data of Celestini et al. (2012) for water (with $\Delta T = 300$ °C).</figcaption></figure>

A first comparison with experiment is undertaken in figure 5(b), where we focus our attention on just the upper layer of experimental points, while the points lower than that are deemed to belong to some transients (see also a remark in § 7 later). It is important to recall that these experiments dealt with small droplets, typically ranging in radius from 1 to 30 µm, which is of the order of or smaller than $\ell_*$ here. For this size range ($\hat{R} \lesssim 1$), specific to take-off observations, the full solution of (4.2) already practically coincides with the asymptotic limit (4.5), cf. figure 5. While this seems to agree well with experiment for 15 µm $\lesssim R \lesssim 30$ µm, an overprediction is nonetheless observed as the droplet size decreases below $R \lesssim 15$ µm. Astonishingly, we observe that it is rather the asymptotic behaviour (4.4) that starts to get closer to the experimental points even if (4.4), with the prefactor it contains, is appropriate in the limit of larger droplets (and not the smaller ones we are discussing right now). We shall come later to what will be the right explanation here.

The droplet radius $R$ decreases over time by evaporation (rather than just being a given constant parameter), and hence $h$, related to $R$ by (4.1), is also a function of time. The steady force balance (4.1) or (4.2) is then assumed to be valid in a quasi-steady sense, and $R(t)$ and $h(t)$ follow the solid curve of figure 5 as time goes on. The mass of the droplet $(4\pi/3)\rho_l R^3$ decreases in time at an evaporation rate given by (3.8) with the scale from table 1. This balance gives rise to the following equation for the droplet radius evolution:

$$ \rho_l R \frac{dR}{dt} = - \frac{\lambda_v \Delta T}{\mathcal{L}} \left[ 1 + \frac{1}{2} \ln \left( 1 + \frac{R}{h} \right) \right]. \quad (4.6) $$

Using the time scale

$$ [\hat{t}] = \frac{\rho_l \mathcal{L} \ell_*^2}{\lambda_v \Delta T} = \frac{\rho_l}{\rho_v} \frac{1}{Pe} \frac{\ell_*^2}{\alpha_v} = \tau_*, \quad (4.7) $$

alongside the scales (4.3), equation (4.6) is rendered free of any parameters similarly to (4.2). Namely, we arrive at

$$ \hat{R} \frac{d\hat{R}}{d\hat{t}} = -1 - \frac{1}{2} \ln \left( 1 + \frac{\hat{R}}{\hat{h}} \right). \quad (4.8) $$
<figure><img src="image_8.png" /><figcaption>Figure 6. Dimensionless evolutions of the radius $\hat{R}$ and the height $\hat{h}$ of a spherical Leidenfrost droplet over time $\hat{t}$ computed by the quasi-steady model (coupled 4.2 and 4.6) for initial droplet radii $\hat{R}_0 = \{1/3, 1/2, 1, 2, 3\}$ (from dark to light blue). The corresponding initial quasi-steady heights are $\hat{h}_0 = \{3.60, 2.89, 1.93, 1.25, 0.97\}$. The inset serves to illustrate the extents to which the $R^2$-law holds and to which the evaporation is accelerated by the superheated substrate, where the time is normalised to the evaporation time of a freely suspended droplet.</figcaption></figure>

Now the dimensionless evolution problem for $\hat{R}(\hat{t})$ and $\hat{h}(\hat{t})$ is defined by a system of two coupled equations, (4.2) and (4.8), for which the initial conditions $\hat{R} = \hat{R}_0$ and $\hat{h} = \hat{h}_0$ at $\hat{t} = 0$ are posed with $\hat{R}_0$ and $\hat{h}_0$ not being independent but rather related by (4.2). Figure 6 illustrates the (numerically obtained) solution for various initial droplet radii $\hat{R}_0 = \{1/3, 1/2, 1, 2, 3\}$. Evidently, the curves demonstrate that $\hat{R}$ decreases over time due to evaporation until extinction, with larger droplets exhibiting longer lifespans. Concurrently, $\hat{h}$ increases over time as the droplet size decreases, larger droplets being closer to the substrate at the initial time in accordance with (4.2). It is important to note that, within the present quasi-steady description, a Leidenfrost droplet takes off reaching an infinite height $\hat{h} \to +\infty$ at the end of its life (as $\hat{R} \to 0$), in accordance with (4.5).

In the inset, $(\hat{R}/\hat{R}_0)^2$ is plotted as a function of $\hat{t}/\hat{t}_{ev}^\infty$ in order to highlight the evaporative behaviour of a spherical Leidenfrost droplet as compared with the well-known limit case of a spherical droplet suspended in an unbounded gas medium. Here, $\hat{t}_{ev}^\infty = (\hat{R}_0)^2/2$ is the dimensionless evaporation time of such a freely suspended droplet (which can be derived from (4.8) in the limit $\delta = \hat{h}/\hat{R} \to +\infty$). Owing to the interaction with the superheated substrate, appearing through the logarithmic term in (4.8), the well-known $R^2$-law is recovered only for large values of $\hat{h}_0$. Thus, $\hat{R}^2$ generally does not linearly decrease in time, while these droplets evaporate faster than their freely suspended counterparts due to the proximity of the superheated substrate.

Needless to note that, by parametrically plotting $\hat{h}(\hat{t})/\hat{R}(\hat{t})$ or $\hat{h}(\hat{t})$ as a function of $\hat{R}(\hat{t})$ ($\hat{t}$ being the parameter), we retrieve the same ‘master’ curve as depicted in figure 5. Within the present quasi-steady approach, such a master curve is just trivially given by an algebraic equation like (4.1) or (4.2). However, this may become less trivial in what follows.

## 4.2. Validity of assumptions

The quasi-steady result that the droplet soars to an infinite height at the end of its life, cf. (4.5) and figure 6(b), looks suspicious from the physical point of view. One can wonder whether the quasi-steadiness criterion (3.20) for the fields in the gas phase is still fulfilled
in view of $\delta \to +\infty$. Furthermore, $h \to +\infty$ also implies infinite velocity and acceleration ($dh/dt \to +\infty$ and $d^2h/dt^2 \to +\infty$), and hence one can wonder whether the steady force balance such as (4.1) is still adequate when neglecting the drag (proportional to $dh/dt$) and inertia (proportional to $d^2h/dt^2$) forces on the droplet. Below, we make an estimation of those effects when the only source of unsteadiness is the evaporation of the droplet, i.e. when the time scale is $\tau = \tau_{ev}$ as given by (3.19).

Disregarding numerical prefactors, the inertia, (Stokes) drag and levitation forces can be estimated as

$$ \text{inertia} \sim \rho_l R^3 \frac{h}{\tau_{ev}^2}, \quad \text{drag} \sim \mu_v R \frac{h}{\tau_{ev}}, \quad \text{levitation} \sim \frac{\mu_v \lambda_v \Delta T}{\rho_v \mathcal{L}} \frac{R^2}{h^2}, \quad (4.9) $$

where the estimation of the levitation force is just based on the left-hand side of (4.1). Using the expression (3.19) for $\tau_{ev}$, one can immediately see that

$$ \frac{\text{inertia}}{\text{drag}} \sim Pr^{-1}Pe. \quad (4.10) $$

As we have $Pr \sim 1$, $Pe \ll 1$ here (cf. § 3.5), inertia can be disregarded against the drag in the present context with $\tau = \tau_{ev}$ (which does not exclude that inertia can be essential in other contexts, cf. § 7 later on). Then, it just remains to compare the drag and levitation forces. Using (4.9) on account of (3.19), one can obtain

$$ \frac{\text{drag}}{\text{levitation}} \sim \frac{\rho_v}{\rho_l} \left(\frac{h}{R}\right)^3. \quad (4.11) $$

Thus, the drag can be neglected in favour of a quasi-steady force balance like (4.1) provided that

$$ \delta^3 \ll \frac{\rho_l}{\rho_v}. \quad (4.12) $$

Given that the liquid density is much greater than the vapour density ($\rho_l/\rho_v \gg 1$), the condition (4.12) leaves quite a considerable margin for the present quasi-steady approach to be valid. It is only for sufficiently small droplets levitating too high (such that $\delta \sim (\rho_l/\rho_v)^{1/3}$) that it breaks down and the drag force should be incorporated (but still not the inertia force, according to the earlier estimations), as intuitively expected. Moreover, one can clearly see that the condition (4.12) is more restrictive in the realm of large $\delta$ than (3.20), which is further reinforced by the fact that $Pe \ll 1$. This means that, even when the drag force becomes important, the temperature and velocity fields between the droplet and the substrate can still be regarded as quasi-steady, and hence expressions like (3.8) and (3.15) are still valid. An analysis aiming at smaller $R$ (and larger $h$) and incorporating the drag force is realised in § 5 and § 6 below.

In the opposite limit of larger $R$ (and smaller $h$), the validity of the quasi-steady approach as used here is therefore not put into question. However, it is rather the full-sphericity assumption that becomes more restrictive, when (even still within $R \ll \ell_c$ and a practical sphericity of the most of the droplet) a small part of the droplet bottom gets flattened by gravity (Pomeau et al. 2012), an essential effect from the Leidenfrost viewpoint. In this regard, the result (4.4) should be understood in an intermediate asymptotic sense, as valid for $R \gg \ell_*$ but $R$ still much smaller than the bottom-flattening scale. This will be considered in more detail in § 8.
<table><caption>Table 2. Scales used to render the various quantities dimensionless in § 5. The square brackets denote the scale of the quantity inside.</caption><thead><tr><th>[r, z, s]</th><th>[S]</th><th>[v]</th><th>[p]</th><th>[F<sub>drag</sub>]</th></tr></thead><tbody><tr><td>R</td><td>R<sup>2</sup></td><td>U</td><td>μ<sub>v</sub>U/R</td><td>-μ<sub>v</sub>RU</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 7. Dimensionless drag force as a function of the relative droplet height $\delta$. Numerical results (blue open circles) agree well with the classical approximation (5.2) (red line). Black dashed lines correspond to the asymptotic behaviours.</figcaption></figure>

## 5. Basic calculations (continued): drag force

As stipulated in § 4.2, the drag force, $F_{drag}$, is required for further analysis. The present section, dedicated to $F_{drag}$, is organised as a continuation of § 3 and mirrors the same style as far as notations, non-dimensionalisation and scales are concerned. The scales relevant here are summarised in table 2 (which is the counterpart of table 1 there), where $U$ is the droplet (translation) velocity in the vertical direction (its only component here).

As the primary need for such a consideration arose in the context of large $\delta$ (cf. § 4.2), a mere use of the (dimensionless) Stokes drag $F_{drag} = 6\pi$ in an unbounded medium could be quite sufficient here (as well as in § 6), where the rigid-sphere prefactor $6\pi$ is used on account of the liquid dynamic viscosity being much larger than the vapour one. Nonetheless, for the sake of generality, we shall here proceed implying $\delta = O(1)$, all the more so as it will be particularly relevant later on in the context of § 7. Thus, the goal is to compute $F_{drag}(\delta)$.

For this purpose, we once again solve the dimensionless Stokes equations (3.9)–(3.10) with the boundary conditions (3.11)–(3.12) (although the dimensional scales are now different and given by table 2). However, the ‘evaporation’ boundary conditions (3.13) are now replaced with

$$v = e_z \qquad (5.1)$$

reflecting droplet translation along $z$. Finally, the same expression as on the right-hand side of (3.14) is used to compute $F_{drag}$.

The computation results are illustrated in figure 7 together with the following approximate expression (Guyon et al. 2012):
$$F_{drag}(\delta) = 6\pi \left(1 + \frac{1}{\delta}\right). \qquad (5.2)$$

While strictly valid in the limits $\delta \to 0$ and $\delta \to +\infty$, the result (5.2) can be seen to deviate from the numerical results by up to 7% for intermediate values of $\delta \sim 1$. A more precise fit of the numerical data is proposed in Appendix B3. Nonetheless, we shall stick in the present study to a simpler and more elegant expression (5.2), as we earlier did with (3.8) and (3.15), which will be sufficient for our purposes.

## 6. Levitation and take-off: (indispensable) dynamic approach for smaller droplets

We now turn to the case of smaller droplets (levitating higher), for which the quasi-steady approach followed through in § 4.1 breaks down. As pointed out in § 4.2, the drag force, depending on the droplet translation velocity as it soars higher while vapourising, becomes essential (although not the inertia force). Supplementing the quasi-steady force balance (4.1) with the drag force, we obtain

$$-6\pi \mu_v R \left(1 + \frac{R}{h}\right) \frac{dh}{dt} + 3\pi \frac{\mu_v \lambda_v \Delta T}{\rho_v \mathcal{L}} \frac{R^2}{h^2} \frac{R + 2h}{R + h} = \frac{4\pi}{3} \rho_l g R^3. \qquad (6.1)$$

Here, the dimensionless expression (5.2) multiplied by the scale from table 2 has been used with $U = dh/dt$. (Strictly speaking, the velocity of the centre of mass is rather given by $(dh/dt) + (dR/dt)$, but we shall disregard such nuances here.)

The natural distinguished scales are such that all terms in (6.1) and in (4.6) are of the same order of magnitude:

$$[\tilde{R}] = \left(\frac{\rho_v}{\rho_l}\right)^{2/9} \ell_*, \quad [\tilde{h}] = \left(\frac{\rho_v}{\rho_l}\right)^{-1/9} \ell_*, \quad [\tilde{t}] = \left(\frac{\rho_v}{\rho_l}\right)^{4/9} \tau_*, \qquad (6.2)$$

$$\tilde{R} = \frac{R}{[\tilde{R}]}, \quad \tilde{h} = \frac{h}{[\tilde{h}]}, \quad \tilde{t} = \frac{t}{[\tilde{t}]}, \quad \epsilon \equiv \frac{[\tilde{R}]}{[\tilde{h}]} = \left(\frac{\rho_v}{\rho_l}\right)^{1/3}, \qquad (6.3)$$

and the dynamical system becomes

$$- \left( 1 + \epsilon \frac{\tilde{R}}{\tilde{h}} \right) \frac{d\tilde{h}}{d\tilde{t}} + \frac{\tilde{R}}{\tilde{h}^2} \frac{\tilde{h} + \frac{1}{2}\epsilon\tilde{R}}{\tilde{h} + \epsilon\tilde{R}} = \frac{2}{9}\tilde{R}^2, \qquad (6.4)$$

$$\tilde{R} \frac{d\tilde{R}}{d\tilde{t}} = -1 - \frac{1}{2} \ln \left( 1 + \epsilon \frac{\tilde{R}}{\tilde{h}} \right). \qquad (6.5)$$

The distinction between the present scales (6.2) and the previously considered ones (4.3) and (4.7) is eventually owing to a small parameter given by the vapour-to-liquid density ratio $\rho_v/\rho_l \ll 1$. We note that the scales $[\tilde{R}]$ and $[\tilde{h}]$ are different here, and the typical relative height of the droplet levitation is now $\delta \sim (\rho_l/\rho_v)^{1/3} \gg 1$. It is exactly at such values of $\delta$ that criterion (4.12) of the quasi-steady force balance breaks down, as expected, which confirms the coherence of the present dynamic approach. At the same time, criterion (3.20) of the quasi-steadiness of the temperature and velocity fields between the droplet and the substrate is still satisfied, as already discussed in § 4.2.

The phase portrait of the dynamical system (6.4)–(6.5) is represented in figure 8(a) (generated using the StreamPlot command in Mathematica). Notably, while a Leidenfrost droplet is theoretically expected to ascend indefinitely within the quasi-steady approach, the presence of the drag force results in a saturation of $\tilde{h}$ towards a finite take-off value.
<figure><img src="image_10.png" /><figcaption>Figure 8. (a) Phase portrait of the dynamical system for $\tilde{R}(\tilde{t})$ and $\tilde{h}(\tilde{t})$ appropriate in the limit of a smaller droplet levitating higher (incorporating the drag force in the quasi-steady force balance) and in the scaling appropriate to that limit (tilded variables). Parameter values used: $\epsilon = 0.076$, corresponding to water droplet experiments by Celestini et al. (2012) at $\Delta T = 300$ °C (cf. Appendix A). The blue solid line shows the earlier obtained quasi-steady ‘master curve’, which coincides with the separatrix (black solid line, ‘dynamic master curve’) for larger $\tilde{R}$ but diverges ($\tilde{h} \to +\infty$) for smaller $\tilde{R}$. (b) Comparison with the mentioned experiment using the quasi-steady and dynamic master curves.</figcaption></figure>

In other words, the Leidenfrost droplets always vanish ($\tilde{R} \to 0$) at a finite height. Physically, the drag force, represented by the first term on the left-hand side of (6.1), becomes so significant towards the end of life of the droplet that it blocks the soaring tendency dictated by the levitation force (the second term ibid). In principle, the droplet can end up at whatever height, depending on the initial conditions (figure 8a). However, there is a distinguished value of the final levitation height valid for most of the droplets, at least for those starting off from a sufficiently large size. In such a case, the phase trajectory is seen to reach the separatrix fast (figure 8a), along which the droplet evolution ensues until the droplet vanishes at the distinguished final height corresponding to the separatrix

$$\tilde{h}_{fin} = 1.69(1 - 0.35 \epsilon), \qquad (6.6)$$

which was computed numerically assuming a linear dependence on the small parameter $\epsilon$. Using the scales (6.2) on account of (4.3), this can be rewritten in dimensional terms

$$h_{fin} = 1.69 \left(\frac{\rho_l}{\rho_v}\right)^{1/9} \left(\frac{\mu_v \lambda_v \Delta T}{\rho_v \rho_l g \mathcal{L}}\right)^{1/3} \left(1 - 0.35 \left(\frac{\rho_v}{\rho_l}\right)^{1/3}\right). \qquad (6.7)$$

For instance, under the conditions of the experiments by Celestini et al. (2012) with water droplets (cf. Appendix A for the parameters), we obtain $\tilde{h}_{fin} = 1.64$ and $h_{fin} = 110.36$ µm. The separatrix now becomes our new, dynamic ‘master curve’. It replaces the quasi-steady one, soaring to infinity ($\tilde{h} \to +\infty$) as $\tilde{R} \to 0$ and corresponding to

$$\frac{1}{\tilde{h}^2} \frac{\tilde{h} + \frac{1}{2}\epsilon\tilde{R}}{\tilde{h} + \epsilon\tilde{R}} = \frac{2}{9}\tilde{R} \qquad (6.8)$$

in terms of the tilded variables, which is also depicted in figure 8 (solid blue line) for comparison. The dynamic master curve is different for smaller droplets ($\tilde{R} \lesssim 3$), whereas for larger droplets the quasi-steady result is recovered, as expected.

Figure 8(b) undertakes a direct comparison with the experiments by Celestini et al. (2012) for water droplets at $\Delta T = 300$ °C. We see that the dynamic model, considered in
the present section, shows a noticeable improvement over the previously considered quasi-
steady approach. The improvement is just manifest for smaller droplets ($R \lesssim 15$ µm).

7. Dynamic approach in general

The dynamic approach considered in § 6 came as indispensable and, in this sense, forced
in the domain of smaller droplets, where the consideration of § 4 based on the quasi-
steady force balance could no longer be valid. At the same time, this made it limited to
that domain, where, in particular, only the drag force was essential while inertia could
be disregarded. In the present section, we take up a full dynamic approach, which would
presumably be valid in the entire range of droplet sizes considered in the present paper. In
§ 4, we considered just the (quasi-)equilibrium positions (heights) of the droplet. Here, we
inquire what the droplet behaviour will be if initially out of that equilibrium.

Complementing the force balance (6.1) with inertia, we arrive at

$$
\frac{4\pi}{3} \rho_l R^3 \frac{d^2h}{dt^2} = -6\pi \mu_v R \left(1 + \frac{R}{h}\right) \frac{dh}{dt} + 3\pi \frac{\mu_v \lambda_v \Delta T}{\rho_v \mathcal{L}} \frac{R^2}{h^2} \frac{R + 2h}{R + h} - \frac{4\pi}{3} \rho_l g R^3. \quad (7.1)
$$

Thus, we end up with a third-order system of ordinary differential equations given by (7.1)
and (4.6). It can in principle be (numerically) solved starting from any initial condition
R = R₀, h = h₀, dh/dt = h'₀ at t = 0 to obtain R(t) and h(t), where R₀ > 0, h₀ > 0 and h'₀
are some initial values.

The solution is illustrated in figure 9 for a water Leidenfrost droplet of an initial radius $R_0 = 30$ µm starting from various initial heights $h_0$ and velocities $h'_0$. As earlier, the parameters of the experiment by Celestini et al. (2012) are used (i.e. $\Delta T = 300$ °C and 1 atm, cf. Appendix A). We note that the chosen value of the initial radius is large enough for the quasi-steady approach to work well (cf. Figure 8b) and the corresponding quasi-steady height, obtained from (4.1) at $R = R_0$, is $h_{QS} = 50.24$ µm. The several initial heights tested are then conveniently expressed in the units of $h_{QS}$. The previously obtained quasi-steady (§ 4) and dynamic (§ 6) master curves are also shown in figure 9 for reference.

For $h_0 = h_{QS}$ with $h'_0 = 0$, we see that the solution adheres to the dynamic master curve, which coincides with the quasi-steady one for larger $R$ and is drag-force-moderated for smaller $R$, and where the incorporation of inertia into the force balance has practically no effect, as expected. For an initial height $h_0$ out of the position $h_{QS}$ (still with $h'_0 = 0$), the droplet approaches the dynamic master curve relatively fast in an oscillatory way (like a damped oscillator), and then the evolution continues along that curve. The droplets initially located too close to the substrate can rebound to considerable heights as propelled by the levitation force. The droplets starting from or propelled to considerable heights rejoin the dynamic master curve at a later time and smaller size. In this case, the rejoining already happens in a monotonic way, quite in accordance with the scenario for smaller droplets described in § 6, where inertia could be disregarded. Furthermore, the droplets finding themselves at a certain moment excessively high vapourise at some finite height before reaching the dynamic master curve, which also forms part of that scenario. Under the conditions explored in figure 9(a,b), the latter scenario occurs whenever $h_0 \lesssim (1/300)h_{QS}(R_0) = 0.18$ µm and $h_0 \gtrsim 72h_{QS}(R_0) = 3.83$ mm.

Exploring considerable initial downward velocities can be relevant in the context of
spraying. This is accomplished in figure 9(c,d) for two values of the initial height
($h_0 = 5h_{QS}$ and $h_0 = 75h_{QS}$). For $h_0 = 5h_{QS}$, finite positive values of $-h'_0$ delay
stabilisation at the master curve. Furthermore, if the initial downward velocity is
sufficiently large, the droplet rebounds so high that it eventually vapourises before reaching
the master curve. In contrast, from a much larger initial height $h_0 = 75h_{QS}$ when the
<figure><img src="image_11.png" /><figcaption>Figure 9. Leidenfrost water droplet trajectories in the plane $h$ versus $R$ computed by means of the general dynamic approach (including drag and inertia) using the parameters of Celestini et al. (2012) (cf. Appendix A). Initial radius $R_0 = 30$ µm. Various initial heights $h_0$ are tested in (a) and (b) for $h'_0 = 0$, which are here represented relative to the equilibrium quasi-steady height $h_{QS} = 50.24$ µm. Mimicking spraying, the role of an initial downward velocity $-h'_0$ comparable to or greater than the droplet fall velocity $2\rho_l g R_0^2 / (9\mu_v) \approx 0.1$ m s$^{-1}$ is explored in (c) for $h_0 = 5h_{QS}$ and in (d) for $h_0 = 75h_{QS}$. The quasi-steady and dynamic master curves (blue and red dashed lines) are also shown for reference.</figcaption></figure>

droplet vapourises before reaching the master curve with $h'_0 = 0$, a finite initial downward velocity can help propel the droplet to its stabilisation at the master curve. Once again, this is unless the velocity is excessively large so that the droplet rebounds too high.

Small rebounds around the master curve can actually reproduce certain oscillatory trends observed in the experimental points by Celestini et al. (2012), as illustrated in figure 10. However, the nature of the experimental points located too close to the substrate remains unclear, the understanding of which may require staging further experiments and thinking of physical factors not included in the present model. The present modelling indicates that there may be a certain dependence on the manner in which the Leidenfrost droplets are deposited in experiment as the initial conditions can be such that the droplet fully evaporates before reaching the master curve.

The consistency of the approach in regard of the oscillatory relaxation obtained here can be assessed as follows. Focusing just on larger droplets with $h \sim R$, the oscillation time scale $\tau_{oscil} = \sqrt{R/g}$ can be compared with the viscous and thermal time scales $\tau_{vis} \sim \rho_v R^2 / \mu_v$ and $\tau_{th} \sim R^2 / \alpha_v$ (the latter two are of the same order on account of $Pr \sim 1$ and can thus be used interchangeably in estimations). As $\tau_{vis} \tau_{th} / \tau_{oscil}^2 = Gr$, cf. (3.17), while it was estimated $Gr \ll 1$ in § 3.5, we see that $\tau_{vis}$, $\tau_{th} \ll \tau_{oscil}$. Thus, the implied
<figure><img src="image_12.png" /><figcaption>Figure 10. Relative height $\delta = h/R$ of a water Leidenfrost droplet as a function of its radius $R$. The experimental data from Celestini et al. (2012) are compared with the quasi-steady and dynamic (drag-moderated) master curves, as well as with the trajectories computed from the general dynamic approach (including both drag and inertia) using the initial conditions $h_0/h_{QS}(R_0) = 3/2$, $R_0 = \{20.8; 23.8; 28.8\} \mu m$ and $h'_0 = 0$.</figcaption></figure>

quasi-steadiness of the temperature and velocity fields does hold during the oscillation cycle, hence the sought consistency. Similarly, the ratio of the inertia and drag forces in (7.1) can be estimated at $\sim (\rho_l/\rho_v)Gr^{1/2}$ taking $\tau_{oscil}$ as the time scale. Even if $Gr \ll 1$, this can be superseded by $(\rho_l/\rho_v) \gg 1$ for larger droplets so that inertia dominates, hence the observed oscillatory relaxation. For smaller droplets, however, as $Gr$ decreases drastically with $R$, it is the drag that comes to dominate, hence a monotonic relaxation and the scenario of § 6.

## 8. Link to a more global picture

The larger the droplet is, the narrower the gap between the spherical droplet and the substrate becomes, as put into evidence by (4.4). As mentioned at the end of § 4.2, a key limitation to the present analysis from the side of larger droplets is expected to be caused by a deviation from sphericity within such a narrow gap, even if the droplet as a whole might still remain largely spherical. As this region is crucial for vapour generation and heat transfer, in spite of its smallness, any morphological change therein can have a significant impact on the Leidenfrost phenomenon. On the other hand, it is through such a morphological change at the bottom of the droplet that a transition from the small spherical to larger non-spherical Leidenfrost droplets is bridged, which is touched upon in the present section.

As established by Pomeau et al. (2012), cf. their (26), a significant deviation from sphericity at the bottom of the droplet in the aforementioned sense occurs for $R$ starting from $R \sim \ell_i$, where the length scale $\ell_i$ is given by

$$ \ell_i = (\ell_*^3 \ell_c^4)^{1/7} \qquad (8.1) $$

(all rewritten in our present notations). For $R \ll \ell_i$, the droplet can be considered fully spherical (even at the bottom) and the analysis of the present paper holds. As $\ell_* \ll \ell_c$, (8.1) implies that $\ell_* \ll \ell_i \ll \ell_c$. For our reference case of a water droplet on a superheated substrate with $\Delta T = 300$ °C (Celestini et al. 2012, cf. Appendix A), we obtain $\ell_i = 367 \mu m$.
<figure><img src="image_13.png" /></figure>

Figure 11. Vapour film thickness *h* under a water Leidenfrost droplet as a function of its radius *R*. Theoretical predictions reproduced from Sobac *et al.* (2014, 2021) for large deformed droplets ($R \gtrsim l_i$), and by the present model for small spherical droplets ($R \ll l_i$) are compared with experimental data from Burton *et al.* (2012) and Celestini *et al.* (2012), respectively. Unlike the previous examples, the computations are here done with $\Delta T = 270$ °C to follow Burton *et al.* (2012) (cf. also Appendix A), whereas the data by Celestini *et al.* (2012) still correspond to $\Delta T = 300$ °C (hence a slight misplacement relative to the theoretical curve as compared with previous figures). For large droplets, the *h* curve splits into two branches $h_{\text{neck}}$ and $h_{\text{centre}}$ at the point where the vapour layer between the droplet and the substrate adopts a ‘pocket-like’ structure edged by a narrow annular neck, such that the minimum thickness no longer corresponds to the centre and switches to the neck.

Figure 11 illustrates the Leidenfrost effect at large, over four decades of the droplet size. It combines the present results for the (small) spherical Leidenfrost droplets to the left of the figure (the dynamic master curve) with the results for the usual (large and deformed) droplets reproduced from Sobac *et al.* (2014, 2021) to the right. We note that for the non-spherical droplets the radius *R* is here defined as the radius of the vertical projection on the substrate (i.e. as the maximum horizontal radius). Two sets of experimental results, corresponding to the two drastically different size domains, are plotted alongside the theoretical curves: the ones by Celestini *et al.* (2012) and by Burton *et al.* (2012). The overlapping occurs at *R* ~ *l*ᵢ, quite as expected. Nevertheless, it does not appear to happen smoothly, which is most definitely due to some accuracy loss towards the limit of small near-spherical droplets within the model employed by Sobac *et al.* (2014, 2021). Notably, it is in this intermediate (overlapping) region that an absolute minimum of the vapour layer thicknesses is attained: from there, *h* increases both towards the small spherical droplets as a take-off precursor and towards the large deformed droplets.

## 9. Conclusions

The dynamics of small spherical Leidenfrost droplets has been investigated theoretically, yielding valuable new insights into their final stages of existence.

After numerically calculating the fluxes, evaporation rate and forces for a spherical Leidenfrost droplet interacting with a superheated flat substrate (under the verifiable assumptions of quasi-stationarity and low Reynolds and Péclet numbers) and coming up with simple formulae approximating these data as a function of the reduced height $\delta = h/R$ (all respecting the asymptotic behaviours as $\delta \to 0$ and $\delta \to +\infty$, also studied
here), a theoretical model has been developed which allows an accurate prediction of the droplet height *h* as a function of the physical parameters without any fitting parameters.

First, by balancing the droplet’s weight against the upward hydrodynamic force induced by evaporation, a ‘quasi-steady master curve’ relating droplet height *h* to droplet radius *R* was derived. This curve follow the *h* ∝ *R*⁻¹/² scaling law formulated by Celestini *et al.* (2012) and Pomeau *et al.* (2012). However, the prefactor is not constant and rather varies from 3/2 when *R* → +∞ to 3/√2 as *R* → 0.

Furthermore, our analysis reveals that the aforementioned classical quasi-steady description, while capturing the general trend of the ‘take-off’ phenomenon, is unable to accurately reproduce the experimental data of Celestini *et al.* (2012). Dynamic effects, especially those related to frictional forces, are crucial to accurately describe the take-off at small scales. Therefore, a ‘dynamic master curve’, drag-moderated, has been also derived and well reproduces the experimental data of Celestini *et al.* (2012). As a consequence of the friction effect, the Leidenfrost droplets disappear at a finite height, the value of which turns out to be universal for sufficiently large initial droplets. This is in contrast to the prediction of the quasi-stationary model, which suggests an infinite final height. A formula that predicts this universal final height has been established.

In addition, a general dynamic model including both drag and inertia effects has been used to investigate the influence of initial conditions on droplet dynamics. For an initial height $h_0$ out of the equilibrium position $h_{QS}$, a sufficiently large droplet approaches the ‘dynamic master curve’ relatively quickly and in an oscillatory manner (like a damped oscillator), and then continues to evolve along this curve. Such small rebounds can indeed reproduce certain oscillatory trends observed in the experimental points of Celestini *et al.* (2012). This scenario holds for initial conditions not too far from equilibrium; otherwise, the droplets that find themselves too high at a given moment evaporate at a finite height before reaching the dynamic master curve. However, even for droplets initially situated that high, an appreciable initial downward velocity (e.g. due to spraying) can help them reach the master curve. All this highlights a certain dependence of the Leidenfrost droplet deposition manner on the result.

Combining the present model (valid when $R < \ell_i$) with the one of Sobac *et al.* (2014, 2021) for larger deformed droplets (valid when $R > \ell_i$), we offer a comprehensive picture of the shape and elevation of Leidenfrost droplets across a the full range of stable axisymmetric shapes, spanning four decades of droplet sizes. These studies also align with the hierarchy of length scales $\ell_* < \ell_i < \ell_c$ pointed out by Pomeau *et al.* (2012), emphasising the dominant physical mechanisms and associated scalings.

We hope this research will stimulate further theoretical investigations, particularly in the intermediate region ($R \approx \ell_i$) where the vapour layer is thinnest, and encourage further experimental studies on Leidenfrost droplets with $R \lesssim \ell_i$, an area that remains underexplored.

**Funding.** B.S. gratefully acknowledges the support from Centre National de la Recherche Scientifique – CNRS, A.R. from BELSPO and ESA PRODEX Evaporation, and P.C. from the Fonds de la Recherche Scientifique – FNRS.

**Declaration of interests.** The authors report no conflict of interest.

## Appendix A. Properties

The properties used in this work are provided in table 3 (Sobac *et al.* 2015, 2017).

To facilitate continuity with some previous studies of the usual (larger non-spherical) Leidenfrost droplets (Sobac *et al.* 2014, 2021; Chantelot & Lohse 2021), apart from the
<table><caption>Table 3. Parameter values at 1 atm (where relevant). Liquid density $ρ_l$, latent heat $L$, surface tension $γ$ and capillary length $ℓ_c = \sqrt{\gamma/ρ_l g}$ ($g$ gravitational acceleration) at the boiling temperature $T_{sat}$; superheat $ΔT = T_w - T_{sat}$, with $T_w$ the substrate temperature; vapour thermal conductivity $λ_v$, viscosity $μ_v$, density $ρ_v$ at the mid-temperature $1/2(T_{sat} + T_w)$; take-off scale $ℓ_*$, cf. (4.3); non-sphericity scale $ℓ_i$, cf. (8.1); $ε = (ρ_v/ρ_l)^{1/3}$; evaporation number $E$, cf. (A1). The first and third cases (rows) correspond to Celestini et al. (2012), while the second corresponds to Burton et al. (2012) and the fourth corresponds to Lyu et al. (2019).</caption><thead><tr><th></th><th><i>T</i><sub>sat</sub><br>(°C)</th><th><i>ρ</i><sub>l</sub><br>(kg/m<sup>3</sup>)</th><th><i>L</i><br>(kJ/kg)</th><th><i>γ</i><br>(mN/m)</th><th><i>ℓ</i><sub>c</sub><br>(mm)</th><th>Δ<i>T</i><br>(°C)</th><th><i>λ</i><sub>v</sub><br>(mW/mK)</th><th><i>μ</i><sub>v</sub><br>(μPa s)</th><th><i>ρ</i><sub>v</sub><br>(kg/m<sup>3</sup>)</th><th><i>ℓ</i><sub>*</sub><br>(μm)</th><th><i>ℓ</i><sub>i</sub><br>(μm)</th><th><i>ε</i><br>(-)</th><th><i>E</i> × 10<sup>6</sup><br>(-)</th></tr></thead><tbody><tr><td>Liquid</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Water</td><td>100</td><td>960</td><td>2555</td><td>59</td><td>2.50</td><td>300</td><td>36.9</td><td>18.6</td><td>0.42</td><td>28.5</td><td>367</td><td>0.076</td><td>1.47</td></tr><tr><td>Water</td><td>100</td><td>960</td><td>2555</td><td>59</td><td>2.50</td><td>270</td><td>35.8</td><td>18.0</td><td>0.43</td><td>26.7</td><td>357</td><td>0.077</td><td>1.21</td></tr><tr><td>Ethanol</td><td>79</td><td>728</td><td>840</td><td>17</td><td>1.56</td><td>321</td><td>25.4</td><td>15.6</td><td>1.10</td><td>26.9</td><td>274</td><td>0.115</td><td>5.12</td></tr><tr><td>Ethanol</td><td>79</td><td>728</td><td>840</td><td>17</td><td>1.56</td><td>217</td><td>22.7</td><td>14.0</td><td>1.22</td><td>21.1</td><td>247</td><td>0.119</td><td>2.48</td></tr></tbody></table>

parameters defined in the main text, we here also follow the (dimensionless) evaporation number

$$ \mathcal{E} = \frac{\mu_v \lambda_v \Delta T}{\gamma \mathcal{L} \rho_v \ell_c} \qquad (A1) $$

We note the formula $ℓ_* = \mathcal{E}^{1/3}ℓ_c$ and $ℓ_i = \mathcal{E}^{1/7}ℓ_c$ resulting from (4.3), (8.1) and (A1).

## Appendix B. Some more precise fits of numerical data

In addition to the simplified fits (approximations) outlined in the main text, more precise fits are proposed herein. These fits are presented alongside and compared with the numerically computed data in figure 12. The asymptotic behaviours derived in Appendix C are also depicted.

### B.1. Fit for $J(\delta)$

The expression (3.8) can further be improved as follows:

$$ J(\delta) = 4\pi \left[ 1 + \frac{1}{2} \ln \left( 1 + \frac{1}{\delta} \right) - \left( 1 - \frac{1}{2} \ln 2 - \gamma \right) \frac{1}{1 + 50.8\delta^2} \right], \quad (B1) $$

where $γ = 0.577216$ is the Euler constant, while the coefficient 50.8 results from fitting to the numerically computed data. The expression (B1) respects the two-term asymptotic expansions both as $δ → 0$ and as $δ → +∞$ obtained in Appendices C1 and C2, respectively.

### B.2. Fit for $F_{ev}(\delta)$

A more precise expression than (3.15) is given by

$$ F_{ev}(\delta) = \frac{3\pi}{\delta^2} \left( 1 + \frac{\delta}{0.924 + \delta} \right), \quad (B2) $$

where the coefficient 0.924 is obtained by fitting to the numerically computed data. As earlier (3.15), the expression (B2) respects the leading-order asymptotic behaviour in the limits $δ → 0$ and $δ → +∞$ considered in Appendices C3 and C4, respectively.
<figure><img src="image_14.png" /><figcaption>Figure 12. (a) Evaporation rate of and forces acting on a spherical Leidenfrost droplet as a function of its reduced height $\delta$. Numerically computed data are presented alongside their asymptotic behaviours and compared with fits of two different levels of precision. (b) Corresponding ratio of the numerical data to the proposed fits as a function of $\delta$.</figcaption></figure>

### B.3. Fit for $F_{drag}(\delta)$

A better approximation than (5.2) is provided by the formula

$$F_{drag}(\delta) = 6\pi \left( 1 + \frac{1}{\delta} + 1.161 \frac{1 + 26.01\delta}{1 + 62.447\delta + 187.12\delta^2 + 2.514\delta^3} \right), \quad (B3)$$

where the coefficients in the last term result from fitting the expression to the numerically computed data. The leading-order asymptotic behaviour as $\delta \to 0$ and the two-term asymptotic expansion as $\delta \to +\infty$ are for (B3) the same as for (5.2).
Appendix C. Asymptotic behaviours

C.1. Case of $J$ as $\delta \to 0$

We now consider the case when the sphere (dimensionless radius unity) is close to the substrate ($\delta \ll 1$). Its profile $z = h(r)$ ($z = 0$ at the substrate) in a small vicinity of the downmost point is well approximated by a parabola $h = \delta + (1/2)r^2$ matching with the outer circular droplet shape. The evaporation flux is given by heat conduction across this thin vapour layer, viz. $j = 1/h$ in dimensionless form. Integrating over such a small vicinity up to a reference point $r = r_1 \ll 1$, one obtains $J_1 = 2\pi \int_0^{r_1} r j dr = 2\pi \ln(\delta + (1/2)r_1^2) - 2\pi \ln \delta$ for the first contribution into the evaporation rate $J$. One can observe a logarithmic divergence and assume the sought asymptotic behaviour in the form $J \sim 2\pi(-\ln \delta + \text{const})$ as $\delta \to 0$. However, to fully determine const, one needs to consider the contribution $J_2$ from the rest of the sphere, such that eventually $J = J_1 + J_2$. For this purpose, it suffices to solve the problem (3.2)–(3.5) for a unit sphere lying on the substrate (with formally $\delta = 0$). This can be done numerically with $J_2$ determined by integrating (3.6) over the rest of the sphere up to the point $r = r_1$ on the lower surface. The result diverges in the limit $r_1 \to 0$: as $J_2 = -4\pi \ln r_1 + \text{const}_1$, where $\text{const}_1$ is determined numerically by choosing sufficiently small $r_1$. In this way, we finally arrive at $J \sim 2\pi(-\ln \delta + 1.85)$ as $\delta \to 0$ (the terms with $r_1$ cancelling out in $J_1 + J_2$). Within the computation precision and with $\gamma$ being the Euler constant, this can be rewritten as $J \sim 2\pi(-\ln \delta + \ln 2 + 2\gamma)$ as $\delta \to 0$, which is taken into account when constructing the fit (B1). This (exact) value of the constant can in principle be obtained from the exact solution of the present heat conduction problem in curvilinear coordinates (similarly to e.g. Lebedev 1972; Popov 2005), although we here limit ourselves to corresponding numerical solutions. The approximation (3.8) respects exactly the logarithmic divergence, but just approximately the constant.

C.2. Case of $J$ as $\delta \to +\infty$

When the sphere is far away from the substrate, the leading-order result is as for the sphere in an unbounded medium

$$T = 1 - 1/r_{sph}, \qquad (C1)$$

where $r_{sph} = \sqrt{r^2 + (z - \delta - 1)^2}$ is the spherical radial coordinate from the centre of the sphere (the difference with the cylindrical radial coordinate $r$ is marked by the subscript). Therefore, $j = \partial_{r_{sph}} T|_{r_{sph}=1} = 1$ and hence $J = 4\pi$. The first correction comes from the sphere reflection in the substrate. The image sphere adds the following primary contribution into the temperature field in the original domain (above the substrate, $z > 0$)

$$T_{im} = 1/r_{sph,im}, \qquad (C2)$$

where $r_{sph,im} = \sqrt{r^2 + (z + \delta + 1)^2}$ and the (additional) subscript ‘im’ is associated with the image. The primary effect of (C2) is to increase the local ambient temperature around the original sphere from $T = 1$ to $T = 1 + 1/(2\delta)$. As $T = 0$ at the sphere surface, cf. (3.5), the values of $j$ and $J$ are then increased in the same proportion. Thus, we arrive at $J \sim 4\pi(1 + 1/(2\delta))$ as $\delta \to +\infty$, which is respected by both the fit (3.8) and (B1).

C.3. Case of $F_{ev}$ as $\delta \to 0$

In dimensionless terms (scales provided in table 1), the lubrication equation in the thin vapour layer between the substrate and the sphere can be written as $(1/12)r^{-1}\partial_r(rh^3\partial_r P_v) + 1/h = 0$ (cf. also Sobac et al. 2021), where $P_v$ is the
vapour pressure excess over the ambient one (hence $P_v \to 0$ far away). Integrating with $h = \delta + (1/2)r^2$ and on account of symmetry $\partial_r P_v|_{r=0} = 0$, one obtains $\partial_r P_v = -12h^{-3}r^{-1}\ln(h/\delta)$. The leading-order contribution into $F_{ev}$ is given by $F_{ev} = 2\pi \int_0^{+\infty} r P_v dr = \pi r^2 P_v|_{r \to +\infty} - \pi \int_0^{+\infty} r^2 \partial_r P_v dr = 3\pi/\delta^2$, which is the sought behaviour as $\delta \to 0$ and is respected in (B2).

### C.4. Case of $F_{ev}$ as $\delta \to +\infty$

Assuming $\delta \gg 1$, we shall distinguish three contributions into the sought asymptotic behaviour: $F_{ev} = F_{ev1} + F_{ev2} + F_{ev3}$, which are all of the same order.

First, the leading-order (spherically symmetric) flow field

$$ \mathbf{v} = \frac{\mathbf{r}_{sph}}{r_{sph}^3} \qquad (C3) $$

from our evaporating sphere ($\mathbf{r}_{sph}$ being the position vector from the sphere's centre) is supplemented by the one from the image sphere

$$ \mathbf{v}_{im} = \frac{\mathbf{r}_{sph,im}}{r_{sph,im}^3} \qquad (C4) $$

(in the original domain $z > 0$). At the location of the original sphere ($r = 0, z = \delta + 1$), in the limit $\delta \gg 1$, the velocity field (C4) gives rise to a (quasi-)uniform streaming velocity $v_0 = 1/(4(\delta + 1)^2) \sim 1/(4\delta^2)$ directed vertically upwards. This in turn gives rise to the Stokes drag $6\pi \mu_v R v_0$ upon the original sphere (in dimensional terms). In our present dimensionless terms (cf. table 1), this amounts to $F_{ev1} = 6\pi v_0 = 3\pi/(2\delta^2)$.

Second, the superposition of the velocity fields (C3) and (C4) does satisfy the impermeability condition at the substrate: $v_z = 0$ at $z = 0$ (hereafter $v_r$ and $v_z$ are the $r$- and $z$-components of the velocity field). However, the no-slip condition $v_r = 0$ at $z = 0$ is not satisfied. To remedy this, we consider another contribution into the velocity field, the addition of which permits to observe the no-slip condition. We proceed in terms of the streamfunction $\psi$

$$ v_r = \frac{1}{r} \partial_z \psi, \quad v_z = -\frac{1}{r} \partial_r \psi. \qquad (C5) $$

The Stokes equation can be written as (cf. e.g. Happel & Brenner 1965)

$$ E^2 E^2 \psi = 0, \quad E^2 = \partial_{rr} - \frac{1}{r} \partial_r + \partial_{zz}, \qquad (C6) $$

which is solved in the domain $z > 0$ with the boundary conditions

$$ \psi = 0, \quad \partial_z \psi = -\frac{2r^2}{[(\delta + 1)^2 + r^2]^{3/2}} \quad \text{at } z = 0, \quad \psi/r_{sph}^2 \to 0 \text{ at infinity.} \qquad (C7) $$

The slip velocity in the second condition (C7) is directed towards the axis and is such as to offset the corresponding contribution from the sum of (C3) and (C4). One can verify that

$$ \psi = -\frac{2r^2 z}{[(\delta + 1 + z)^2 + r^2]^{3/2}} \qquad (C8) $$

is an exact solution of the problem (C6), (C7). Eventually, we are just interested in the velocity field value $v_0$ at the location of the original sphere ($r = 0, z = \delta + 1$). Using (C5) and (C8), one obtains $v_0 = 1/(2(\delta + 1)^2) \sim 1/(2\delta^2)$ directed vertically upwards. As with the first contribution, the Stokes drag considerations lead to the result $F_{ev2} = 6\pi v_0 = 3\pi/\delta^2$.
Third, the temperature field (C2) from the image sphere gives rise not only to an effective uniform temperature increase in the original sphere surrounding, already taken into account in Appendix C2, but also to a (dimensionless) temperature gradient $\partial_z T \sim -1/(4\delta^2)$. This breaks down the spherical symmetry of the evaporation flux (j no longer constant along the sphere surface) and of the evaporative flow (a correction upon (C3)). Hydrodynamically, this can engender an additional force contribution $F_{ev3}$. We proceed with the analysis using the spherical coordinates $\{r_{sph}, \theta\}$ related to the original sphere, such that $r = r_{sph} \sin \theta$ and $z = \delta + 1 + r_{sph} \cos \theta$. Then the problem for the mentioned (gradient-related) part of the temperature field around the original sphere can be formulated as

$$ \nabla^2 T = 0, \quad \nabla^2 = \partial_{r_{sph}r_{sph}} + \frac{2}{r_{sph}} \partial_{r_{sph}} + \frac{1}{r_{sph}^2 \sin \theta} \partial_{\theta} \sin \theta \partial_{\theta}, \quad (C9) $$

$$ T = 0 \text{ at } r_{sph} = 1, \quad T \sim -\frac{1}{4\delta^2} r_{sph} \cos \theta \text{ as } r_{sph} \to +\infty, \quad (C10) $$

to be solved in the domain $r_{sph} > 1$. The infinity in (C10) formally corresponds to $1 \ll r_{sph} \ll \delta$. The solution of (C9) with (C10) is

$$ T = -\frac{1}{4\delta^2} \left( r_{sph} - \frac{1}{r_{sph}^2} \right) \cos \theta. \quad (C11) $$

Therefore,

$$ j = \partial_{r_{sph}} T |_{r_{sph}=1} = -\frac{3}{4\delta^2} \cos \theta, \quad (C12) $$

which shows that the present contribution corresponds to evaporation reduction at the upper part of the sphere ($j < 0$ for $0 \le \theta < \pi/2$) and intensification at the lower part of the sphere ($j > 0$ for $\pi/2 < \theta \le \pi$), closer to the substrate, as expected. To calculate the flow induced by (C12), we work once again in terms of the streamfunction, now in the spherical coordinates:

$$ v_r = \frac{1}{r_{sph}^2 \sin \theta} \partial_{\theta} \psi, \quad v_{\theta} = -\frac{1}{r_{sph} \sin \theta} \partial_{r_{sph}} \psi \quad (C13) $$

for the $r_{sph}$- and $\theta$-components of the velocity field. The problem is formulated in the domain $r_{sph} > 1$. The Stokes equation (cf. e.g. Happel & Brenner 1965) and the boundary conditions can be written as

$$ E^2 E^2 \psi = 0, \quad E^2 = \partial_{r_{sph}r_{sph}} + \frac{\sin \theta}{r_{sph}^2} \partial_{\theta} \frac{1}{\sin \theta} \partial_{\theta}, \quad (C14) $$

$$ \psi = -\frac{3}{8\delta^2} \sin^2 \theta, \quad \partial_{r_{sph}} \psi = 0 \text{ at } r_{sph} = 1, \quad \psi/r_{sph}^2 \to 0 \text{ as } r_{sph} \to +\infty. \quad (C15) $$

The first condition (C15) corresponds to $v_r = j$ (in our dimensionless terms) on account of (C12) and (C13), while the second condition (C15) to no slip ($v_{\theta} = 0$). The solution of (C14) with (C15) is

$$ \psi = -\frac{3}{8\delta^2} \left( r_{sph} + \frac{1}{r_{sph}} \right) \frac{\sin^2 \theta}{2}. \quad (C16) $$

The force acting on the sphere in the z-direction is determined by ($-4\pi$) times the ‘Stokeslet’ prefactor, i.e. the one at the term $r_{sph} \sin^2 \theta/2$ (cf. Happel & Brenner 1965). Thus, from (C16), one obtains $F_{ev3} = 3\pi/(2\delta^2)$.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Relative height $\delta = h/R$ of ethanol droplets as a function of $R$. The dynamic and quasi-steady master curves are here compared with the experimental data by Celestini *et al.* (2012) in (a) and Lyu *et al.* (2019) in (b), cf. table 3 for the parameter values.</figcaption></figure>

Summing up the three contributions, one finally obtains $F_{ev} = 6\pi/\delta^2$ (as $\delta \to +\infty$). It is noteworthy that the asymptotic behaviour remains $O(\delta^{-2})$ in both limits: as $\delta \to +\infty$ and as $\delta \to 0$ (cf. Appendix C3). Yet, the prefactors are twice different.

## Appendix D. Case of ethanol

Since Celestini *et al.* (2012) and Lyu *et al.* (2019) have conducted experiments with ethanol droplets, it is worthwhile to compare the present theoretical predictions with those experimental results. Although the predictions match the form and order of magnitude of the experimental results, the agreement is less satisfactory than in the case of water. Specifically, we overestimate the data from Celestini *et al.* (2012) and underestimate the data from Lyu *et al.* (2019), cf. figure 13. Such multidirectional discrepancies make it challenging to propose a hypothesis that could explain the differences, which might be attributed to experimental errors or some physical ingredient missing in the model that is important for ethanol, but not for water. Further investigation (including experimental) is needed to address this issue.

## REFERENCES

* AJAEV, V.S. & KABOV, O.A. 2021 Levitation and self-organization of droplets. *Annu. Rev. Fluid Mech.* **53** (1), 203–225.

* BAIER, T., DUPEUX, G., HERBERT, S., HARDT, S. & QUÉRÉ, D. 2013 Propulsion mechanisms for Leidenfrost solids on ratchets. *Phys. Rev. E* **87** (2), 021001.

* BERGEN, J.E., BASSO, B.C. & BOSTWICK, J.B. 2019 Leidenfrost drop dynamics: exciting dormant modes. *Phys. Rev. Fluids* **4** (8), 083603.

* BIANCE, A.-L., CLANET, C. & QUÉRÉ, D. 2003 Leidenfrost drops. *Phys. Fluids* **15** (6), 1632–1637.

* BOERHAAVE, H. 1732 *Elementae Chemiae*, vol. 1. Lugduni Batavorum.

* BOUILLANT, A., COHEN, C., CLANET, C. & QUÉRÉ, D. 2021a Self-excitation of Leidenfrost drops and consequences on their stability. *Proc. Natl Acad. Sci. USA* **118** (26), e2021691118.

* BOUILLANT, A., LAFOUX, B., CLANET, C. & QUÉRÉ, D. 2021b Thermophobic Leidenfrost. *Soft Matter* **17** (39), 8805–8809.

* BOUILLANT, A., MOUTERDE, T., BOURRIANNE, P., LAGARDE, A., CLANET, C. & QUÉRÉ, D. 2018 Leidenfrost wheels. *Nat. Phys.* **14** (12), 1188–1192.

* BRUNET, P. & SNOEIJER, J.H. 2011 Star-drops formed by periodic excitation and on an air cushion – a short review. *Eur. Phys. J. Special Topics* **192** (1), 207–226.

* BRUTIN, D. (ed.) 2015 *Droplet Wetting and Evaporation: From Pure to Complex Fluids*. Academic Press.
BURTON, J.C., SHARPE, A.L., VAN DER VEEN, R.C.A., FRANCO, A. & NAGEL, S.R. 2012 Geometry of the vapor layer under a Leidenfrost drop. *Phys. Rev. Lett.* **109** (7), 074301.

CELESTINI, F., FRISCH, T. & POMEAU, Y. 2012 Take off of small Leidenfrost droplets. *Phys. Rev. Lett.* **109** (3), 034501.

CHAKRABORTY, I., CHUBYNSKY, M.V. & SPRITTLES, J.E. 2022 Computational modelling of Leidenfrost drops. *J. Fluid Mech.* **936**, A12.

CHANTELOT, P. & LOHSE, D. 2021 Drop impact on superheated surfaces: short-time dynamics and transition to contact. *J. Fluid Mech.* **928**, A36.

DODD, L.E., AGRAWAL, P., GERALDI, N.R., XU, B.B., WELLS, G.G., MARTIN, J., NEWTON, M.I., MCHALE, G. & WOOD, D. 2020 Planar selective Leidenfrost propulsion without physically structured substrates or walls. *Appl. Phys. Lett.* **117** (8), 081601.

DODD, L.E., AGRAWAL, P., PARNELL, M.T., GERALDI, N.R., XU, B.B., WELLS, G.G., STUART-COLE, S., NEWTON, M.I., MCHALE, G. & WOOD, D. 2019 Low-friction self-centering droplet propulsion and transport using a Leidenfrost herringbone-ratchet structure. *Phys. Rev. Appl.* **11** (3), 034063.

DUPEUX, G., LE MERRER, M., CLANET, C. & QUÉRÉ, D. 2011 Trapping Leidenfrost drops with crenelations. *Phys. Rev. Lett.* **107** (11), 114503.

GRAEBER, G., REGULAGADDA, K., HODEL, P., KÜTTEL, C., LANDOLF, D., SCHUTZIUS, T.M. & POULIKAKOS, D. 2021 Leidenfrost droplet trampolining. *Nat. Commun.* **12** (1), 1727.

GUYON, E., HULIN, J.-P. & PETIT, L. 2012 *Hydrodynamique physique 3e édition*. EDP Sciences.

HAPPEL, J. & BRENNER, H. 1965 *Low Reynolds Number Hydrodynamics*. Prentice-Hall.

LEBEDEV, N.N. 1972 *Special Functions and Their Applications*. Revised edition. Dover Publications.

LEIDENFROST, J.G. 1756 *De Aquae Communis Nonnullis Qualitatibus Tractatus*. Duisburg: Ovenius.

LI, A., LI, H., LYU, S., ZHAO, Z., XUE, L., LI, Z., LI, K., LI, M., SUN, C. & SONG, Y. 2023 Tailoring vapor film beneath a Leidenfrost drop. *Nat. Commun.* **14** (1), 2646.

LINKE, H., ALEMÁN, B.J., MELLING, L.D., TAORMINA, M.J., FRANCIS, M.J., DOW-HYGELUND, C.C., NARAYANAN, V., TAYLOR, R.P. & STOUT, A. 2006 Self-propelled Leidenfrost droplets. *Phys. Rev. Lett.* **96** (15), 154502.

LIU, D. & TRAN, T. 2020 Size-dependent spontaneous oscillations of Leidenfrost droplets. *J. Fluid Mech.* **902**, A21.

LYU, S., MATHAI, V., WANG, Y., SOBAC, B., COLINET, P., LOHSE, D. & SUN, C. 2019 Final fate of a Leidenfrost droplet: explosion or takeoff. *Sci. Adv.* **5** (5), eaav8081.

MA, X. & BURTON, J.C. 2018 Self-organized oscillations of Leidenfrost drops. *J. Fluid Mech.* **846**, 263–291.

MA, X., LIÉTOR-SANTOS, J.-J. & BURTON, J.C. 2017 Star-shaped oscillations of Leidenfrost drops. *Phys. Rev. Fluids* **2** (3), 031602.

MARENGO, M. & DE CONINCK, J. (eds) 2022 *The Surface Wettability Effect on Phase Change*. Springer.

MARÍN, Á.G., ARNALDO DEL CERRO, D., RÖMER, G.R.B.E., PATHIRAJ, B., HUIS IN 'T VELD, A. & LOHSE, D. 2012 Capillary droplets on Leidenfrost micro-ratchets. *Phys. Fluids* **24** (12), 122001.

POMEAU, Y., LE BERRE, M., CELESTINI, F. & FRISCH, T. 2012 The Leidenfrost effect: from quasi-spherical droplets to puddles. *Comptes Rendus Mécanique* **340** (11-12), 867–881.

POPOV, Y.O. 2005 Evaporative deposition patterns: spatial dimensions of the deposit. *Phys. Rev. E* **71** (3), 036313.

QUÉRÉ, D. 2013 Leidenfrost dynamics. *Annu. Rev. Fluid Mech.* **45** (1), 197–215.

SNOEIJER, J.H., BRUNET, P. & EGGERS, J. 2009 Maximum size of drops levitated by an air cushion. *Phys. Rev. E* **79** (3), 036307.

SOBAC, B., REDNIKOV, A., DORBOLO, S. & COLINET, P. 2014 Leidenfrost effect: accurate drop shape modeling and refined scaling laws. *Phys. Rev. E* **90** (5), 053011.

SOBAC, B., REDNIKOV, A., DORBOLO, S. & COLINET, P. 2017 Self-propelled Leidenfrost drops on a thermal gradient: a theoretical study. *Phys. Fluids* **29** (8), 082101.

SOBAC, B., REDNIKOV, A., DORBOLO, S. & COLINET, P. 2021 Erratum: Leidenfrost effect: accurate drop shape modeling and refined scaling laws [Phys. Rev. E 90, 053011 (2014)]. *Phys. Rev. E* **103** (3), 039901(E).

SOBAC, B., TALBOT, P., HAUT, B., REDNIKOV, A. & COLINET, P. 2015 A comprehensive analysis of the evaporation of a liquid spherical drop. *J. Colloid Interface Sci.* **438**, 306–317.

SOTO, D., LAGUBEAU, G., CLANET, C. & QUÉRÉ, D. 2016 Surfing on a herringbone. *Phys. Rev. Fluids* **1** (1), 013902.

STEWART, S.M. 2022 Leidenfrost drop dynamics: a forgotten past and modern day rediscoveries. *Eur. J. Phys.* **43** (2), 023001.