2.1 Physical model

We consider the levitation of a stationary rigid sphere (radius $a^*$ , density $\rho_s^*$ ) above an infinite liquid bath (density ${\rho}_b^*$ ). The levitation is sustained by a layer formed of the vapour of either the liquid bath or the sphere (thermal conductivity $k^*$ , density $\rho^*$ , viscosity $\mu^*$ ). Assuming axial symmetry, we denote by $h^*(r^*)$ the height of the liquid bath measured relative to the bottom of the sphere, with $r^*$ the radial distance from the symmetry axis; we denote by $h^*_{\infty}$ the unperturbed height of the liquid bath and refer to it as the sphere displacement (Figure 1(a)). The liquid–vapour interfacial tension is $\gamma^*$ and the gravitational acceleration is $g^*$ . We use a superscript asterisk to indicate a dimensional quantity.

Figure 1.

(a) Dimensional schematic of the physical problem formulated in Section 2.1. (b) Dimensionless schematic corresponding to the ‘unified small-slope’ model developed in Section 2.3.

Vapour production in this system is associated with a temperature difference (absolute magnitude $\Delta T^*$ ) between the sphere and bath, both of which assumed isothermal on account of their thermal conductivities being sufficiently high. Specifically, we assume that phase change occurs at the boundary of the cooler object, which is at its temperature for either evaporation or sublimation. This translates into the following condition on the vapour velocity at the phase-changing boundary [Reference Pomeau, Le Berre, Celestini and Frisch26]:

(2.1) \begin{equation}\mathbf{v}^* = \frac{k^*}{\rho^*l^*}(\hat{\mathbf{n}}\boldsymbol{\cdot}\boldsymbol{\nabla}^* T^*)\hat{\mathbf{n}},\end{equation}

where $l^*$ is the pertinent latent heat, $\hat{\mathbf{n}}$ is a normal unit vector pointing into the vapour layer and $\boldsymbol{\nabla}^* T^*$ is the local vapour temperature gradient. We assume that heat is transported solely by conduction through the vapour layer. We additionally assume that all of the physical quantities mentioned above can be considered constant on the time scale of interest and that the dynamical effect of the liquid bath flow can be neglected on the assumption that the liquid is much more viscous than the vapour.

The problem thus consists of (i) the Stokes equations governing the flow in the vapour layer (supposed to fill the entire fluid domain), subject to (2.1) on the phase-changing surface and zero velocity on the other; (ii) the dynamical boundary condition at the liquid–vapour interface, in which the dynamical effect of the liquid bath enters through capillarity and the hydrostatic pressure

(2.2) \begin{equation}-\rho_b^*g^*\left(h^*(r^*)-h^*_{\infty}\right),\end{equation}

which is taken to vanish at the unperturbed bath height; (iii) the condition that the gravity force on the sphere is balanced by the vertical hydrodynamic force exerted by the vapour layer on the sphere; and (iv) the steady heat conduction problem in the vapour layer which determines the distribution of the temperature gradient in (2.1).

The above problem involves three independent dimensionless groups:

(2.3a–c) \begin{equation}B=\frac{{\rho}_b^*g^* {a^*}^2}{\gamma^*}, \quad E=\frac{\mu^*\lambda^*}{\gamma^*a^*}, \quad D=\frac{{\rho_s^*}}{{\rho}_b^*},\end{equation}

where $\lambda^*=k^*\Delta T^* / \rho^*l^*$ . Here $B={a^*}^2/{\ell_c^*}^2$ is the Bond number, namely the ratio squared of the sphere radius to the liquid bath capillary length $\ell_c^*=(\gamma^*/\rho_b^*g^*)^{1/2}$ , E is an evaporation number and D is a density ratio. (It will be convenient to refer to E as an evaporation number regardless of whether vapour is being produced by evaporation or sublimation.) In the inverse-Leidenfrost scenario where an inert sphere levitates above an evaporating liquid nitrogen bath, typical values are $\ell_c^* \approx 1$ mm and $a^* E \approx 10$ nm [Reference Adda-Bedia, Kumar, Lechenault, Moulinet, Schillaci and Vella2].

2.2 Unified small-slope model

Our interest is in the small-sphere and weak-evaporation regime where B and E are both small, while D is of order unity. Below, we intuitively formulate a ‘unified small-slope’ lubrication model to describe this regime. This model will serve as a starting point for a systematic asymptotic analysis in subsequent sections.

We assume that the slope of the bath height $dh^*/dr^*$ is everywhere small. The dynamical condition at the liquid–vapour interface can then be approximated as

(2.4) \begin{equation}p^*=-\rho_b^*g^*(h^*-{h_{\infty}}^*)+\frac{{\gamma^*}}{r^*}\frac{d}{dr^*}\left(r^*\frac{dh^*}{dr^*}\right),\end{equation}

where $p^*$ is the vapour pressure field evaluated at the liquid–vapour interface. Considering (2.4) as an equation for the bath height $h^*$ , we have the boundary condition

(2.5) \begin{equation}\frac{dh^*}{dr^*}=0 \quad \text{at} \quad r^*=0,\end{equation}

which follows from axial symmetry, and the far-field condition

(2.6) \begin{equation}h^*\to h_{\infty}^* \quad \text{as} \quad r^*\to\infty.\end{equation}

Since (2.4) involves $p^*$ , we need to add equations governing the vapour flow in order to close the problem for $h^*$ . In particular, we anticipate a thin-film region whose radial extent, say of order $\ell_f^*$ , is large relative to the thickness of the film $t^*$ yet small relative to $a^*$ . In that region, the sphere boundary can be approximated by a paraboloid such that

(2.7) \begin{equation}t^*=\frac{1}{2a^*}{r^*}^2 - h^*,\end{equation}

where $dt^*/dr^*\ll1$ follows as a consequence of (2.7) and the above assumptions.

Given that both $h^*$ and $t^*$ are slowly varying, we can employ a standard lubrication approximation to describe the vapour flow in the thin-film region. Consistent with this approximation, we use $\hat{\mathbf{n}}\boldsymbol{\cdot} \boldsymbol{\nabla} T^*\approx \Delta T^*/t^* $ in the boundary condition (2.1) since the vapour temperature varies approximately linearly across the thin film (for a detailed derivation, see [Reference Pomeau, Le Berre, Celestini and Frisch26]). We thereby find the Reynolds, or flux conservation, equation

(2.8) \begin{equation}\frac{1}{r^*}\frac{d}{dr^*}\left(r^*q^*\right)=\frac{\lambda^*}{t^*},\end{equation}

where the flux density $q^*$ is defined as

(2.9) \begin{equation}q^*=-\frac{{t^*}^3}{12\mu}\frac{dp^*}{dr^*},\end{equation}

and the source term on the right-hand side of (2.8) accounts for vapour production at the phase-changing boundary. (Note that the net volumetric flux out of a cylindrical control volume of radius $r^*$ , which is collinear with the symmetry axis, is given by $2 \pi r^* q^*$ .) From axial symmetry, we have the boundary condition

(2.10) \begin{equation}q^*=0 \quad \text{at} \quad r^*=0.\end{equation}

Recall that (2.7)–(2.10) were meant to describe the thin-film region. We anticipate, however, that $p^*$ attenuates fast enough for $r\gg \ell_f^*$ such that the vapour pressure term in the dynamical condition (2.4) can be safely neglected in that region. This suggests extending the lubrication equations (2.7)–(2.10) to all $r^*$ , to be solved in conjunction with the far-field condition

(2.11) \begin{equation}p^*\to0 \quad \text{as} \quad r^*\to\infty\end{equation}

and the approximate vertical-force balance

(2.12) \begin{equation}\int_0^{\infty} p^*r^*\,dr^*=\frac{2}{3}\rho_b^* g^* {a^*}^3.\end{equation}

In the latter balance, the integral represents the vertical pressure force contributed by the thin-film region, assuming that $p^*$ attenuates sufficiently fast that the upper range of integration can be taken to infinity; viscous stresses as well as the slope of the sphere boundary have been neglected consistently with previous assumptions.

We will refer to (2.4)–(2.12) as the unified small-slope model. We use the term ‘unified’ rather than ‘uniform’ as we do not expect the solutions to constitute uniform asymptotic approximations in any given limit. Rather, the model is meant to provide solutions that constitute, for any fixed $r^*$ , leading-order approximations for $t^*$ in the small-sphere and weak-evaporation regime. It is a simple matter to verify a posteriori that the asymptotic solutions we shall derive based on the unified small-slope model are indeed consistent with the assumptions on which that model is predicated.

2.3 Dimensionless equations

Dimensional analysis of the unified small-slope model furnishes just two independent dimensionless groups, which can be chosen as the modified Bond and evaporation numbers

(2.13a,b) \begin{equation}\mathcal{B}=D^{1/2}B, \quad \mathcal{E}=\frac{E}{(BD)^3},\end{equation}

respectively. Note that $\mathcal{B}\propto {a^*}^2$ and $\mathcal{E}\propto {1/{a^*}^7}$ , when varying the sphere radius $a^*$ with all other physical parameters fixed. In what follows, we adopt a dimensionless convention where the radial coordinate $r^*$ is normalised by $(BD)^{1/2}a^*$ , the height of the liquid bath $h^*$ , the sphere displacement $h^*_{\infty}$ and the thickness $t^*$ by $BDa^*$ , the vapour pressure $p^*$ by $\gamma^*/a^*$ and the flux density $q^*$ by $(BD)^{5/2}\gamma^{*} a^{*}/\mu^*$ ; the associated normalised quantities are denoted similarly but with the asterisk omitted. In particular, the relation (2.7) between the thickness and the bath height becomes

(2.14) \begin{equation}t=\frac{1}{2}r^2-h.\end{equation}

A dimensionless schematic is shown in Figure 1(b). We note that the modified evaporation number $\mathcal{E}$ differs from those used in preceding analyses of the Leidenfrost effect [Reference Maquet, Sobac, Darbois-Texier, Duchesne, Brandenbourger, Rednikov, Colinet and Dorbolo24, Reference Sobac, Rednikov, Dorbolo and Colinet30, Reference van Limbeek, Sobac, Rednikov, Colinet and Snoeijer31].

In dimensionless form, the unified small-slope model consists of the lubrication equations

(2.15a,b) \begin{equation}\frac{1}{r}\frac{d}{dr}\left(rq\right)=\frac{\mathcal{E}}{t}, \quad q=-\frac{{t}^3}{12}\frac{dp}{dr};\end{equation}

the dynamic condition

(2.16) \begin{equation}p=-\mathcal{B}^2\left(h-h_{\infty}\right)+\frac{1}{r}\frac{d}{dr}\left(r\frac{dh}{dr}\right),\end{equation}

which in terms of thickness can be written

(2.17) \begin{equation}p=2-\frac{\mathcal{B}^2}{2}r^2+\mathcal{B}^2\left(t+h_{\infty}\right)-\frac{1}{r}\frac{d}{dr}\left(r\frac{dt}{dr}\right);\end{equation}

the boundary conditions

(2.18a,b) \begin{equation}\frac{dh}{dr}=0 \,\, \left(\text{or }\frac{dt}{dr}=0\right), \quad q=0 \quad \text{at} \quad r=0;\end{equation}

the far-field conditions

(2.19a,b) \begin{equation}p\to0, \quad h\to h_{\infty} \quad \text{as} \quad r\to\infty;\end{equation}

and the vertical force balance

(2.20) \begin{equation}\int_0^{\infty}pr\,dr=\frac{2}{3}.\end{equation}

In the remainder of the paper, we asymptotically analyse the unified small-slope model in the small-sphere regime $\mathcal{B}\ll1$ . We shall see that this regime is comprised of formally infinite sequences of distinguished limits, the limits being defined by the smallness of $\mathcal{E}$ relative to $\mathcal{B}$ . We will initially assume that the modified evaporation number $\mathcal{E}$ is order unity and find that the unified small-slope model breaks down for large enough $\mathcal{E}$ . Subsequently, we shall consider limits corresponding to relatively smaller and smaller $\mathcal{E}$ .

We note that the unified small-slope model is fairly straightforward to solve numerically, for example using Matlab routines for boundary-value problems [25]. (In order to accurately impose the far-field conditions, it is useful to develop local asymptotics of the fields for large r.) Numerical solutions obtained in this manner will be used throughout the paper for the sake of comparison with asymptotic approximations. At this stage, the reader may wish to refer to Figures 2(a), (b), 6 and 10 showing film profiles calculated for $\mathcal{B}=0.1$ and $\mathcal{E}=0.1,10^{-6},10^{-22}$ and $10^{-29}$ , respectively. With respect to the morphological evolution discussed in the introduction, Figure 2 demonstrates the initial stage of the evolution where a neck–bubble structure forms, whereas Figure 10 demonstrates the final stage of the evolution where the film consists mostly of a uniform-thickness region which expands radially. Our analysis will theoretically illuminate the linkage between these stages.

Figure 2.

The height profile h(r) for $\mathcal{B}=0.1$ and (a) $\mathcal{E}=0.1$ ; (b) $\mathcal{E}=10^{-6}$ . The numerical solution of the unified small-slope model (lower solid line) is compared with asymptotic approximations derived in Section 3 in the first distinguished limit for the thin film region (dashed line) and hydrostatic region (dash-dotted line). The sphere, locally approximated by a paraboloid, is represented by the upper solid line. As $\mathcal{E}$ decreases, the sphere levitates closer to the bath and the thin film decomposes into a bubble region bounded by a narrow and thin neck region (see Section 3.3).

While values of $\mathcal{E}$ smaller than around $10^{-6}$ are unrealistic (see estimates below (2.3)), the need to consider exceedingly small $\mathcal{E}$ in order to observe the late stages of the evolution is a consequence of fixing $\mathcal{B}$ in the above numerical examples. We do this so that the formally small parameters in our analysis, $\mathcal{B}$ and $\mathcal{E}$ , remain numerically small throughout the evolution. Alternatively, a similar evolution can be observed by increasing the sphere radius $a^*$ with all other parameters fixed, in which case $\mathcal{B}$ increases while $\mathcal{E}$ decreases (as shown in [Reference Maquet, Sobac, Darbois-Texier, Duchesne, Brandenbourger, Rednikov, Colinet and Dorbolo24] for drops levitated above a liquid substrate). While in that scenario $\mathcal{E}$ and $\mathcal{B}$ maintain physically reasonable values, $\mathcal{B}$ is not so small during the later stages of the evolution. Thus, we expect only a qualitative correspondence between our asymptotic analysis and the later stages of the evolution as they would be observed numerically or experimentally for typical values of the physical parameters. Moreover, in that scenario our unified small-slope model is unlikely to be accurate since, as will become clear, the assumptions underlying it imply $\mathcal{B}\ll1$ .

3.1 Thin-film region

In this section we consider the unified small-slope model in the limit $\mathcal{B}\to0$ , with $\mathcal{E}\simeq 1$ . (Henceforth, we use $\simeq$ to denote asymptotic order.) We shall later refer to this regime as the first distinguished limit.

We start by posing the straightforward expansions

(3.1a–d) \begin{equation} h=h_0(r)+\cdots, \quad t=t_0(r)+\cdots, \quad q= q_0(r) + \cdots, \quad p= p_0(r) + \cdots,\end{equation}

where, from (2.14),

(3.2) \begin{equation}t_0=\frac{1}{2}r^2-h_0.\end{equation}

We now list the equations satisfied by the leading-order fields. From (2.15) and (2.16),we have

(3.3a–c) \begin{equation}\frac{1}{r}\frac{d}{dr}\left(rq_0\right)=\frac{\mathcal{E}}{t_0}, \quad q_0=-\frac{t_0^3}{12}\frac{dp_0}{dr}, \quad p_0=\frac{1}{r}\frac{d}{dr}\left(r\frac{dh_0}{dr}\right),\end{equation}

which are identical to the equations of the unified small-slope model except for the absence of gravity terms in the approximate dynamic condition (3.3c). From (2.18), we have the boundary conditions

(3.4a,b) \begin{equation}\frac{dh_0}{dr}=0, \quad q_0=0 \quad \text{at} \quad r=0,\end{equation}

while (2.19a) gives the attenuation condition

(3.5) \begin{equation}p_0\to0 \quad \text{as} \quad r\to\infty.\end{equation}

We have yet to refer to the far-field condition (2.19b) on the bath height and to the force constraint (2.20). With respect to the former, it is noted that the attenuation of $p_0$ together with (3.3c) implies the far-field behaviour

(3.6) \begin{equation}h_0 \sim A\ln r \quad \text{as} \quad r\to\infty,\end{equation}

where A is a constant. With $h_0$ growing logarithmically, the far-field condition (2.19b) cannot be satisfied. The issue is that at large distances, comparable to the capillary length, or $r\simeq 1/\mathcal{B}$ in our normalisation scheme, the gravity terms neglected in (2.16) become important. Expansions (3.1) are accordingly interpreted as ‘outer’ approximations, i.e., for r fixed, which which describe the thin-film region. In principle, the value of A in (3.6) comes from matching with an ‘inner’ hydrostatic region. It is clear, however, that the pressure integral in the force balance (2.20) is dominated by the thin-film region, whereby (2.20) yields

(3.7) \begin{equation}\int_0^\infty p_0 r\,dr=\frac{2}{3}.\end{equation}

Hence, A can be readily calculated by substituting (3.3c) and (3.6) into (3.7). This gives

(3.8) \begin{equation}A=\frac{2}{3}.\end{equation}

The leading-order thin-film problem is now closed; it is, in fact, mathematically identical to that formulated in [Reference Pomeau, Le Berre, Celestini and Frisch26] for a small Leidenfrost drop levitated above a flat solid substrate. In the present problem, however, it remains to satisfy the far-field condition (2.19b) on the bath height; this does not imply an overdetermined problem as that condition involves the yet unknown sphere displacement $h_{\infty}$ . In the next subsection, we calculate $h_{\infty}$ by matching the thin-film region with the hydrostatic region.

It is straightforward to solve the leading-order thin-film problem numerically. In particular, unlike the unified small-slope model, the present reduced problem does not involve multiple length scales, at least for values of $\mathcal{E}$ on the order of unity (see Section 3.3). A comparison between the two schemes is presented in Figure 2, where the bath-height profile is plotted for $\mathcal{B}=0.1$ and two values of $\mathcal{E}$ , $0.1$ and $10^{-6}$ . Further comparison is presented in Figure 3, which shows the variation with $\mathcal{E}$ of the film thickness at the origin.

Figure 3.

Scaled film thickness at the origin as a function of $\mathcal{E}$ , for $\mathcal{B}=0.1$ . Square symbols: numerical solution of the unified small-slope model (Section 2.3). Dash-dotted curve: approximation in the first distinguished limit (Section 3.1). Dashed curve: approximation in the second distinguished limit (Section 4.1). Dotted curve: approximation in the third distinguished limit (Section 4.3). Solid curve: approximation valid throughout the second formally infinite sequence of distinguished limits (Section 6).

We will see that the specific output of the thin-film problem that is required for matching with the hydrostatic region is the limit

(3.9) \begin{equation}\beta = \lim_{r\to\infty}\left(h_0-\frac{2}{3}\ln r\right),\end{equation}

which corresponds to the constant order-unity correction to the far-field expansion (3.6). Like the thin-film problem itself, $\beta$ is a function of $\mathcal{E}$ alone. Analytical approximations for $\beta$ can be deduced for small and large $\mathcal{E}$ . In the latter limit, corresponding to strong evaporation, or relatively small spheres, the bath height is approximately constant and equal to $\beta$ . The geometry is then simply that of a sphere above a vertically displaced flat bath. In that case, theapproximation

(3.10) \begin{equation}\beta \sim -\frac{3}{2}\mathcal{E}^{1/2} \quad \text{as} \quad \mathcal{E}\to\infty\end{equation}

follows from the analysis in [Reference Pomeau, Le Berre, Celestini and Frisch26] of an analogous limit (there appears to be a typo in the approximation obtained therein). The diametric limit $\mathcal{E}\to0$ , which corresponds to weak evaporation, or relatively large spheres, will be discussed in Section 3.3 and in a wider context for the remainder of this paper. For the moment, we simply state the result

(3.11) \begin{equation}\beta \to \frac{1}{3}\left(\ln\frac{3}{2}+1\right) \quad \text{as} \quad \mathcal{E}\to0,\end{equation}

which will fall out from the analysis in Section 4. (As explained in [Reference Sobac, Rednikov, Dorbolo and Colinet30], the weak-evaporation analysis in [Reference Pomeau, Le Berre, Celestini and Frisch26] is erroneous. Thus, (3.11) cannot be extracted from that analysis.) The function $\beta(\mathcal{E})$ , as obtained from numerically solving the thin-film problem, is shown in Figure 4 along with the analytical approximations (3.10) and (3.11).

Figure 4.

Sphere displacement $h_{\infty}$ as a function of $\mathcal{E}$ , after subtracting the $\mathcal{B}$ dependence anticipated from the asymptotic approximation (3.18). Numerical solutions of the unified small-slope model for several values of $\mathcal{B}$ (symbols) collapse on the asymptotic prediction based on (3.18), which following the subtraction equals the function $\beta(\mathcal{E})$ defined by the thin-film problem of Section 3.1. The large- $\mathcal{E}$ approximation (3.10) and small- $\mathcal{E}$ approximation (3.11) of that function are depicted by the dashed and dash-dotted curves, respectively. The latter approximation also describes a small non-wetting sphere floating on a liquid bath.

3.2 Hydrostatic region, matching and sphere displacement

To study the hydrostatic region we define the strained coordinate

(3.12) \begin{equation}\hat{r}=\mathcal{B}r.\end{equation}

Given (3.6), it is clear that $h\simeq 1$ in the hydrostatic region so that the gravity and capillarity terms in the dynamic condition (2.16) are $\simeq \mathcal{B}^2$ . In contrast, since $p_0\simeq 1/r^6$ as $r\to\infty$ , as readily follows from local analysis of the thin-film problem, the pressure term in (2.16) is $\simeq \mathcal{B}^6$ . (These scalings are accurate to within logarithmic factors.) Thus, in the hydrostatic region the dynamic condition provides an equation for the bath height which is uncoupled from the vapour flow, except through asymptotic matching.

In light of the above, we pose the expansions

(3.13a,b) \begin{equation}h = \hat{h}(\hat{r}) + \cdots, \quad h_{\infty} = \hat{h}_{\infty} + \cdots.\end{equation}

The leading-order field $\hat{h}$ satisfies

(3.14) \begin{equation}\frac{1}{\hat{r}}\frac{d}{d\hat{r}}\left(\hat{r}\frac{d\hat{h}}{d\hat{r}}\right)-\hat{h} = -\hat{h}_{\infty},\end{equation}

which follows from the dynamic condition (2.16), the far-field condition

(3.15) \begin{equation}\hat{h}\to \hat{h}_{\infty} \quad \text{as} \quad \hat{r}\to\infty,\end{equation}

which follows from (2.19b), and matching conditions in the limit $\hat{r}\to0$ .

The solution to (3.14) and (3.15) consists of a superposition of the constant particular solution $\hat{h}_{\infty}$ and a homogeneous solution that decays as $\hat{r}\to\infty$ . The latter solution is proportional to the modified Bessel function of the second kind $K_0(\hat{r})$ that possesses the asymptotic behaviour

(3.16) \begin{equation}K_0(\hat{r})\sim -\ln \hat{r} + \ln 2 - \gamma_E \quad \text{as} \quad \hat{r}\searrow 0,\end{equation}

where $\gamma_E\doteq 0.5772\ldots$ is the Euler–Gamma constant [Reference Abramowitz and Stegun1]. Using this behaviour, straightforward matching between the thin-film and hydrostatic regions implies the solution

(3.17) \begin{equation}\hat{h} = \hat{h}_{\infty} -\frac{2}{3} K_0(\hat{r}),\end{equation}

with the leading-order sphere displacement $\hat{h}_{\infty}$ provided as

(3.18) \begin{equation}\hat{h}_{\infty}=\frac{2}{3}\ln\frac{2}{\mathcal{B}}-\frac{2}{3}\gamma_E + \beta ,\end{equation}

wherein $\beta$ is the function of $\mathcal{E}$ calculated from the thin-film problem (cf. (3.9)).

In Figure 4, we validate the asymptotic prediction (3.18) against numerical solutions of the unified small-slope model. The numerical values for the sphere displacement, for several small values of $\mathcal{B}$ , are seen to collapse on the asymptotic prediction once the $\mathcal{B}$ dependence anticipated from (3.18) is subtracted. Analytical approximations for the sphere displacement in the limits of large and small $\mathcal{E}$ can be obtained by substituting the corresponding approximations (3.10) and (3.11) for $\beta$ into (3.18). In particular, substituting the latter gives the limit

(3.19) \begin{equation}\hat{h}_{\infty}\to \frac{2}{3}\ln\frac{{\sqrt{6}}}{\mathcal{B}}+\frac{1}{3}-\frac{2}{3}\gamma_E \quad \text{as} \quad \mathcal{E}\to0,\end{equation}

which corresponds to the displacement of a small non-wetting sphere floating on a liquid bath [Reference Cooray, Cicuta and Vella11, Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris15].

We conclude this subsection with the following remarks:

1. (1) Note the appearance in (3.18) of the logarithm of the small parameter $\mathcal{B}$ . Except where stated otherwise, we shall collect together terms that are of the same algebraic order in $\mathcal{B}$ . When technically possible, this convention allows obtaining asymptotic approximations that are algebraically, rather than logarithmically, accurate, without the need to carry out the analysis at successive logarithmic orders [Reference Hinch19]. Consistently with this convention, the symbol $\simeq$ should be interpreted as asymptotic order ignoring logarithmic factors.

2. (2) At this stage, it may not be clear why we introduced a unified small-slope model, only to asymptotically decompose that model at the first opportunity into thin-film and hydrostatic regions. We could have also arrived at this picture by a direct matched asymptotics analysis of the original physical problem formulated in Section 2.1. Nonetheless, the unified formulation will prove useful in the remainder of the paper, where we will study a large number of distinguished limits associated with the relative smallness of $\mathcal{B}$ and $\mathcal{E}$ , the analysis in this section corresponding to just the first of these limits. Analysing each of these limits directly from the original formulation is technically possible but very tedious and repetitive.

3. (3) As discussed in Section 3.4, the reduced model derived in the present section loses validity at sufficiently small $\mathcal{E}$ . Nonetheless, the result (3.18) for the sphere displacement, including its limiting value (3.19), remains valid. This is confirmed by the analysis in the subsequent sections and is also evident from the comparison presented in Figure 4.

4. (4) From (3.10) it can be deduced that $h_0\simeq\mathcal{E}^{1/2}$ for large $\mathcal{E}$ . It then follows from (3.2) that $dt_0/dr\simeq\mathcal{E}^{1/4}$ in that limit. Using the scaling factors in Section 2.3, the corresponding scaling of $dt^*/dr^*$ is $(BD)^{1/2} \mathcal{E}^{1/4}$ . For D of order unity, the latter scaling implies that the unified small-slope model breaks down for $ \mathcal{E}\simeq 1/\mathcal{B}^{2}$ .

3.3 Formation of neck–bubble morphology

Consider the thin-film problem in the limit $\mathcal{E}\to0 $ . As seen in Figure 2, the vapour film decomposes into a very thin central bubble region bounded by a narrow annular neck region which is even thinner. This morphology is analogous to that familiar in the classical Leidenfrost effect [Reference Sobac, Rednikov, Dorbolo and Colinet30]. We note, in particular, that for small $\mathcal{E}$ the origin becomes a maximum of the thickness profile.

Let $L_B$ and $T_B$ be the orders of the radial extent and thickness of the vapour bubble, respectively. Assuming that $T_B\ll L_B^2$ , namely that the bubble is thin compared to the vertical variation of the sphere boundary over that region, the pressure in the bubble region must approximately equal the capillary pressure associated with the sphere’s curvature, i.e.,

(3.20) \begin{equation}p_0 \sim 2.\end{equation}

Assuming that the bubble region dominates the pressure integral in the force balance (3.7), we have that $L_B=1$ ; specifically, the vapour bubble must terminate at $r=l$ , where

(3.21) \begin{equation}l= \sqrt{2/3}.\end{equation}

Let $Q_B$ , P $^{\prime}$ and $\Delta P$ be the orders of the flux density, pressure correction and pressure variation, respectively, in the bubble region. With reference to the flux-conservation equation (3.3a), it is clear that the vapour production term on the right-hand side of that equation cannot be negligible in the bubble region; otherwise, the flux density would be singular at the origin. Thus, the thin-film equations (3.3) imply

(3.22a–c) \begin{equation}Q_B=\frac{\mathcal{E}}{T_B}, \quad P'=T_B, \quad \Delta P=\frac{\mathcal{E}}{T_B^5},\end{equation}

where the scaling $T_B\ll1$ remains to be determined.

Consider now the neck region at the edge of the bubble region. Denote the orders of the neck’s radial width, thickness and flux density by $L_N$ , $T_N$ and $Q_N$ , respectively. Since $p_0$ is clearly order unity and not constant in the neck region, (3.3b) and (3.3c) imply

(3.23a,b) \begin{equation}T_N=L_N^2, \quad Q_N=L_N^5.\end{equation}

It is also evident that $L_N\ll1$ and $T_N\ll T_B$ .

We now relate the bubble and neck scalings based on matching considerations. First, using $L_N\ll1$ , the dynamical condition (3.3c) simplifies in the neck region as

(3.24) \begin{equation}\frac{d^2t_0}{dr^2}\approx 2-p_0.\end{equation}

Since $p_0\to 2$ towards the bubble region (cf. (3.20)), $t_0$ grows linearly in that direction. Matching $dt_0/dr$ thus gives $T_B/L_B=T_N/L_N$ , i.e., $T_B=L_N$ . Second, we argue that $Q_B=Q_N$ . Without vapour production, this would be obvious from continuity. With vapour production, inspection of the flux-conservation equation (3.3a) together with the fact that the bubble thickness vanishes linearly towards the neck means that $Q_B$ and $Q_N$ can differ only by a factor of logarithmic order in $\mathcal{E}$ . Thus, up to such logarithmic factors, we find

(3.25a–c) \begin{equation}L_N = \mathcal{E}^{1/6}, \quad T_N = \mathcal{E}^{1/3}, \quad Q_N = \mathcal{E}^{5/6},\end{equation}

and

(3.26a–d) \begin{equation}T_B=\mathcal{E}^{1/6}, \quad Q_B=\mathcal{E}^{5/6}, \quad P^{\prime}_B=\mathcal{E}^{1/6}, \quad \Delta P_B = \mathcal{E}^{1/3}.\end{equation}

It follows from the above scalings that the pressure correction $p_0-2$ in the bubble region is approximately uniform. Furthermore, vapour production in the neck is not negligible.

In principle, the decomposition of the thin-film problem at small $\mathcal{E}$ includes a third region, exterior to the neck region. We shall refer to it as the transition region as it bridges the neck region and the far-field behaviour (3.6). We denote the orders of the radial extent, thickness, flux density and pressure in the transition region by $L_T,T_T,Q_T$ and $P_T$ , respectively. Anticipating that $p_0$ decreases when crossing the neck in the outward radial direction, the approximate dynamical condition (3.24) implies that $t_0$ grows quadratically in that direction away from the neck. Order-of-magnitude matching of $d^2t_0/dr^2$ accordingly gives $T_T/L_T^2=T_N/L_N^2$ , or $T_T=L_T^2$ . Together with the height-thickness relation (3.2) and far-field condition (3.6), this implies $L_T=T_T=1$ . Similarly, matching the flux density gives $Q_T=\mathcal{E}^{5/6}$ , whereby the flux-pressure relation (3.3b) gives $P_T=\mathcal{E}^{5/6}$ . For later reference, we summarise the transition-region scalings:

(3.27a–d) \begin{equation}L_T=1, \quad T_T=1, \quad Q_T=\mathcal{E}^{5/6}, \quad P_T=\mathcal{E}^{5/6}.\end{equation}

Note that the pressure in the transition region is small, whereby at order unity the pressure in the neck region attenuates rapidly from its value (3.20) in the bubble to zero. This is consistent with our tacit assumption that the force condition (3.7) is dominated by the bubble region. Furthermore, substituting the scalings (3.27) into the flux conservation equation (3.3a) confirms that vapour production is negligible in the transition region, whereby the order $\mathcal{E}^{5/6}$ flux density originates from the vapour production in the bubble and neck regions.

3.4 Breakdown of first distinguished limit

Recall that the above scalings, which imply the formation of a neck–bubble morphology for small $\mathcal{E}$ , were derived from the leading-order thin-film problem of Section 3.1. The latter problem, in turn, was derived from the unified small-slope model by considering the limit $\mathcal{B}\to0$ with $\mathcal{E}\simeq1$ . This suggests returning to the unified small-slope model to check, based on the scalings derived, whether new physics not included in the thin-film problem of Section 3.1 become important for sufficiently small $\mathcal{E}$ .

In particular, consider the dynamical condition (2.17) written in terms of the thickness profile t. In the bubble region, the small- $\mathcal{E}$ scalings imply a dominant balance at order $\mathcal{E}^{1/6}$ which includes a uniform pressure correction and the curvature of the thickness profile. Given the smallness in $\mathcal{E}$ of these leading terms, it is plausible that the gravity terms in (2.17), neglected in the leading-order thin-film problem on account of their smallness in $\mathcal{B}$ , enter this dominant balance for sufficiently small $\mathcal{E}$ . Indeed, using the scaling $L_B=1$ and the fact that $\hat{h}_{\infty}\simeq 1$ as $\mathcal{E}\to0$ , we have the estimate $\mathcal{B}^2(h_{\infty}-r^2/2)\simeq \mathcal{B}^2$ ; since $T_B\ll1$ , the remaining gravity term $\mathcal{B}^2 t$ is relatively negligible. Thus, the order $\mathcal{B}^{2}$ gravity terms in (2.17) are not negligible in the order $\mathcal{E}^{1/6}$ balance of the dynamic condition for

(3.28) \begin{equation}\mathcal{E}\simeq \mathcal{B}^{12}.\end{equation}

For such small $\mathcal{E}$ , the thin-film problem of Section 3.1 no longer holds – the first distinguished limit of the unified small-slope model breaks down. This process is schematically shown in Figure 5. We note that a similar breakdown, wherein gravity terms become important in the vapour film, does not occur in the case of small Leidenfrost drops levitating above a flat solid substrate.

Figure 5.

Depiction of the first four distinguished limits (see Sections 3 and 4). The bold letters C, P, E and G indicate the importance of capillarity, pressure-driven flow, vapour production and variations in gravitational hydrostatic pressure, respectively, for determining the leading-order thickness profile; furthermore, bold L indicates a radially localised region. The blue arrows and dashed curves indicate evolution of the film morphology with decreasing $\mathcal{E}$ .

4.1 Second distinguished limit

We are led to analyse a second distinguished limit of the unified small-slope model, where both $\mathcal{B}$ and $\mathcal{E}$ are small in accordance with (3.28). We expect the thin-film region to decompose into bubble, neck and transition regions as described in the preceding section, only with gravity now having a leading-order effect in the bubble region. Accordingly, we expect the small- $\mathcal{E}$ scalings derived in Section 3.3, based on the first limit, to remain valid in the second. Furthermore, their analysis in the second limit should overlap with the small- $\mathcal{E}$ behaviour of the first limit, which is why we did not go beyond a scaling analysis in Section 3.3.

In analysing the second distinguished limit, it will be convenient to work with the single small parameter $\mathcal{B}$ . We accordingly define the rescaled evaporation number

(4.1) \begin{equation}\mathcal{\tilde{E}}= \mathcal{E}/\mathcal{B}^{12}.\end{equation}

Similarly, we use (3.28) to rewrite the scalings derived in Section 3.3 in terms of $\mathcal{B}$ . Thus, the bubble scalings become

(4.2a–e) \begin{equation}L_B=1, \quad T_B=\mathcal{B}^2, \quad Q_B=\mathcal{B}^{10}, \quad P_B^{\prime}=\mathcal{B}^2, \quad \Delta P_B=\mathcal{B}^4,\end{equation}

the neck scalings become

(4.3a–c) \begin{equation}L_N = \mathcal{B}^2, \quad T_N = \mathcal{B}^4, \quad Q_N=\mathcal{B}^{10}\end{equation}

and the transition scalings become

(4.4a–d) \begin{equation}L_T=1, \quad T_T=1, \quad Q_T=\mathcal{B}^{10}, \quad P_T=\mathcal{B}^{10}.\end{equation}

4.1.1 Bubble region

We start with the bubble region, where the scalings (4.2) suggest the expansions

(4.5a–c) \begin{equation}t= \mathcal{B}^{2}\kern2pt\tilde{\kern-2pt t}(r) + \cdots, \quad q = \mathcal{B}^{10}\tilde{q}(r) + \cdots,\quad p =2+ \mathcal{B}^{2}\tilde{p} + \cdots,\end{equation}

where $\tilde{p}$ is constant and $0\le r < l$ (cf. (3.21)). Expanding the sphere displacement as in (3.13b), the leading non-trivial balance of the dynamical condition (2.17) yields

(4.6) \begin{equation}\frac{1}{r}\frac{d}{dr}\left(r\frac{d\kern2pt\tilde{\kern-2pt t}}{dr}\right)=\hat{h}_{\infty}-\frac{1}{2}r^2-\tilde{p},\end{equation}

which can be considered an equation for $\kern2pt\tilde{\kern-2pt t}$ . The first two terms on the right-hand side together account for the effect of gravity owing to the displacement between the local film height and the unperturbed bath. This equation is supplemented by the boundary conditions

(4.7a,b) \begin{equation}\frac{d\kern2pt\tilde{\kern-2pt t}}{dr}=0 \quad \text{at} \quad r=0, \qquad \kern2pt\tilde{\kern-2pt t}=0 \quad \text{at} \quad r=l,\end{equation}

where the first follows from (2.18a) and (2.14), while the second represents matching with the relatively thin neck region. Solving (4.6) together with (4.7) gives

(4.8) \begin{equation}\kern2pt\tilde{\kern-2pt t}=\frac{1}{32}\left(l^2-r^2\right)\left(r^2+\alpha\right),\end{equation}

in which

(4.9) \begin{equation}\alpha = l^2+8\left(\tilde{p}-\hat{h}_{\infty}\right)\end{equation}

is a constant to be determined together with $\tilde{p}$ and $\hat{h}_{\infty}$ . Note that $\alpha$ must be positive as otherwise $\kern2pt\tilde{\kern-2pt t}$ would vanish within the bubble domain.

Next, the leading-order balance of the flux-conservation equation (2.15a), together with (4.1), gives

(4.10) \begin{equation}\frac{1}{r}\frac{d}{dr}\left(r\tilde{q}\right)=\frac{\mathcal{\tilde{E}}}{\kern2pt\tilde{\kern-2pt t}},\end{equation}

to be solved for $\tilde{q}$ in conjunction with the boundary condition

(4.11) \begin{equation}\tilde{q}=0 \quad \text{at} \quad r=0,\end{equation}

which follows from (2.18b). Substituting (4.8) and integrating using (4.11) yields

(4.12) \begin{equation}\tilde{q} = \frac{16\mathcal{\tilde{E}}}{(l^2+\alpha) r}\ln \frac{(1+ r^2/\alpha)}{1-r^2/l^2}.\end{equation}

For later reference, we note the behaviour

(4.13) \begin{equation}\tilde{q}\sim -\frac{8\sqrt{6}\mathcal{\tilde{E}}}{l^2+\alpha}\ln\frac{\sqrt{6}\left(l-r\right)\alpha}{l^2+\alpha} \quad \text{as} \quad r\nearrow l,\end{equation}

where the asymptotic error is algebraic in the indicated limit.

4.1.2 Neck region

Consider next the neck region near $r=l$ . With reference to the scalings (5.8), we introduce the strained coordinate

(4.14) \begin{equation}R=(r-l)/\mathcal{B}^2\end{equation}

and pose the expansions

(4.15a–c) \begin{equation}t=\mathcal{B}^4T(R)+\cdots, \quad q=\mathcal{B}^{10}Q(R)+\cdots, \quad p=P(R)+\cdots.\end{equation}

From (2.15) and (2.17), we find the governing equations

(4.16a–c) \begin{equation}\frac{d{Q}}{dR}=\frac{\mathcal{\tilde{E}}}{{T}}, \quad {Q}=-\frac{{T}^3}{12}\frac{d{P}}{d R}, \quad{P}=2-\frac{d^2{T}}{d R^2}.\end{equation}

These are to be solved in conjunction with matching conditions as $R\to\pm\infty$ to be discussed below.

It turns out that the parameter $\tilde{\mathcal{E}}$ can be factored out from both (4.16) and the matching conditions by use of the transformations

(4.17a–c) \begin{equation}T(R)=\tilde{\mathcal{E}}^{1/3}T'(R'), \quad Q(R)=\tilde{\mathcal{E}}^{5/6}Q'(R'), \quad P(R)=P'(R'),\end{equation}

where

(4.18) \begin{equation}R=R_s+\tilde{\mathcal{E}}^{1/6}R'\end{equation}

and $R_s$ is a constant shift discussed below. In particular, (4.16) become

(4.19a–c) \begin{equation}\frac{dQ'}{dR'}=\frac{1}{T'}, \quad Q'=-\frac{T'^3}{12}\frac{dP'}{dR'}, \quad P'=2-\frac{d^2T'}{d R'^2}.\end{equation}

Consider now the conditions as $R\to\pm\infty$ . Straightforward matching of the pressure field between the neck and bubble regions gives

(4.20) \begin{equation}P'\to2 \quad \text{as} \quad R'\to-\infty,\end{equation}

while the asymptotically small magnitude of the pressure in the transition region implies

(4.21) \begin{equation}P'\to0 \quad \text{as} \quad R'\to\infty.\end{equation}

With (4.20) and (4.21), local analysis of (4.19) as $R\to\pm\infty$ yields the behaviours

(4.22a,b) \begin{equation}T'\sim -\theta R' , \quad Q'\sim -\frac{1}{\theta}\ln |R'| + \tau \quad \text{as} \quad R'\to-\infty\end{equation}

and

(4.23a,b) \begin{equation}T'= R'^2 + O(1), \quad Q'\sim \nu \quad \text{as} \quad R'\to\infty,\end{equation}

where $\theta, \tau$ and $\nu$ are constants to be determined, and we have used the displacement invariance associated with the constant shift $R_s$ to eliminate the linear term in (4.23a); that choice defines the constant $R_s$ , which we will calculate later. The transformed neck problem consisting of (4.19), (4.22) and (4.23) determines two relations between the three unknowns $\theta, \tau$ and $\nu$ , which we conceptually represent by functions: $\tau(\theta)$ and $\nu(\theta)$ . In terms of these functions, asymptotic matching of the flux density between the bubble and neck regions provides the relations

(4.24a,b) \begin{equation}\theta=\frac{2+3\alpha}{24\sqrt{6}\tilde{\mathcal{E}}^{1/6}}, \quad \theta \tau =\ln \frac{8 \theta}{\alpha \mathcal{B}^2}.\end{equation}

Numerical solution of the transformed neck problem together with the conditions (4.24) determines $\alpha$ and $\theta$ (and hence also $\tau$ and $\nu$ ) as functions of $\tilde{\mathcal{E}}$ and $\ln \mathcal{B}$ .

As an alternative to the above numerical procedure, further asymptotic reduction can be achieved by expanding the solutions in inverse logarithmic powers (recall our convention to group together terms of the same algebraic order). We do not pursue this here as the resulting approximations are inaccurate while not adding significant insight.

4.1.3 Transition and hydrostatic regions

We continue with the transition region, where the scalings (3.27) suggest the expansions

(4.25) \begin{equation}h= \acute{h}(r)+\cdots, \quad t =\acute{t}(r) + \cdots,\end{equation}

where $r>l$ . From (2.14), the leading-order thickness and height profiles are related as

(4.26) \begin{equation}\acute{t}=\frac{1}{2}r^2-\acute{h}.\end{equation}

Given the smallness of p (cf. (3.27d)), the dynamical condition (2.16) gives

(4.27) \begin{equation}\frac{1}{r}\frac{d}{dr}\left(r\frac{d\acute{h}}{dr}\right)=0.\end{equation}

Boundary conditions on $\acute{h}$ can be derived by asymptotic matching between the neck and transition regions. Thus, matching of the thickness profile gives the conditions

(4.28a,b) \begin{equation}\acute{t}=0, \quad \frac{d\acute{t}}{dr}=0 \quad \text{at} \quad r=l,\end{equation}

which are intuitive given the quadratic widening of the vapour film outward from the thin neck. In terms of $\acute{h}$ , (4.28) read as

(4.29a,b) \begin{equation}\acute{h}=\frac{1}{3}, \quad \frac{d\acute{h}}{dr}=\sqrt{\frac{2}{3}} \quad \text{at} \quad r=l.\end{equation}

Solving (4.27) together with (4.29), we find

(4.30) \begin{equation}\acute{h}=\frac{2}{3}\ln r + \frac{1}{3}\left(1+\ln\frac{3}{2}\right).\end{equation}

The analysis of the hydrostatic region is similar to that in the first distinguished limit. In particular, the bath height h is expanded as in (3.13a), with the leading-order field $\hat{h}$ satisfying the differential equation (3.14) and far-field condition (3.15). Together with asymptotic matching with the transition region (cf. (4.30)), we find the same solution (3.17) but with the leading-order sphere displacement given by

(4.31) \begin{equation}\hat{h}_{\infty} = \frac{2}{3}\ln\frac{{\sqrt{6}}}{\mathcal{B}}+\frac{1}{3}-\frac{2}{3}\gamma_E.\end{equation}

The value (4.31) is just the small- $\mathcal{E}$ limit (3.19) of the formula (3.18) obtained for $\hat{h}_{\infty}$ in the analysis of the first distinguished limit. Conversely, the limit (3.11) of the function $\beta$ appearing in that formula is readily deduced from (4.30). We have already noted that (3.19) corresponds to the displacement of a non-wetting sphere floating at a liquid–bath interface.

4.1.4 Perturbed bubble pressure and shift of the neck radius

With $\hat{h}_{\infty}$ and $\alpha$ determined, the latter’s definition (4.9) yields the constant correction $\tilde{p}$ to the pressure in the bubble region as a function of $\mathcal{\tilde{E}}$ and $\ln\mathcal{B}$ . Furthermore, carefully enforcing the force constraint (2.20) beyond leading order allows calculating the constant $R_s$ , which corresponds to a small shift of the neck radius (cf. (4.18)). The details are given in Appendix A, where we find

(4.32) \begin{equation}R_s=\frac{1-6\hat{h}_{\infty}}{12\sqrt{6}}\end{equation}

with $\hat{h}_{\infty}$ provided by (4.31). While the value of $R_s$ does not affect the leading-order approximations obtained to this point, knowledge of it is necessary in order to compare our asymptotic approximations for the neck region with numerical solutions of the unified small-slope model.

In Figure 6, we compare the neck and bubble approximation obtained in the second distinguished limit with a numerical solution of the unified small-slope model for $\mathcal{B}=0.1$ and $\mathcal{E}=10^{-22}$ . In plotting the neck profile, we use the result (4.32) for the displacement of the neck position from $r=l$ . Furthermore, the approximation for the film thickness at the origin is depicted in Figure 3 as a function of $\mathcal{E}$ by the dashed curve.

Figure 6.

The film-thickness profile for $\mathcal{B}=0.1$ and $\mathcal{E}=10^{-22}$ . The numerical solution of the unified small-slope model is depicted by the black solid line. Asymptotic predictions in the second distinguished limit (Section 4.1) are depicted by the blue dotted line (bubble-region approximation) and dash-dotted line (neck-region approximation). The inset shows the vicinity of the neck near $r=l$ .

4.2 Breakdown of second distinguished limit

4.2.1 Variation of bubble morphology with $\tilde{\mathcal{E}}$

From the above calculation scheme, we numerically observe that $\alpha$ decreases monotonically with $\mathcal{\tilde{E}}$ . Consider the corresponding change in the morphology of the bubble region, which is described by (4.8). For $\alpha > l^2$ , the thickness profile has a local maximum at $r=0$ given by

(4.33) \begin{equation}\kern2pt\tilde{\kern-2pt t}(0)=\frac{\alpha}{48}.\end{equation}

For $\alpha < l^2$ , however, the thickness profile has a local minimum at $r=0$ , again given by (4.33), and an interior local maximum at $r=(l^2-\alpha)^{1/2}/2^{1/2}$ given by $(l^2+\alpha)/128$ .

Thus, as $\tilde{\mathcal{E}}$ is decreased, the bubble region transitions from a dimple-like morphology where the thickness is maximal at the origin, to an annular morphology where the thickness is maximal at a positive radius smaller than the neck radius; following this transition, the latter maximal thickness continues to increase while the minimal thickness at the origin continues to decrease. This process is schematically shown in Figure 5. Note that the variation with $\mathcal{E}$ of the thickness at the origin follows an opposite trend relative to that predicted in the first distinguished limit (Section 3.3).

In what follows, we look to the neck region to extract the asymptotic dependence of $\alpha$ upon $\tilde{\mathcal{E}}$ . This will allow us to identify the breakdown of the second distinguished limit as the thickness at the origin becomes so small that the approximations underlying the description of the bubble region collapse.

4.2.2 Limit neck problem

Consider the small $\tilde{\mathcal{E}}$ limit of the neck problem formulated in Section 4.1. From (4.24), we have that

(4.34) \begin{equation}\theta \sim \frac{1}{12\sqrt{6}}\tilde{\mathcal{E}}^{-1/6} \quad \text{as} \quad \tilde{\mathcal{E}}\to 0.\end{equation}

We will use this behaviour together with relation (4.24b) to extract the corresponding behaviour of $\alpha$ .

Note that the relation (4.24b) involves the function $\tau(\theta)$ , which is governed by the transformed neck problem (4.19)–(4.23). Thus, an intermediate step is to study that problem as $\theta\to\infty$ in order to extract the behaviour of $\tau(\theta)$ in that limit. As a first step, a dominant-balance argument hints to the leading-order expansions

(4.35a–c) \begin{equation}T'(R')\sim \theta^2T''(R''), \quad Q'(R')\sim \theta^5 Q''(R''), \quad P'(R') \sim P''(R''),\end{equation}

as $\theta\to\infty$ , where we define the double-primed strained coordinate

(4.36) \begin{equation}R''=R'/\theta.\end{equation}

Using these scalings, we see from the flux-conservation equation (4.19a) that vapour production is negligible at leading order, whereby

(4.37) \begin{equation}Q''(R'')\equiv \chi,\end{equation}

where $\chi$ is a constant; from (4.22) and (4.23), the functions $\tau$ and $\nu$ are then seen to possess the leading-order behaviours

(4.38a,b) \begin{equation}\tau \sim \chi \theta^5, \quad \nu \sim \chi \theta^5,\end{equation}

as $\theta\to\infty$ . With Q $^{\prime\prime}$ being constant, the neck problem reduces to the differential equation

(4.39) \begin{equation}\frac{T''^3}{12}\frac{d^3T''}{d{R''}^3}=\chi,\end{equation}

which is obtained by combining (4.19a) and (4.19b), together with the far-field conditions

(4.40a,b) \begin{equation}T''\sim -{R''} \quad \text{as} \quad R''\to - \infty, \quad T''= {R''}^2 + O(1) \quad \text{as} \quad R''\to\infty,\end{equation}

which follow from (4.22) and (4.23), respectively. The differential equation (4.39) is third order, while the far-field expansions (4.40) effectively provide four auxiliary conditions. The extra condition serves to determine the constant $\chi$ .

The constant-flux neck problem consisting of (4.39) and (4.40) is essentially equivalent to that considered in Appendix D of [Reference Jones and Wilson21]. Adopting the numerical result given there to our notation, we have

(4.41) \begin{equation}\chi \doteq 0.051.\end{equation}

In Appendix B, we formulate a slightly generalised constant-flux neck problem which can be reduced to the present one via straightforward rescaling. It will be convenient to refer to this generalised problem later in the paper.

Combining (4.38a) and (4.41), we have a leading-order approximation for $\tau(\theta)$ as $\theta\to\infty$ . It turns out, however, that in order to extract from (4.24b) the small $\tilde{\mathcal{E}}$ behaviour of $\alpha$ , we actually need $\tau$ to higher order. To that end, we substitute the far-field behaviour (4.40a) into the leading nontrivial balance of the primed flux-conservation equation (4.19a), written in terms of the double-primed coordinate and thickness profile; comparison with the primed far-field condition (4.22), which defines $\tau$ , then gives

(4.42) \begin{equation}\tau \sim \chi \theta^5+\frac{\ln \theta-\ln k}{\theta} \quad \text{as} \quad \theta\to\infty,\end{equation}

where k is a pure constant whose actual value is not important for our purposes. (k can be determined by continuing the analysis of the neck problem to higher orders as $\theta\to\infty$ .)

Using (4.34) and (4.42), we find from (4.24b) the behaviour

(4.43) \begin{equation}\ln \frac{1}{\alpha} = \frac{\chi}{(12\sqrt{6})^6}\frac{1}{\tilde{\mathcal{E}}}+\ln\frac{\mathcal{B}^2}{8} + \ln\frac{1}{k} +o(1) \quad \text{as} \quad \tilde{\mathcal{E}}\to0.\end{equation}

It is this behaviour of the logarithm of $\alpha$ that will be needed later, rather than the behaviour of $\alpha$ itself. For the sake of completeness, we also give the latter as

(4.44) \begin{equation}\alpha \sim \frac{8k}{\mathcal{B}^2}\exp\left\{{-\frac{\chi}{(12\sqrt{6})^6}\frac{1}{\tilde{\mathcal{E}}}}\right\} \quad \text{as} \quad \tilde{\mathcal{E}}\to0,\end{equation}

showing more explicitly that $\alpha$ vanishes exponentially as $\tilde{\mathcal{E}}\to0$ . Note that terms up to order unity in (4.43) are crucial to determining the leading-order behaviour (4.44).

4.2.3 Formation of sub-film region and breakdown of second distinguished limit

From (4.33), the rescaled bubble thickness at the origin vanishes linearly as $\alpha\to0$ . Accordingly, the quadratic behaviour of the rescaled thickness profile $\kern2pt\tilde{\kern-2pt t}(r)$ (cf. (4.8)) hints to the formation of a new sub-film region near the origin where $r\simeq \alpha^{1/2}$ and $\kern2pt\tilde{\kern-2pt t}\simeq\alpha$ . Furthermore, (4.12) and (4.43) together imply that $\tilde{q}\simeq 1/\alpha^{1/2}$ in that region. Using (4.5), we see that in terms of the originaldimensionless variables introduced in Section 2.3, the sub-film radial extent, thickness and flux density, respectively, possess the scalings

(4.45a–c) \begin{equation}L_S=\alpha^{1/2}, \quad T_S=\mathcal{B}^2\alpha, \quad Q_S=\mathcal{B}^{10}\alpha^{-1/2}.\end{equation}

The above scalings are based on the asymptotic model derived in the second distinguished limit. Returning to the original pressure-flux relation (2.15b), however, we find using (4.45) that the pressure variation across the sub-film region scales like

(4.46) \begin{equation}\Delta P_S=\frac{\mathcal{B}^4}{\alpha^3}.\end{equation}

There is therefore the possibility that for $\alpha$ sufficiently small pressure variations may no longer be small compared to the constant pressure correction $\mathcal{B}^2\tilde{p}$ appearing in the pressure expansion (4.5c). Indeed, it follows from (4.9) and (4.31) that $\tilde{p}\simeq 1$ as $\alpha\to0$ . Accordingly, the second distinguished limit breaks down for

(4.47) \begin{equation}\alpha \simeq \mathcal{B}^{2/3}.\end{equation}

With (4.44), condition (4.47) implies that the second distinguished limit breaks down when $\tilde{\mathcal{E}}$ is just logarithmically small, i.e., of order $1/\ln \mathcal{B}$ . Continuing to group terms that are of the same algebraic order in $\mathcal{B}$ , this breakdown is seen to occur at the same algebraic scaling (3.28) of $\mathcal{E}$ with $\mathcal{B}$ which was used to define the second distinguished limit in the first place. Nonetheless, the exponential behaviour of $\alpha$ as $\tilde{\mathcal{E}}\to0$ , as expressed by (4.44), provides us with the magnifying glass required to identify the breakdown of the second distinguished limit and next to continue our analysis to even smaller $\mathcal{E}$ , beyond the domain of validity of the approximation obtained in that limit.

4.3 Third distinguished limit

Guided by the above discussion, we shall investigate a third distinguished limit where the algebraic scaling (3.28) of $\mathcal{E}$ with $\mathcal{B}$ still holds despite $\alpha$ being small as in (4.47); as shown schematically in Figure 5, we expect a sub-film region near the origin which is asymptotically distinct from the bubble region. With (3.28) still valid, the asymptotic analyses in Section 4.1 of the bubble, neck, transition and hydrostatic regions need to be only slightly modified. Below, we outline these modifications and subsequently analyse the new sub-film region.

4.3.2 Bubble region

Consider next the bubble region, with the fields expanded as in Section 4.1. The leading-order film thickness $\kern2pt\tilde{\kern-2pt t}$ is still governed by the differential equation (4.6) and boundary condition (4.7b). In contrast, the symmetry boundary condition (4.7a) no longer holds; instead, matching with the sub-film region implies

(4.48) \begin{equation}\kern2pt\tilde{\kern-2pt t}=0 \quad \text{at} \quad r=0,\end{equation}

since the radial extent and thickness of the sub-film region are both small compared to the corresponding properties of the bubble region. Solving (4.6) together with (4.7b) and (4.48) yields

(4.49) \begin{equation}\kern2pt\tilde{\kern-2pt t}=\frac{1}{32}\left(l^2-r^2\right) r^2,\end{equation}

which as expected is the same as (4.8) for $\alpha=0$ . It follows that expression (4.9) for the constant pressure correction $\tilde{p}$ reduces to

(4.50) \begin{equation}\tilde{p} = \hat{h}_{\infty} - \frac{1}{12}.\end{equation}

Following the analysis of the bubble region in Section 4.1, we next consider the leading-order flux density $\tilde{q}$ . Thus, integrating (4.10), with $\kern2pt\tilde{\kern-2pt t}$ given by (4.49), we find

(4.51) \begin{equation}\tilde{q} = \frac{24\mathcal{\tilde{E}}}{r}\left(\ln \frac{l^2r^2}{l^2-r^2}+\zeta\right),\end{equation}

where $\zeta$ is a constant to be determined. (Note that the symmetry condition (4.11) no longer applies.) For later reference, we note the behaviours

(4.52a,b) \begin{equation}\tilde{q}\sim 12\sqrt{6}\tilde{\mathcal{E}}\left\{\ln \frac{l}{3\left(l-r\right)}+\zeta\right\} \; \text{as} \quad r\to l, \quad\tilde{q} \sim \frac{24\tilde{\mathcal{E}}}{r}\left(2\ln r+\zeta\right) \quad \text{as} \quad r\to0,\end{equation}

where for both expansions the error is algebraic in the indicated limit.

Comparing with the analysis in Section 4.1, we see that the constant $\zeta$ effectively replaces $\alpha$ as the single unknown appearing in the leading-order description of the bubble region. Analogously with that analysis, we expect that $\zeta$ could be determined via matching with the neck region; the value of $\zeta$ should then influence the solution in the sub-film region via matching. Thus, in the present regime, information propagates in the inward radial direction.

4.3.3 Neck region

The leading-order problem in the neck region is essentially the same as in Section 4.1, only that the two matching conditions (4.24) need to be updated based on revised matching of flux between the bubble and neck regions. First, instead of (4.24a), we find

(4.53) \begin{equation} \theta = \frac{1}{12\sqrt{6}\tilde{\mathcal{E}}^{1/6}}, \end{equation}

which is just (4.24a) for $\alpha=0$ . Unlike (4.24a), (4.53) gives $\theta$ directly in terms of $\tilde{\mathcal{E}}$ , thence also the function $\tau(\theta)$ defined by solution of the transformed neck problem (cf. (4.19)–(4.23)). Second, instead of (4.24b), we find

(4.54) \begin{equation}\zeta = \theta \tau - \ln \frac{8\theta}{\mathcal{B}^2},\end{equation}

which determines $\zeta$ in terms of $\mathcal{\tilde{E}}$ and $\ln \mathcal{B}$ .

4.3.4 Sub-film region

Consider now the sub-film region near the origin. The scalings (4.45) and (4.46), which characterise the sub-film region, can be expressed solely in terms of $\mathcal{B}$ by substituting the condition (4.47), which defines the third distinguished limit. This suggests defining the strained coordinate

(4.55) \begin{equation}\check{r} = r/\mathcal{B}^{1/3}\end{equation}

and posing the expansions

(4.56a–c) \begin{equation} t= \mathcal{B}^{8/3}\check{t}+ \cdots, \quad q=\mathcal{B}^{29/3}\check{q}+ \cdots,\quad p =2+ \mathcal{B}^{2}\check{p} + \cdots.\end{equation}

The sub-film problem then consists of the differential equations (cf. (2.15) and (2.17))

(4.57a–c) \begin{equation}\frac{1}{\check{r}}\frac{d}{d \check{r}}\left(\check{r}\check{q}\right)=\frac{\mathcal{\tilde{E}}}{\check{t}}, \quad \check{q}=-\frac{\check{t}^3}{12}\frac{d\check{p}}{d \check{r}}, \quad \check{p}= \hat{h}_{\infty}-\frac{1}{\check{r}}\frac{d}{d\check{r}}\left(\check{r}\frac{d \check{t}}{d\check{r}}\right),\end{equation}

the symmetry boundary conditions (cf. (2.18))

(4.58a,b) \begin{equation}\frac{d\check{t}}{d\check{r}}=0, \quad\check{q}=0 \quad \text{at} \quad \check{r}=0,\end{equation}

and the far-field conditions

(4.59a,b) \begin{equation}\check{t} \sim \frac{1}{48}\check{r}^2,\quad\check{q} \sim \frac{48\mathcal{\tilde{E}}}{\check{r}} \ln \check{r} + \frac{24\tilde{\mathcal{E}}}{\check{r}}\left(l^2\ln \mathcal{B}+\zeta\right) \quad \text{as} \quad \check{r}\to\infty,\end{equation}

which follow from matching the film thickness and flux density between the sub-film and bubble regions. Note that the leading logarithmic term in (4.59b) is coupled to the leading quadratic term in (4.59a). Thus, it is the constant term in (4.59b) which serves as the fourth auxiliary boundary condition that closes the above fourth-order problem. In line with our earlier comment regarding information propagating radially inwards, the sub-film problem does not involve any undetermined constants.

In Figure 3, the film thickness at the origin, as predicted by the above sub-film problem, is depicted by the dash-dotted line. (As explained below, the sub-film problem can be solved using the scheme devised for solving the film problem formulated in Section 3.1.) The agreement with the numerical solution of the unified small-slope model is noticeably not as good as in the first and second distinguished limits, though this is to be expected (for the moderately small value of $\mathcal{B}$ used in the numerical simulations) given that the scalings characterising the second and third limits are very close. Nonetheless, the asymptotic approximation obtained in the third distinguished limit qualitatively agrees with the numerical data.

4.4 Second neck–bubble pair and breakdown of third distinguished limit

The sub-film problem is mathematically similar to the film problem found in the first distinguished limit (see Section 3.1), despite these regions being characterised by markedly differing scales and the matching here being with the bubble region rather than directly with the hydrostatic region. In fact, we show in Appendix C that, despite the different form of far-field conditions, the sub-film problem can be exactly mapped to the film problem in Section 3.1. The mapping is nontrivial, however, in the sense that it is defined by stretching factors that are governed by transcendental equations. Using that mapping, we show in Appendix C that the behaviour of the solutions to the sub-film problem as $\mathcal{\tilde{E}}\to0$ is analogous to the behaviour discussed in Section 3.4 of the solutions to the film problem as $\mathcal{E}\to0$ .

The mathematical analogy with the film problem in the first distinguished limit implies that, as $\tilde{\mathcal{E}}$ is decreased, $\check{t}(\check{r})$ transitions from having a local maximum to a local minimum at the origin. Moreover, $\check{t}(0)$ vanishes as $\tilde{\mathcal{E}}\to0$ , leading to the creation of a neck–bubble pair near the origin. Combining the mapping developed in Appendix C, the scalings characterising the the sub-film and film regions and the discussion in Section 3.4 of the creation of the original neck–bubble pair prior to the breakdown of the first distinguished limit, we find that the new bubble region is characterised by the scalings

(4.60a–c) \begin{equation}L_{B,2}=\mathcal{B}^{1/3}\tilde{\mathcal{E}}^{-5/6}, \quad T_{B,2}=\mathcal{B}^{8/3}\tilde{\mathcal{E}}^{-2/3}, \quad Q_{B,2}=\mathcal{B}^{29/3}\tilde{\mathcal{E}}^{5/6}\end{equation}

while the new neck region is characterised by the scalings

(4.61a–c) \begin{equation}L_{N,2}=\mathcal{B}^{1/3}\tilde{\mathcal{E}}^{1/6}, \quad T_{N,2}=\mathcal{B}^{8/3}\tilde{\mathcal{E}}^{1/3}, \quad Q_{N,2}=\mathcal{B}^{29/3}\tilde{\mathcal{E}}^{5/6},\end{equation}

where as usual L, T and Q refer to radial extent, thickness and flux density. These scalings show that the new bubble region is asymptotically smaller than the original one, both in thickness and radial extent (as long as $\tilde{\mathcal{E}} \gg \mathcal{B}^{2/5}$ , but this condition turns out to be irrelevant given the breakdown discussed below).

Analogously to the breakdown of the first distinguished limit, we find that variations in gravitational hydrostatic pressure in the unified small-slope model, negligible in the sub-film region of the third distinguished limit, become non-negligible in the new bubble region for sufficiently small $\tilde{\mathcal{E}}$ . With reference to the interfacial dynamical condition (2.17), this occurs when the order $T_{B,2}/L_{B,2}^2$ curvature of the thickness profile becomes as small as the order $\mathcal{B}^2 L_{B,2}^2$ hydrostatic pressure associated with the quadratic variation in the height of the sphere. It follows that the third distinguished limit breaks down for

(4.62) \begin{equation}\tilde{\mathcal{E}}\simeq \mathcal{B}^{1/4} \Rightarrow \mathcal{E}\simeq \mathcal{B}^{12+1/4}.\end{equation}

Thus, there is a fourth distinguished limit where $\mathcal{E}$ scales as in (4.62). Substituting (4.62) into the scalings (4.60) and (4.61), and referring back to the the unified small-slope model, we see that in this new distinguished limit vapour production is non-negligible only in the new neck–bubble pair. See Figure 5 for a schematic depiction of the fourth distinguished limit.

5.1 Chains of neck–bubble pairs

Consider the scenario hypothesised above, where the asymptotic structure of the thin film includes n pairs of neck and bubble regions. We use the index $m=1,2,\ldots,n$ to identify the pairs, starting from the outermost pair adjacent to the transition region and ending at the innermost one closest to the origin. Let $L_{N,m}$ , $T_{N,m}$ and $Q_{N,m}$ , respectively, be the scalings of the radial extent, thickness and flux density for the mth neck, and $L_{B,m}$ , $T_{B,m}$ and $Q_{B,m}$ the corresponding scalings for the mth bubble. In determining these scalings, we shall specifically assume a scenario analogous to the second (and fourth) distinguished limits, wherein hydrostatic-pressure variations are important in all of the bubble regions; we shall accordingly be skipping the corresponding analogues of the first (and third) distinguished limits. This assumption effectively slaves $\mathcal{E}$ to $\mathcal{B}$ , whereby the requisite scalings will be determined solely in terms of $\mathcal{B}$ .

Based on what we have seen in the first four distinguished limits, we anticipate that:

1. (1) The thickness and radial extent of the bubbles attenuate inwards.

2. (2) Successive bubbles are separated by a relatively narrow neck, i.e., $L_{N,m} \ll L_{B,m}$ . It follows that

   (5.1) \begin{equation}\text{radial position of } m\text{th neck } \simeq L_{B,m}.\end{equation}

3. (3) The radial extent of the first bubble region remains the same when there are multiple neck–bubble pairs:

   (5.2) \begin{equation}L_{B,1}=1.\end{equation}
   In particular, this was the case for the secondary neck–bubble pair that emerged in the third and fourth distinguished limits (the latter corresponding to ‘ $n=2$ ’ here).

   Alternatively, we could assume that the pressure is order unity across all neck and bubble regions. Together with the first and last of this list, the scaling (5.2) would in that case follow from a leading-order balance of the force constraint (2.20).

4. (4) The scalings characterising the transition and hydrostatic regions remain the same as in the previous distinguished limits.

In the neck regions, the pressure variation associated with the capillarity term in (2.17) (i.e., the curvature of the thickness profile) drives a flux in accordance with the flux-pressure relation (2.15b). Thus, if $\Delta P_{N,m}$ denotes the scaling of the pressure variation across that neck, we find from (2.17) and (2.15b) that

(5.3a,b) \begin{equation}\Delta P_{N,m}= \frac{T_{N,m}}{L_{N,m}^2}, \quad Q_{N,m}=T_{N,m}^3\frac{\Delta P_{N,m}}{L_{N,m}}.\end{equation}

These two relations can be combined into

(5.4) \begin{equation}T_{N,m}^4 = Q_{N,m}L_{N,m}^3.\end{equation}

In the bubble regions, the capillarity term and the quadratic hydrostatic term in (2.17) are expected to be comparable (analogously to the case in the second distinguished limit), i.e.,

(5.5) \begin{equation}T_{B,m} = \mathcal{B}^2 L_{B,m}^4.\end{equation}

To relate the scalings that describe successive regions, we derive further scaling relations based on asymptotic matching. We verified in Section 4.1 that, for the first neck region, the thickness profile increases quadratically to the right (cf. (4.23a)) and linearly to the left (cf. (4.22a)). For a chain of neck–bubble pairs where the bubbles attenuate inwards, we expect this to be the case for all the necks. (This can be understood a posteriori from the discussion below regarding the pressure field.) Accordingly, the leading-order thickness profile for each bubble region vanishes quadratically towards the origin and linearly towards its outward bounding neck. (The quadratic vanishing at the origin holds both for bubbles $m=1,\ldots,n-1$ , on which scale the position of the inward bounding neck is indistinguishable from the origin, and the final nth bubble, where axial symmetry applies at the origin.) Matching therefore implies the relations

(5.6a,b) \begin{equation} \frac{T_{B,m}}{L_{B,m}^2} = \frac{T_{N,m+1}}{L_{N,m+1}^2},\qquad \frac{T_{B,m}}{L_{B,m}} = \frac{T_{N,m}}{L_{N,m}}.\end{equation}

We next relate the scalings of the flux density between the different regions. In regions where the vapour production term in (cf. 2.15a) can be neglected, continuity dictates that q attenuates like $1/r$ . In that case, (5.1) and $L_{B,m+1}\ll L_{B,m}$ imply the relations

(5.7a,b) \begin{equation}Q_{B,m} L_{B,m} = Q_{B,m-1} L_{B,m-1}, \quad Q_{B,m}=Q_{N,m}.\end{equation}

A more detailed scaling argument shows that (5.7) also hold for regions where vapour production is non-negligible.Footnote ¹

Solving the recurrence relations (5.4)–(5.7) together with the initial condition $L_{B,1}=1$ ,we find

(5.8a–c) \begin{align}L_{B,m} = \mathcal{B}^{s_m}, \quad T_{B,m} = \mathcal{B}^{2 + 4 s_m},\quad Q_{B,m}=\mathcal{B}^{10-s_m},\quad \end{align}

(5.8d–f) \begin{align}L_{N,m} = \mathcal{B}^{2 - 13 s_m},\quad T_{N,m} = \mathcal{B}^{4 - 10 s_m}, \quad Q_{N,m} = \mathcal{B}^{10-s_m}, \end{align}

wherein

(5.9) \begin{equation}s_m=\frac{8}{7}\left(\frac{1}{8}-\frac{1}{8^m}\right), \quad m=1,2,\ldots,n,\end{equation}

is a monotonically increasing sequence that approaches $1/7$ as $m\to\infty$ (assuming n can be arbitrarily large, see Section 5.3). Note that the scalings (5.8) are consistent with those obtained in Section 4 for the single neck–bubble pair in the second distinguished limit and the two neck–bubble pairs in the fourth distinguished limit.

We can check that the scalings (5.8) are consistent with the assumptions listed at the beginning of this subsection. This is evident for the first two assumptions. Next, substituting (5.8) into the dynamic condition (2.17) shows that $p\simeq 1$ across all neck–bubble pairs, therefore our third assumption is also consistent. In fact, the same argument suggests that $p\sim 2$ in all of those regions, excluding the first neck region across which the pressure at order unity must drop to zero; thus, the radial extent of the first bubble remains $\sim l$ (cf. (5.1)). The fourth assumption is evident given that the scalings of the first neck–bubble pair remain as before.

In those regions where $p\sim 2$ , the leading-order film thickness and flux density are influenced by higher-order terms in an expansion for p. This motivates us to introduce a reduced pressure $p'=p-2-\mathcal{B}^2 h_{\infty}$ , which the dynamic condition (2.17) gives as

(5.10) \begin{equation}p'=-\frac{\mathcal{B}^2}{2}r^2+\mathcal{B}^2t-\frac{1}{r}\frac{d}{dr}\left(r\frac{dt}{dr}\right).\end{equation}

The reduced pressure includes only non-uniform contributions. In particular, the scalings (5.8) show that the second term in (5.10) is negligible compared to the first (quadratic-hydrostatic) and last (capillarity) contributions for all neck and bubble regions. Furthermore, let $P^{\prime}_{N,m}$ and $P^{\prime}_{B,m}$ be the scalings of the reduced pressure in the mth neck and bubble, respectively. For the neck regions, the capillarity contribution in (5.10) dominates the reduced pressure, hence (5.8) gives

(5.11) \begin{equation}P^{\prime}_{N,m} =\mathcal{B}^{16s_m}.\end{equation}

From (5.3), this is also the scaling $\Delta P_{N,m}$ of the pressure variation across the mth neck. For the bubble regions, the by-design balance between the capillarity and quadratic-hydrostatic contributions in the dynamic condition (2.17) implies that these terms together dominate the reduced pressure; it therefore follows from (5.8) and (5.10) that $P^{\prime}_{B,m}=\mathcal{B}^2L_{B,m}^2$ . The flux-pressure equation (2.15b), however, suggests a generally different scaling $\Delta P_{B,n}=Q_{B,m}L_{B,m}/T_{B,m}^3$ for the pressure variation. Thus, using (5.8), we find that

(5.12) \begin{equation}P^{\prime}_{B,m}=\mathcal{B}^{2+2s_m}\quad\gg\quad\Delta P_{B,m}=\mathcal{B}^{4-12s_m}.\end{equation}

We see that the reduced pressure in the bubble regions is uniform, to leading order, consistently with the second distinguished limit (cf. (4.5c)). From these comments, we can deduce that, in the bubble regions, the pressure possesses the expansions

(5.13a) \begin{align}p= 2 + \mathcal{B}^2\tilde{p}_1+\cdots, \quad \text{for} \quad m=1,\end{align}

(5.13b) \begin{align}p= 2 + \mathcal{B}^2\hat{h}_{\infty}+\mathcal{B}^{2+2s_m}\tilde{p}_m +\cdots, \quad \text{for} \quad m=2,3,\ldots,n, \\ \nonumber \end{align}

wherein $\tilde{p}_m$ are constants; similarly, for the neck regions,

(5.14a) \begin{align}p= P_1 +\cdots, \quad \text{for} \quad m=1,\end{align}

(5.14b) \begin{align}p= 2+\mathcal{B}^2P_2 +\cdots, \quad \text{for} \quad m=2,\end{align}

(5.14c) \begin{align}p= 2+\mathcal{B}^2\hat{h}_{\infty}+\mathcal{B}^{16s_m}P_m +\cdots, \quad \text{for} \quad m=3,4,\ldots,n, \\ \nonumber\end{align}

wherein $P_m$ are functions of the corresponding neck coordinates $(r-L_{B,m}l_m)/L_{N,m}$ , $L_{B,m}l_m$ being the approximate position of the mth neck (in particular, $l_1=l$ ). In (5.13) and (5.14), we have approximated the sphere deflection $h_{\infty}$ by its leading-order value $\hat{h}_{\infty}$ (3.18); it can be shown that higher-order corrections do not affect the above expansions up to the orders shown.Footnote ²

5.2 Sequence of distinguished limits

In the preceding scaling analysis, the number n of neck–bubble pairs appears to be arbitrary. We now show that n is effectively determined by the relative smallness of $\mathcal{E}$ and $\mathcal{B}$ . To that end, consider the innermost (i.e., nth) bubble region, which is unique in that axial symmetry applies and dictates that the flux density q must vanish at the origin (cf. (2.18b)). For that to be possible, however, the left- and right-hand sides of the flux-conservation equation (2.15a) must balance to avoid a leading-order solution that is singular at the is singular at the origin. With reference to the scalings (5.8), that balance implies a sequence of distinguished limits:

(5.15) \begin{equation}\mathcal{E} \simeq \mathcal{B}^{12 + 2 s_n},\quad n=1,2,\ldots,\end{equation}

respectively associated with a vapour-film morphology exhibiting a chain of n neck–bubble pairs as described in Section 5.1.

As expected, (5.15) reduces to the second distinguished limit (3.28) for $n=1$ , and the fourth distinguished limit (4.62) for $n=2$ . Note that substitution of (5.8) and (5.15) into the flux-conservation equation (2.15a) shows that the vapour-production term therein is negligible for all neck–bubble pairs but the innermost one, consistently with what we found in the second and fourth distinguished limits. A schematic of the distinguished limit corresponding to the nth term in the sequence (5.15) is depicted in Figure 7(a).

Figure 7.

Depiction of the nth element of the sequence (5.15) of distinguished limits (a) and the distinguished limit (5.17) corresponding to the accumulation point of that sequence (b). The bold letters C, P, E, G and L have the same meaning as in Figure 5.

In order to further corroborate the above picture, in Appendix D we build on the scalings derived in this section to construct a detailed asymptotic description of the vapour film in the distinguished limit corresponding to the nth term in the sequence (5.15). We find that, for each n, the solution in each neck–bubble pair can be recursively determined starting from the outermost pair ( $m=1$ ). Thus, the solution for pair m is independent of that for pair $m+1$ and subsequent pairs. In particular, the flux density in pairs $m<n$ , where vapour production is negligible, can be determined from continuity, in conjunction with the knowledge of the flux density through the outermost neck; we thereby find in Appendix D that the approximation

(5.16) \begin{equation}q \sim \mathcal{B}^{10} \frac{\chi l}{(12\sqrt{6})^5}\frac{1}{r}\end{equation}

holds, for $n>1$ , everywhere except in the innermost pair (cf. (D20)).

5.3 Accumulation of distinguished limits and formation of central region

Formally, taking the limit as $n\to\infty$ in (5.15) gives the scaling

(5.17) \begin{equation}\mathcal{E}\simeq B^{12+2/7}.\end{equation}

In that limit, the scalings (5.8) and (5.11)–(5.12) of the extent, thickness and flux density for the innermost neck–bubble pair give

(5.18a–c) \begin{equation}L_{B,n},L_{N,n} \to \mathcal{B}^{1/7},\quad T_{B,n},T_{N,n} \to \mathcal{B}^{18/7},\quad Q_{B,n},Q_{N,n} \to \mathcal{B}^{69/7},\end{equation}

while the corresponding scalings of the pressure variation and reduced pressure give

(5.19) \begin{equation}\Delta P_{B,n},\Delta P_{N,n},P^{\prime}_{B,n},P^{\prime}_{N,m}\to\mathcal{B}^{16/7}.\end{equation}

The present picture, however, where the vapour film features an increasing number n of neck–bubble pairs, namely gravity capillary waves, as $\mathcal{E}$ is decreased or $\mathcal{B}$ is increased, does not hold for arbitrarily large n. Indeed, large n implies that there are numerous neck–bubble pairs localised near the origin. But for m large like $\ln\ln 1/\mathcal{B}$ , scrutiny of (5.8) shows that the that the corresponding neck and bubble regions are, in fact, characterised by the same algebraic scalings; additionally, neighbouring neck–bubble pairs are characterised by the same algebraic scalings. Thus, for any fixed value of $\mathcal{B}$ , no matter how small, the maximal number of neck–bubble pairs formed as $\mathcal{E}$ is decreased is not infinite, but rather large like $\ln\ln1/\mathcal{B}$ . Incidentally, the latter very mild scaling rationalises why numerically we observe at most a few oscillations of the film’s thickness profile.

We shall refer to the scaling (5.17) as the ‘accumulation point’ of the sequence (5.15) of distinguished limits, since the limits defined by that sequence which are associated with n of order $\ln\ln1/\mathcal{B}$ and larger are effectively inseparable from the accumulation-point limit (5.17). In that limit, we expect the picture of a chain of neck–bubble regions terminating at an innermost bubble to fail sufficiently close to the origin; we instead hypothesise an outer, formally infinite, chain of neck–bubble regions that matches a new central region localised at the origin (as depicted schematically in Figure 7(b)). From (5.18) and (5.19) we deduce the scalings

(5.20a–d) \begin{equation}L_{C} = \mathcal{B}^{1/7},\quad T_{C} = \mathcal{B}^{18/7},\quad Q_{C} = \mathcal{B}^{69/7},\quad P_C^{\prime},\Delta P_{C} = \mathcal{B}^{16/7},\end{equation}

of the radial extent, thickness, flux, reduced pressure and pressure variation, respectively, which characterise that central region.

Recall that the finite chains of neck–bubble pairs of Sections 5.1 and 5.2 were determined via a recursive process in which information propagates inwards. Accordingly, in the accumulation-point limit (5.17) we expect the formally infinite chain of neck–bubble regions to retain the same scalings as derived in Section 5.1 (and asymptotic description as derived in Appendix D) for finite chains, only that now there is no innermost pair and vapour production is negligible for all neck–bubble pairs. Indeed, substitution of the scalings (5.8) and (5.20) into (2.15a) confirms that vapour production is appreciable only in the central region.

5.4 Central region

We now study the central region anticipated in the accumulation-point limit (5.17). Assuming that limit, and given the scalings (5.20), we define the rescaled evaporation number

(5.21) \begin{equation}\breve{\mathcal{E}}=\mathcal{E}/\mathcal{B}^{12+2/7},\end{equation}

the strained radial coordinate

(5.22) \begin{equation}\breve{r}=r/\mathcal{B}^{1/7}\end{equation}

and pose the expansions

(5.23a–c) \begin{equation}t= \mathcal{B}^{18/7}\breve{t}(\breve{r})+\cdots, \quad q= \mathcal{B}^{69/7}\breve{q}(\breve{r}) +\cdots, \quad p = 2 + \mathcal{B}^2\hat{h}_{\infty} + \mathcal{B}^{16/7}\breve{p}(\breve{r}) +\cdots,\end{equation}

where we have approximated the sphere displacement $h_{\infty}$ by its leading-order value $\hat{h}_{\infty}$ , for the same reasons noted below expansions (5.13) and (5.14).

From (2.15) and (2.17), the leading-order fields satisfy

(5.24a–c) \begin{equation}\frac{1}{\breve{r}}\frac{ d (\breve{r}\breve{q})}{d \breve{r}}=\frac{\breve{\mathcal{E}}}{\breve{t}}, \quad \breve{q}=-\frac{\breve{t}^3}{12}\frac{d\breve{p}}{d \breve{r}}, \quad \breve{p} = -\frac{1}{2}\breve{r}^2 -\frac{1}{\breve{r}}\frac{d}{d\breve{r}}\left(\breve{r}\frac{d\breve{t}}{d\breve{r}}\right),\end{equation}

while (2.18) provide the symmetry conditions

(5.25a,b) \begin{equation}\frac{d\breve{t}}{d\breve{r}}=0, \quad \breve{q} = 0 \quad \text{at} \quad \breve{r}=0.\end{equation}

Note that the governing equations (5.24) involve all the physical mechanisms present in the neck and bubble regions described in Sections 5.1 and 5.2, including vapour production, which there was non-negligible only for the innermost pair. This is perhaps unsurprising given that the central region arises from the the accumulation of indistinguishable neck and bubble regions localised near the origin.

The above formulation of the problem governing the central region is missing far-field conditions as $\breve{r}\to\infty$ , which should follow from asymptotic matching between the central region and the formally infinite chain of neck–bubble pairs. One such condition,

(5.26) \begin{equation}\breve{q}\sim \frac{\chi l}{(12\sqrt{6})^5}\frac{1}{\breve{r}},\end{equation}

readily follows from approximation (5.16) for the flux density, which holds throughout the neck–bubble chain given that the vapour-production term in (2.15a) is only non-negligible in the central region. Given (5.26), for large $\breve{r}$ we find that (5.24b) and (5.24c) reduce to

(5.27a,b) \begin{equation}\frac{\chi l}{(12\sqrt{6})^5}\frac{1}{\breve{r}}\sim-\frac{\breve{t}^3}{12}\frac{d\breve{p}}{d\breve{r}} , \quad \breve{p}\sim-\frac{1}{2}\breve{r}^2-\frac{1}{\breve{r}}\frac{d}{d\breve{r}}\left(\breve{r}\frac{d\breve{t}}{d\breve{r}}\right).\end{equation}

These equations contain precisely the leading-order physics governing the neck–bubble chain, including both the neck and bubble regions. They therefore constitute a ‘composite’ approximation that could, in principle, be used to describe that chain as a whole (as an alternative to the matched asymptotics analysis in Appendix D).

This suggests that the central region and neck–bubble chain asymptotically overlap and can therefore be systematically matched. Traditional matching principles [Reference Hinch19], however, are tailored to matching asymptotically distinguished regions. They are therefore inadequate in the present scenario, where the matching is between a distinguished central region and a region comprised of a formally infinite chain of distinguished neck and bubble regions. Such an unconventional relation between asymptotic regions has been encountered in many other problems involving capillary waves [Reference Benilov and Benilov3–Reference Benilov, Chapman, McLeod, Ockendon and Zubkov5, Reference Bowles7, Reference Cuesta and Velázquez12–Reference Duchemin, Lister and Lange14, Reference Hewitt, Balmforth and De Bruyn18, Reference Jensen20, Reference Jones and Wilson21, Reference van Limbeek, Sobac, Rednikov, Colinet and Snoeijer31, Reference Wilson and Jones32]. To the best of our knowledge, detailed matching has never been systematically attempted in these situations, although ad hoc numerical schemes have been devised [Reference Jones and Wilson21].

In the present problem, we envisage this matching could be carried out by first using (5.27) to describe the large- $\breve{r}$ behaviours of the central-region fields in terms of a formally infinite neck–bubble chain, and subsequently comparing the latter chain and the actual one, whose detailed asymptotic description is derived in Appendix D. We assume that such a comparison would allow extracting sufficient information pertaining to the large- $\breve{r}$ chain in order to formulate a numerical scheme for solving the central-region problem. Given the two boundary conditions (5.25) and the far-field condition (5.26), only one additional independent condition would be needed.

We shall not attempt this matching procedure here, partly because a detailed solution of the central region is of little interest and partly because the development of new asymptotic-matching techniques lies beyond the scope of this paper. Instead, we shall exploit the formulation of the central-region problem to extend our asymptotic description of the vapour-film morphology beyond the distinguished scaling (5.17), i.e., to $\mathcal{E}\ll \mathcal{B}^{12+2/7}$ .

5.5 Formation of a uniform-film region

Consider the behaviour of the solution to the central-region problem for small $\breve{\mathcal{E}}$ . Notwithstanding the latter smallness, there must still be a region containing the origin wherein the vapour-production term in (2.15a) remains non-negligible. Otherwise, the leading-order flux density would be singular at the origin. The vapour produced in this region should, moreover, be compatible with the flux density demanded by (5.26). A dominant-balance argument of (5.24) then shows that, to leading order as $\breve{\mathcal{E}}\to0$ , the dynamic condition (5.24c) degenerates to a balance between pressure and gravity,

(5.28) \begin{equation}\breve{p}\sim -\frac{1}{2}{\breve{r}}^2,\end{equation}

whereas (5.24a) and (5.24b) remain valid in that limit. Solving the resulting leading-order problem in conjunction with the boundary condition (5.25), we find the explicit approximations

(5.29a,b) \begin{align}\nonumber \\[-30pt] \breve{t} \sim \left(6 \breve{\mathcal{E}}\right)^{1/4},\quad \breve{q} \sim \frac{1}{12}\left(6\breve{\mathcal{E}}\right)^{3/4}\breve{r},\end{align}

representing a region containing the origin which consists of a film of uniform thickness, with the flux density growing in proportion to the radius. Furthermore, since vapour production is supposed to be negligible except in this uniform-film region, comparison of (5.29b) with the demanded flux density (5.26) suggests that the latter region terminates abruptly at a value of $\breve{r}$ given by

(5.30) \begin{equation}\breve{r}_f= \frac{1}{(12\sqrt{6})^2} \sqrt{\frac{\chi}{3}}\left(6 \breve{\mathcal{E}}\right)^{-3/8},\end{equation}

in a manner which will be clarified in a more general context in the next section.

Note that $\breve{r}_f$ is asymptotically large as $\breve{\mathcal{E}}\to0$ . Hence, the central region expands in that limit, with its innermost part morphing into a uniform-thickness film. For later reference, we give the scalings characterising this uniform-film region in terms of $\mathcal{B}$ and $\mathcal{E}$ . Thus, using (5.20), (5.21) and (5.28)–(5.30), we find for the radial extent, thickness, flux density and pressure variation the respective scalings

(5.31a–d) \begin{equation} L_U=\mathcal{B}^{19/4}\mathcal{E}^{-3/8}, \; T_U=\mathcal{B}^{-1/2}\mathcal{E}^{1/4}, \; Q_U=\mathcal{B}^{21/4}\mathcal{E}^{3/8},\; \Delta P_{U} = \mathcal{B}^{23/2}\mathcal{E}^{-3/4}.\end{equation}

6.1 Second formally infinite sequence of distinguished limits

In the previous section, we identified the formally infinite sequence (5.15) of distinguished limits and subsequently the accumulation-point limit (5.17). In the latter limit, we found a central region near the origin which matches unconventionally with a formally infinite chain of neck–bubble regions. In the present section, we track the film morphology as $\mathcal{E}$ is further decreased, or $\mathcal{B}$ is further increased.

Our starting point is the analysis in Section 5.5 of the central-region problem as $\breve{\mathcal{E}}\to0$ , which revealed that the inner part of that region expands into a uniform film. Eventually, the radial extent of the uniform film may become comparable to the radial extent of the kth $(k=1,2,\ldots)$ bubble in the neck–bubble chain, which order $L_{B,{k}}$ is given by (5.8a). In that scenario, we expect the larger neck–bubble pairs ( $m<k$ ) to remain intact, the kth bubble region to be reduced in size and modified by the expansion of the uniform-film region, and the smaller pairs ( $m > k$ ) to be ‘squeezed’, in a sense to be described, into localised oscillations sandwiched between the uniform-film region and the modified kth bubble region.

The above discussion implies a second sequence of distinguished limits defined by the condition $L_U = L_{B,{k}}$ . Using (5.8a) and (5.31a), the latter condition explicitly yields

(6.1) \begin{equation}\mathcal{E} \simeq \mathcal{B}^{12 + (2 -8 s_{{k}})/3},\quad k=1,2,\ldots,\end{equation}

where $s_k$ is the sequence given by (5.9). The accumulation point (5.17) of the first sequence of distinguished limits (5.15) is also the accumulation point of the second sequence of distinguished limits (6.1). Here, however, decreasing $\mathcal{E}$ corresponds to the index k decreasing from infinity to unity, rather than the other way around. Analogously to the first sequence of distinguished limits, the limits represented by the terms of the sequence (6.1) are closely spaced; in practice, only the ‘last’ limit corresponding to $k=1$ can be clearly observed in numerical solutions of the unified small-scale model.

In the remainder of this section, we consider the distinguished limit corresponding to the kth term in (6.1). It is convenient for that purpose to define the rescaled evaporation number

(6.2) \begin{equation}\hat{\mathcal{E}}_k =\mathcal{E} /\mathcal{B}^{12 + (2 -8 s_{{k}})/3}.\end{equation}

In that limit, we anticipate that the vapour domain decomposes into the following asymptotic regions (schematically depicted in Figure 8):

1. (1) A uniform-film region containing the origin, similar to the one discussed in Section 5.5. This innermost region is the only region where the vapour-production term in (2.15a) is expected to be non-negligible.

2. (2) $k-1$ neck–bubble pairs whose leading-order description remains the same as pairs $m<n$ described in Section 5 and Appendix D, for which vapour production is negligible, with $n=k$ and the index $m =1,2,\ldots, k-1$ used as therein to count the pairs starting from the largest and outermost one. Note that these bubble regions are large relative to the radial extent of the uniform-film region.

3. (3) The kth neck–bubble pair. The radial extent of the kth bubble region is comparable to that of the uniform-film region; hence, the former region is expected to be substantially modified by the expansion of the uniform-film region.

4. (4) A localised sequence of narrow edge regions, which as we will see is necessary to bridge the uniform-film region and the modified kth bubble region.

5. (5) Transition and hydrostatic regions as in Section 4.1.3.

Figure 8.

Depiction of the kth element of the second sequence (6.1) of distinguished limits. The bold letters C, P, E, G and L have the same meaning as in Figure 5.

In this section, our goal is to describe the above decomposition in more detail, stopping short of a detailed analysis of all regions. Thus, in Section 6.2 we briefly revisit the uniform-film region and in Section 6.3 we qualitatively describe the edge regions alluded to above. In the subsequent Section 7, however, we carry out a detailed analysis of the last limit ( $k=1$ ) and obtain accurate approximations describing the ‘squishing’ of the outermost bubble by the expanding uniform-film region.

6.2 Uniform film

The scalings characterising the uniform-film region in the present limit readily follow from substitution of (6.1) into (5.31):

(6.3a–d) \begin{equation}L_{U}=\mathcal{B}^{s_k}, \quad T_{U}=\mathcal{B}^{(8-2s_k)/3}, \quad Q_{U}=\mathcal{B}^{10-s_k},\quad \Delta P_{U} = \mathcal{B}^{2 + 2 s_k}.\end{equation}

The analysis of this region closely follows that in Section 5.5. Thus, we introduce the strained coordinate

(6.4) \begin{equation}r_k = r/\mathcal{B}^{s_k},\end{equation}

pose the expansions

(6.5a,b) \begin{align}t= \mathcal{B}^{(8-2s_k)/3} \hat{t}(r_k) \cdots, \quad q= \mathcal{B}^{10-s_k}\hat{q}(r_k) + \cdots, \end{align}

(6.5c) \begin{align}p \sim 2 + \mathcal{B}^2\hat{h}_{\infty} +\mathcal{B}^{2+2s_k}\hat{p}(r_k) +\cdots \\ \nonumber \end{align}

and derived, from the leading-order balances of (2.15), (2.17) and (2.18), the approximations

(6.6a–c) \begin{equation}\hat{t} = \left(6\hat{\mathcal{E}}_k \right)^{1/4} , \quad \hat{q} = \frac{1}{12} \left(6 \hat{\mathcal{E}}_k \right)^{3/4} r_k, \quad \hat{p} =-\frac{1}{2}r_k^2.\end{equation}

(These are identical to the solutions (5.28) and (5.29) obtained in Section 5.5.)

As in Section 5.5, we expect the uniform region to be finite. We accordingly write

(6.7) \begin{equation}\text{extent of uniform film } \sim \mathcal{B}^{s_k} l_*,\end{equation}

namely that the uniform-film region ends at $\hat{r}=l_*$ , where $l_*$ depends solely on $\hat{\mathcal{E}}_k$ . For $k>1$ , the outermost neck–bubble pair remains unmodified and hence the approximation (5.16) holds everywhere excluding the uniform-film region, where that result is perturbed by appreciable vapour production. Matching the latter approximation with (6.6b) at the edge of the uniform film gives

(6.8) \begin{equation} l_* = \frac{1}{(12\sqrt{6})^2} \sqrt{\frac{\chi}{3}}\left(6 \hat{\mathcal{E}}_k \right)^{-3/8}\quad \text{for} \quad k > 1,\end{equation}

which is identical to the result (5.30) already obtained in Section 5.5. In the case $k=1$ , in which the uniform-film region does modify the outermost neck–bubble pair, a more detailed analysis is necessary to determine $l_*$ (see Section 7).

Recall Figure 3, where we plot the film thickness at the origin as a function of $\mathcal{E}$ , comparing numerical solutions of the unified small-slope model with asymptotic approximations derived in several distinguished limits. In that figure, the present approximation (6.6a) for the thickness of the uniform-film region is depicted by the solid red line, showing excellent agreement over a wide range of $\mathcal{E}$ .

6.3 Localised capillary waves at the edge of the uniform-film region

Clearly, it is not possible to directly match the uniform-film region and the kth bubble region. This hints to the existence of an edge region localised near the outer rim of the uniform film. Let $L_E$ , $T_E$ , $Q_E$ and $\Delta P_E$ be the scalings of the extent, thickness, flux density and pressure variation in that region, respectively. Matching between the finite-extent uniform-film region and the edge region implies the scalings $T_E = T_U$ and $Q_E = Q_U$ . Moreover, the pressure variation in the uniform film is dominated by the quadratic–hydrostatic gravitational contribution, whereas capillarity is also important in the kth bubble; this hints that both of these physical mechanisms should be present in the edge region. With reference to the uniform-film scalings (6.3) and the dynamic condition (2.17), these conditions yield

(6.9a–d) \begin{equation}L_E=\mathcal{B}^{(2-5s_k)/9}, \quad T_E=\mathcal{B}^{(8-2s_k)/3}, \quad Q_E=\mathcal{B}^{10-s_k},\quad \Delta P_E = \mathcal{B}^{(20+4s_k)/9}.\end{equation}

Based on the above scalings, we introduce a strained coordinate centred at the edge of the uniform film,

\begin{equation*}\check{r}= \left(r - \mathcal{B}^{s_k}l_*\right)/\mathcal{B}^{(2-5s_k)/9}\nonumber\end{equation*}

and pose the expansions

(6.10a,b) \begin{align}t= \mathcal{B}^{(8 - 2s_k)/3}\check{t}(\check{r})+\cdots, \quad q= \mathcal{B}^{10-s_k}\check{q}(\check{r})+\cdots, \end{align}

(6.10c) \begin{align}p =2 + \mathcal{B}^2\hat{h}_{\infty} -\frac{\mathcal{B}^{2 + 2 s_k}}{2}l_*^2 + \mathcal{B}^{(20+4s_k)/9}\check{p}(\check{r})+\cdots. \\ \nonumber \end{align}

The dynamic condition (2.17) at leading order gives

(6.11) \begin{equation}\check{p}=- l_* \check{r} - \frac{d^2 \check{t}}{d\check{r}^2}.\end{equation}

Furthermore, substitution of (6.2) and (6.10a,b) into the flux-conservation equation (2.15a) confirms that $\check{q}$ is uniform. Matching with the uniform-film flux approximation (6.6b) therefore gives

(6.12) \begin{equation}\check{q} = \frac{1}{12}\left(6 \hat{\mathcal{E}}_k\right)^{3/4} l_*.\end{equation}

Using (6.10b)–(6.12), the flux-pressure relation (2.15b) gives the differential equation

(6.13) \begin{equation}\check{t}^3\left( l_*+\frac{d^3 \check{t}}{d\check{r}^3}\right)= \left(6 \hat{\mathcal{E}}_k\right)^{3/4} l_*\end{equation}

governing the leading-order thickness profile. Matching with the uniform-film approximation (6.6a) gives the far-field condition

(6.14) \begin{equation}\check{t} \to \left(6 \hat{\mathcal{E}}_k\right)^{1/4}\quad \text{as}\quad \check{r}\to-\infty.\end{equation}

To close the problem governing $\check{t}$ , two additional far-field conditions are required in the diametric limit $\check{r}\to\infty$ .

The differential equation (2.15) has been studied by several authors, in particular in the context of gravity-capillary waves at the interface of a two-dimensional thin film which flows at constant flux [Reference Cuesta and Velázquez12, Reference Wilson and Jones32]. It has been shown that solutions of that equation which approach a constant value as $\check{r}\to-\infty$ , as required by (6.14), exhibit oscillations of growing amplitude as $\check{r}\to\infty$ . Moreover, the latter oscillations are asymptotically comprised of an infinite chain of pairs of neck and bubble regions. Specifically, in the respective neck regions, (2.15) reduces to

(6.15) \begin{equation}\check{t}^3 \frac{d^3 \check{t}}{d\check{r}^3}\approx \left(6 \hat{\mathcal{E}}_k\right)^{3/4} l_*,\end{equation}

which we saw describes neck regions where vapour production is negligible (cf. (4.39)); furthermore, in the respective bubbles regions, (2.15) reduces to

(6.16) \begin{equation}\frac{d^3 \check{t}}{d\check{r}^3} \approx - l_*,\end{equation}

which, upon integration, gives a radially localised version of the equations we found to describe the unmodified bubble regions (cf. (4.6) and (D5)). It is therefore seen that the chain of neck–bubble pairs associated with the large $\check{r}$ behaviour of (6.13) constitutes a simplified version of the chain of neck–bubble pairs studied in Section 5 and Appendix D.

While the oscillatory growth of $\check{t}$ forbids direct matching between the edge region and the kth-bubble region, the above observations suggest that it may be possible to match the latter bubble region with a formally infinite sequence of edge-localised (i.e., two-dimensional) neck–bubble pairs, which, in turn, could be matched with the edge region, similarly to the unconventional matching described in Section 5.4. We expect that such matching could be employed to close the edge-region problem (6.13) and (6.14).

The scalings of the edge-localised pairs can be obtained in a manner similar to that presented in Section 5. (Accordingly, the solution in each of these regions can be determined from an analysis similar to that carried out in Appendix D). In particular, scalings for the outermost and largest edge-localised pair are shown in Figure 8. We note that the chain of edge-localised pairs in the present limit evolved from pairs $m>k$ that existed in the distinguished limit (5.17). Similarly, in the present limit we expect the kth pair to be continuously deformed from a non-localised to an edge-localised pair as $\mathcal{\hat{E}}_k$ decreases.

7.1 Outermost bubble

The analysis of the outermost bubble region closely follows that of the bubble region in the second distinguished limit (Section 4.1), the main differences being the location and conditions satisfied at the inner boundary, as well as the negligibility of vapour production. (The analysis of the outermost bubble region in Appendix D, where vapour production is negligible, is even closer.) In particular, the pressure, thickness and flux-density fields are expanded as in (4.5), where it is understood that the constant $\tilde{p}$ remains to be determined in the present limit. The leading-order thickness profile $\kern2pt\tilde{\kern-2pt t}$ is governed by the differential equation (4.6). The bubble region is now bounded by the uniform-film edge $r=l_*$ and the outermost neck position $r=l>l_*$ . (Since $p\sim 2$ in all regions inwards from the outermost neck, it follows from a leading-order balance of the force condition (2.20) that l still equals $\sqrt{2/3}$ .) At the inner boundary, indirect matching with the the uniform-film region, as discussed in Section 6.3 for arbitrary k, implies the conditions

(7.1a,b) \begin{equation}\kern2pt\tilde{\kern-2pt t}=0, \quad \frac{d\kern2pt\tilde{\kern-2pt t}}{dr}=0 \quad \text{at} \quad r=l_*.\end{equation}

At the outer boundary, matching with the outermost neck implies the condition

(7.2) \begin{equation}\kern2pt\tilde{\kern-2pt t}=0 \quad \text{at} \quad r=l.\end{equation}

Solving (4.6) together with (7.1), we find

(7.3) \begin{equation}\kern2pt\tilde{\kern-2pt t}=\frac{1}{32}\left(l_*^2-r^2\right)\left(r^2-\eta\right) + \frac{1}{16}{l_*}^2(l_*^2-\eta)\ln\frac{r}{l_*},\end{equation}

where we have introduced the auxiliary parameter

(7.4) \begin{equation}\eta=8(\hat{h}_{\infty}-\tilde{p})-l_*^2,\end{equation}

with $\hat{h}_{\infty}$ given by (4.31). Applying the boundary condition (7.2), we obtain the relation

(7.5) \begin{equation}\eta = l^2 + \left(l^2-l_*^2\right)\frac{2\ln\frac{l}{l_*}}{l^2/l_*^2-1 - 2\ln\frac{l}{l_*}}.\end{equation}

Next, expansion of the flux-conservation equation (2.15a) in the present limit confirms that vapour production is negligible. The leading-order balance of that equation therefore reduces to a statement of continuity, whereby the leading-order flux density $\tilde{q}$ attenuates like $1/r$ . In conjunction with matching with the uniform-film region, we find

(7.6) \begin{equation}\tilde{q}=\frac{l_*\hat{q}(l_*)}{r},\end{equation}

wherein

(7.7) \begin{equation}\hat{q}(l_*) = \frac{l_*}{12}\left(6\hat{\mathcal{E}_1}\right)^{3/4}\end{equation}

is the uniform-film approximation for the flux density (6.6b) evaluated at the edge (6.7).

The solutions (7.3) and (7.6) depend on the radial extent $l_*$ of the uniform film. To determine this quantity we need to go beyond the matching condition (7.2) and investigate the outermost neck in more detail.

7.2 Outermost neck

The analysis of the outermost neck region closely follows that of the neck region in the second distinguished limit (Section 4.1), the main difference being the negligibility of vapour production here. (The analysis of the outermost neck region in Appendix D, where vapour production is negligible, is essentially identical except from the different matching conditions here owing to changes to the outermost bubble region. In particular, we employ the strained coordinate R defined by (4.14) and expand the thickness, flux density and pressure fields as in (4.15).

Expanding the flux-conservation equation (2.15a) in the present limit confirms that vapour production is negligible, as mentioned above. It readily follows that the leading-order flux density Q(R) is constant, thence matching with the outermost bubble region gives

(7.8) \begin{equation}Q(R)= \tilde{q}(l),\end{equation}

where $\tilde{q}(l)$ is obtained by evaluating (7.6) at $r=l$ .

Substitution of (7.8) into a combined leading-order balance of the flux-pressure relation (2.15b) and dynamical condition (2.17) yields the differential equation

(7.9) \begin{equation}\frac{T^3}{12}\frac{d^3T}{dR^3}=\tilde{q}(l),\end{equation}

which governs the leading-order thickness profile. Matching with the outermost bubble region, using (7.3), yields the far-field condition

(7.10) \begin{equation}T\sim -\vartheta R \quad \text{as} \quad R\to-\infty,\end{equation}

wherein

(7.11) \begin{equation}\vartheta =\frac{1}{16l}\left(l^2-{l_*}^2\right)\left(2l^2+l_*^2-\eta\right).\end{equation}

Furthermore, matching with the transition region yields the far-field condition

(7.12) \begin{equation} T \sim R^2 - 2 R_s R \quad \text{as} \quad R\to\infty,\end{equation}

wherein $R_s$ corresponds to a constant shift of the neck radius. That constant was calculated in Appendix A assuming the second distinguished limit, giving the closed-form result (4.32). It is straightforward to revisit that derivation in the present limit to confirm that (4.32) remains valid.

The neck problem consisting of (7.9)–(7.12) is of the form considered in Appendix B. It follows that the constants $\vartheta$ and $\tilde{q}(l)$ are related as

(7.13) \begin{equation}\vartheta=\left\{\tilde{q}(l)/\chi\right\}^{1/5},\end{equation}

wherein $\chi$ is the numerically evaluated constant (4.41). Combining (7.13) with (7.6), (7.7) and (7.11), we find the relation

(7.14) \begin{equation}\left(l^2-l_*^2\right)\left(2l^2+l_*^2-\eta\right)=\frac{16 l_*^{2/5}l}{(12\chi l)^{1/5}}\left(6\hat{\mathcal{E}}_1\right)^{3/20}.\end{equation}

Numerical solution of the pair of equations (7.5) into (7.14) furnishes $l_*$ as a function of $\hat{\mathcal{E}}_1$ . With $l_*$ determined, the leading-order asymptotic approximations describing the outermost bubble and neck regions are closed.

In Figure 10, we plot the thickness profile for $\mathcal{B}=0.1$ and $\mathcal{E}=10^{-29}$ , with the solid line representing a numerical solution of the unified small-slope model and the dashed, dotted and dash-dotted lines representing the leading-order approximations obtained in the present limit for the uniform-film, outermost bubble and outermost neck regions, respectively. The uniform-film, bubble and neck regions are clearly visible in the numerical solution, as are the anticipated secondary oscillations at the edge of the uniform-film region. The particularly good agreement between the numerical solution and the uniform-film approximation suggests that the error in the latter is exponentially small (in contrast to the numerically appreciable algebraic errors in the bubble and neck regions). This can indeed be confirmed by higher order scrutiny of the uniform-film problem and the unconventional matching between that region and the outermost bubble region.

7.3 Small- $\hat{\mathcal{E}}_1$ limit

From (7.5) and (7.14), we find the behaviour

(7.15) \begin{equation}l_* - l \sim-\frac{\sqrt{6}\left(6\hat{\mathcal{E}}_1\right)^{3/40}}{(12\chi l^4)^{1/10}} \quad \text{as} \quad \hat{\mathcal{E}}_1\to0.\end{equation}

This expression shows that the uniform film expands across the entire outermost bubble region. For sufficiently small $\hat{\mathcal{E}}_1$ , the remainder of the outermost bubble region joins with the capillary waves at the edge of the uniform film. The resulting morphology is that of a uniform film bounded by capillary waves, as we observe numerically. We expect a similar morphology for moderately sized spheres, i.e., for $\mathcal{B}\simeq 1$ , in the weak-evaporation limit, though this cannot be confirmed here as the unified small-slope lubrication model does not hold in that regime. In that sense, we expect the ‘last’ limit considered in this section to be not just the last of the second sequence (6.1) of distinguished limits, but also the last limit relevant to the description of the evolution of the vapour film across the small-sphere regime. This is consistent with the fact, already noted in the introduction, that a similar morphology was found in the case of large drops levitated above a heated liquid bath [Reference van Limbeek, Sobac, Rednikov, Colinet and Snoeijer31], although in that extreme scenario the gravity-dominated film region is not asymptotically uniform as it is here.