2.1. Full minimal model and its control parameters

Binary droplet evaporation is a complex multi-physics problem involving interfacial mass transfer, which has been successfully analysed numerically using transient (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ) and quasi-stationary models (Diddens et al. Reference Diddens, Li and Lohse2021). We study the stability of quasi-stationary flow solutions, focusing on solutal Marangoni instabilities when the most volatile component has lower surface tension.

Typically, rapid evaporation in experiments with water–ethanol droplets induces significant transient and thermal effects. Our minimal model isolates solutal Marangoni effects, making it directly applicable only to slow-evaporating systems where thermal effects are negligible. Limitations and implications for real systems are discussed in § 6.

When modelling the evolution of each component $\alpha$ in the liquid phase, it is sufficient to track the mass fraction $y_A$ of the more volatile component $A$ , considering that $y_A + y_B = 1$ , where $y_B$ is the mass fraction of the component $B$ . Typically, the liquid’s velocity $\boldsymbol{u}$ is much greater than the interface velocity $\boldsymbol{u}_I$ (Diddens et al. Reference Diddens, Li and Lohse2021). For slowly evaporating droplets, we assume only minimal compositional variations in the liquid phase during evaporation, making it reasonable to express $y_A$ as $y_A = y_{A,0} + y$ , where $y_{A,0}$ is the spatially averaged composition and $y$ is a small perturbation. Naturally, the averaged mass fraction $y_{A,0}$ changes over time; however, under the specified conditions, the process can be considered quasi-stationary at each instant of the drying time, as shown in § 6.

The composition-dependent liquid properties can then be approximated using a first-order Taylor expansion around $y_{A,0}$ . We neglect the dependence of the liquid’s properties on temperature for simplicity, but we acknowledge these can often be relevant (see § 6), therefore limiting the model to cases where thermal effects are less pronounced. While variations in the liquid density can influence the flow pattern within binary droplets (Edwards et al. Reference Edwards, Atkinson, Cheung, Liang, Fairhurst and Ouali2018; Li et al. Reference Li, Diddens, Lv, Wijshoff, Versluis and Lohse2019; Diddens et al. Reference Diddens, Li and Lohse2021), this effect is neglected in the present study for simplicity, in order to focus instead on pure Marangoni-induced flows. The droplet is assumed to perfectly form a spherical-cap shape throughout the calculations in this section, implying that the capillary (and Bond) number(s) are effectively zero. The contact angle with the substrate is called $\theta$ . All material properties are assumed to be constant regardless of the composition (and temperature), except for the liquid–gas surface tension, given by $\sigma (y_A) =\sigma (y_{A,0})+y\partial _{y_{A}}\sigma$ . We introduce non-dimensionalised scales (marked with tildes) for space, time, velocity and pressure:

(2.1) \begin{align} x = V^{1/3}\tilde {x}, \qquad t = \dfrac {V^{2/3}}{D_0}\tilde {t}, \qquad \boldsymbol{u} = \dfrac {D_0}{V^{1/3}}\tilde {\boldsymbol{u}},\qquad p = \dfrac {\mu _0 D_0}{V^{2/3}}\tilde {p}, \end{align}

where $V$ is the volume of the droplet and $D_0$ and $\mu_0$ are, respectively, the diffusion coefficient and dynamic viscosity evaluated at the average composition $y_{A,0}$ . Similarly, one can obtain the non-dimensionalised gas concentration $\tilde {c}$ , whose profile around the droplet determines the evaporation rate of each component, assuming the droplet evaporates under ambient conditions (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997, Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Popov Reference Popov2005). As shown by Diddens et al. (Reference Diddens, Li and Lohse2021), when the droplet is well mixed, which is achieved when there is sufficient mixing induced by solutal Marangoni flow, the effect of local liquid composition fluctuations on the vapour pressure – and hence on the evaporation rate – can be safely neglected. In this framework, the vapour concentration is assumed to depend solely on the spatially averaged composition of the droplet, meaning that higher-order terms in the Taylor expansion of Raoult’s law about this average are omitted. Although incorporating these first-order corrections could yield a more precise description, doing so would introduce additional parameters and complicate the model without significantly altering the results. This simplification allows us to write the vapour concentration as

(2.2) \begin{equation} \tilde {c} = (c_\alpha - c_\alpha ^\infty )/(c^{{eq}}_\alpha - c_\alpha ^\infty ), \end{equation}

where $c_\alpha ^\infty$ is the ambient vapour concentration far from the droplet and $c^{{eq}}_\alpha = c_\alpha ^{{pure}} \gamma _\alpha x_\alpha$ . Here, $c_\alpha ^{{pure}}$ is the saturation concentration of the pure component $\alpha$ ( $A$ or $B$ ) at its average composition, $\gamma _\alpha$ is the activity coefficient, and $x_\alpha$ is the liquid mole fraction. For simplicity, the saturation concentration $c_\alpha ^{{eq}}$ is assumed to be constant. Following this scaling, $\tilde {c}=1$ at the liquid–gas interface and $\tilde {c}=0$ far from the droplet. As originally validated by Deegan et al. (Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000) and Hu & Larson (Reference Hu and Larson2002) for pure liquids and by Diddens (Reference Diddens2017) for multicomponent droplets, the solution of the Laplace equation

(2.3) \begin{equation}\tilde{\nabla} ^2 \tilde {c} = 0 \end{equation}

in the gas phase gives $\tilde {c}$ . The solution of (2.2) depends only on the geometry of the droplet. We disregard Stefan’s flow (Brutin & Starov Reference Brutin and Starov2018) and natural convection in the gas phase based on the results of Diddens et al. (Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ), who have shown that these effects are irrelevant specifically for evaporating water–ethanol droplets at ambient conditions. The dimensionless evaporation rate is given by the diffusive mass flux

(2.4) \begin{equation} \tilde {j} = - \tilde {\partial }_n \tilde {c}, \end{equation}

with $\partial_n$ being the gradient in the normal direction $\boldsymbol{n}$ of the interface, which induces variations in the liquid composition through the flux boundary condition

(2.5) \begin{equation} -\tilde {\boldsymbol{\nabla }} y \boldsymbol{\cdot } \boldsymbol{n} = {Ev} \tilde {j}. \end{equation}

Here, ${Ev}$ is the evaporation number, which measures the strength of the concentration gradient created in the liquid due to the preferential evaporation of one of the components, defined as (Diddens et al. Reference Diddens, Li and Lohse2021)

(2.6) \begin{align} {Ev} = \dfrac {(1 - y_{A,0})D_A^{\textit{vap}}(c_{A,0}^{{eq}} - c_A^\infty )-y_{A,0} D_B^{\textit{vap}}(c_{B,0}^{{eq}} - c_B^\infty )}{\rho _0 D_0}. \end{align}

where $D^{\textit{vap}}_\alpha$ is the vapour diffusivity of the component $\alpha$ and $\rho_0$ is the mass density evaluated at the average composition $y_{A,0}$ . In the liquid phase, the conservation of mass and momentum is described by the following Stokes equations, respectively, given by

(2.7) \begin{align}&\qquad \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } \tilde {\boldsymbol{u}} = 0, \end{align}

(2.8) \begin{align}& -\tilde {\boldsymbol{\nabla }} \tilde {p} + {\tilde{\boldsymbol{\nabla}}}^2 \tilde{\boldsymbol{u}} = 0, \end{align}

where gravity is neglected for simplicity. For convenience, we define the modified mass fraction variations $\xi = y / {Ev}$ , which allows us to reduce the number of control parameters of the problem. When substituting $y = {Ev} \xi$ in (2.4), the liquid–gas interface boundary condition becomes

(2.9) \begin{equation} -\tilde {\boldsymbol{\nabla }} \xi \boldsymbol{\cdot } \boldsymbol{n} = \tilde {j}. \end{equation}

Since the diffusion coefficient in the liquid phase is typically small, we consider the full convection–diffusion equation for the variations of mass fraction $\xi$ :

(2.10) \begin{equation} \partial _{\tilde {t}} \xi + \tilde {\boldsymbol{u}} \boldsymbol{\cdot } \tilde {\boldsymbol{\nabla }} \xi = \tilde {\nabla }^2 \xi . \end{equation}

Equation (2.8) contains the only relevant transient term that can affect the stability of this system, since the large Schmidt number ( ${Sc} = {\mu _0}/{(\rho _0 D_0)} \sim 10^3$ ) indicates that velocity is ‘enslaved’ to compositional gradients. The interface is subject to Marangoni stresses:

(2.11) \begin{equation} \boldsymbol{n} \boldsymbol{\cdot } (\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}} + (\tilde {\boldsymbol{\nabla }}\tilde {\boldsymbol{u}})^t) \boldsymbol{\cdot } \boldsymbol{t} = {Ma} \tilde {\boldsymbol{\nabla }}_t \xi \boldsymbol{\cdot } \boldsymbol{t}, \end{equation}

where $\tilde {\boldsymbol{\nabla }}_t$ represents the surface gradient operator. The Marangoni number $\textit{Ma}$ is defined as

(2.12) \begin{equation} {Ma} = \dfrac {V^{1/3} \partial _{y_A} \sigma }{D_0 \mu _0} {Ev} , \end{equation}

and $\mu _0$ is the dynamic viscosity evaluated at the average composition $y_{A,0}$ .

We consider a sphericalcap for the droplet shape (i.e. capillary number $\textit{Ca} \rightarrow 0$ ). We restrict our analysis to the case where the droplet is in a quasi-stationary state, such that there is no total mass transfer across the interface. This scenario can be realised theoretically by controlling the ambient humidities of components $A$ and $B$ so that evaporative mass loss of $A$ is compensated by the mass gain of $B$ . In this case, the total mass of the droplet remains constant and thermal gradients are inexistent. We emphasise that this assumption is generally not fully achievable experimentally. We can therefore enforce $\tilde {\boldsymbol{u}} \boldsymbol{\cdot } \boldsymbol{n} = 0$ at the droplet–gas interface by a normal traction. At the substrate, a no-slip boundary condition is imposed, i.e. $\tilde {\boldsymbol{u}} = 0$ , and the Neumann boundary condition $\tilde {\boldsymbol{\nabla }} \xi \boldsymbol{\cdot } \boldsymbol{n} = 0$ is applied to prevent normal mass flux through the substrate. To eliminate the null space of the pressure field, an average pressure constraint is imposed, i.e. $\int \tilde {p} \, \mathrm{d}V = 0$ . This is enforced numerically via a Lagrange multiplier. Similarly, a zero-average constraint is imposed to remove the constant shift for $\xi$ .

The system of equations is solved using a finite element method implemented in the software package pyoomph (Diddens & Rocha Reference Diddens and Rocha2024), which is based on oomph-lib (Heil & Hazel Reference Heil and Hazel2006) and GiNaC (Bauer, Frink & Kreckel Reference Bauer, Frink and Kreckel2002). All domains are discretised using an axisymmetric mesh composed of triangular elements. Linear basis functions are employed for $\xi$ and $\tilde {p}$ , while quadratic basis functions are used for the $\tilde {\boldsymbol{u}}$ field. The solution process begins with the construction of the residual system of equations, followed by the computation of its Jacobian and mass matrices. Quasi-stationary solutions are obtained by providing an initial value for the degrees of freedom and setting the mass matrix contributions to zero. Due to hysteresis, the quasi-stationary solutions can vary depending on the initial values assigned to the degrees of freedom. Additionally, the stability of a solution depends on the sign of the real part of the eigenvalues $\lambda$ of the Jacobian matrix. Stability analysis is conducted using the shift-inverted Arnoldi method to evaluate the eigenvalues. If a bifurcation is detected – where the real part of an eigenvalue crosses zero – the bifurcation curve is tracked using the method outlined by Diddens & Rocha (Reference Diddens and Rocha2024). These bifurcations, classified as either fold or Hopf (in our particular system), lead to distinct types of instabilities. Fold bifurcations result in a change in solution stability, whereas Hopf bifurcations give rise to oscillatory behaviour, distinguishable by the presence of an imaginary component in the eigenvalue. For a more detailed explanation of the numerical methods used, we refer to Diddens & Rocha (Reference Diddens and Rocha2024).

2.2. Instabilities close to $\theta =90^\circ$

In the limit of zero capillary number, there are two control parameters in the problem of this study: $\textit{Ma}$ and the contact angle $\theta$ . For glycerol–water mixtures, where $\textit{Ma}$ is positive, the quasi-stationary solutions for the flow and composition fields within the droplet are stable, unique, axisymmetric, and not subject to hysteresis for all $\theta$ (Diddens et al. Reference Diddens, Li and Lohse2021). Exactly at $\theta = 90^\circ$ , $\tilde {j}$ is uniform along the liquid–gas interface, resulting in a diffusive profile for $\xi$ and absence of flow within the droplet. On hydrophilic substrates, i.e. forming contact angles $\theta \lt 90^\circ$ , the larger evaporation rate at the contact line reduces the surface tension locally, inducing a Marangoni flow from the contact line towards the top of the droplet. On hydrophobic substrates, i.e. forming contact angles $\theta \gt 90^\circ$ , the Marangoni flow is in the opposite direction, since the evaporation rate is larger at the top of the droplet.

In contrast, water–ethanol mixtures, where $\textit{Ma}$ is negative, exhibit much more complex behaviour. Strong transient effects and asymmetric, seemingly chaotic flows have been reported both experimentally (Machrafi et al. Reference Machrafi, Rednikov, Colinet and Dauby2010; Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ) and numerically (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ). Here, we examine the stationary solutions at negative Marangoni numbers $\textit{Ma}$ with an initially small modulus, before the onset of these instabilities, and gradually decrease $\textit{Ma}$ to more negative numbers while monitoring the linear stability of the solutions.

Figure 2. Phase portrait of solution branches as a function of $\textit{Ma}$ , for $\theta = 89^\circ$ (a), $\theta = 90^\circ$ (b) and $\theta = 91^\circ$ (c). The $x$ -axis represents the average tangential velocity $U_t$ at the liquid–gas interface, a key measure of the flow field, indicating the flow direction within the droplet when a single vortex is present. The flow field and $\xi$ profile for selected A $\rightarrow$ R, R $\rightarrow$ A and 2V solutions are also depicted. Here, the diverging blue to red colours in the contour plots represent increasing $\xi$ of the fluid with the lowest surface tension, e.g. ethanol in water–ethanol mixtures. At $\theta = 90^\circ$ , the $\xi$ profile is purely diffusive for low $\textit{Ma}$ , becoming unstable at a critical $\textit{Ma}_{{cr}}$ in an imperfect pitchfork bifurcation (green). The inset in (b) shows the $\theta = 90^\circ$ phase portrait in a range of $\textit{Ma}$ from 150 to 190, where the imperfect pitchfork is clearly visible. All other solution branches lose their stability in either fold (red) or Hopf (blue) bifurcations.

At exactly $\theta = 90^\circ$ , we observe a purely diffusive profile for $\xi$ , but this solution is only stable at sufficiently low $|{Ma}|$ . Here, a non-zero Rayleigh number ( ${Ra}$ ), i.e. buoyant forces, would initiate flow, but we emphasise again that our analysis here is for ${Ra}=0$ , i.e. negligible density contrasts. Above the critical $|{Ma}|$ value, $-{Ma}\gt -{Ma}_{{cr}}\approx 167$ , multiple solution branches coexist (see phase portrait in figure 2 b). Depending on the initial conditions, the flow can develop either a vortex from the apex to the rim (A $\rightarrow$ R), a vortex from the rim to the apex (R $\rightarrow$ A), or two competing vortices (2V). We do not rule out the possibility of additional stable solution branches, such as multiple vortices, but we have not observed them in our axisymmetric numerical simulations. At $\textit{Ma}_{{cr}}$ , the flow field’s symmetry is broken in an imperfect pitchfork bifurcation (see e.g. Golubitsky & Schaeffer Reference Golubitsky and Schaeffer1979; Strogatz Reference Strogatz2018), corresponding to the normal form $\dot {x}=rx+hx^2-x^3$ , where $r$ is the control parameter and $h$ the imperfection parameter. If $h=0$ , we would retrieve a perfect pitchfork bifurcation. In other words, the diffusive profile becomes unstable through a transcritical bifurcation, where it exchanges stability with the R $\rightarrow$ A solution branch, while the A $\rightarrow$ R branch loses stability via a fold bifurcation at a slightly lower $\textit{Ma}$ (see inset in figure 2 b). Our numerical tests suggest that the no-slip boundary condition at the substrate and the axisymmetric coordinate system disrupt the physical symmetry that would otherwise be present, leading to the imperfect pitchfork bifurcation. At higher $|{Ma}|$ , the R $\rightarrow$ A, 2V and A $\rightarrow$ R solution branches each lose stability through either a fold or a Hopf bifurcation.

For slightly lower $\theta$ , e.g. $\theta = 89^\circ$ , the diffusive state no longer exists and the flow exhibits an A $\rightarrow$ R solution for low $|{Ma}|$ . Despite the induced Marangoni flow pushing the fluid from the top towards the contact line (i.e. in the A $\rightarrow$ R direction), an R $\rightarrow$ A solution can still be found, which shows strong hysteresis. The stability of these solutions is limited by fold bifurcations at different $\textit{Ma}$ (see figure 2 a). For slightly larger $\theta$ , e.g. $\theta = 91^\circ$ , the flow exhibits an R $\rightarrow$ A solution for low $|{Ma}|$ . Curiously, this solution branch loses its stability in a fold bifurcation at a lower $|{Ma}|$ than its A $\rightarrow$ R counterpart for the same set of parameters (see figure 2 c), which can again be attributed to the three-dimensional geometry and the presence of the no-slip condition. While Marangoni flow pushes the fluid from the contact line towards the top of the droplet (i.e. in the R $\rightarrow$ A direction), an A $\rightarrow$ R solution presents, counterintuitively, for large enough $\textit{Ma}$ , better stability than the R $\rightarrow$ A one.

This intriguing introduction to the vast array of bistable solutions within the narrow $\textit{Ma}$ – $\theta$ phase space investigated here sets the stage for the topic to be discussed in the following subsection, which covers the stability of the quasi-stationary solutions on an extensive $\textit{Ma}$ – $\theta$ phase diagram.

2.3. Ma vs $\theta$ phase diagram

In this section, we utilise the bifurcation tracking tool outlined by Diddens & Rocha (Reference Diddens and Rocha2024). Building upon the initial findings from the preceding subsection, we employ arclength continuation on $\theta$ to trace the bifurcation curves along the $\textit{Ma}$ – $\theta$ phase diagram. Arclength continuation on $\theta$ requires remeshing the domain whenever the grid deformation surpasses a certain threshold. We used highly refined meshes near the liquid–gas interface, with a total of $\sim$ 40 000 degrees of freedom, to accurately capture the Marangoni flow. The curves in figure 3 remained unchanged with further grid refinement, confirming mesh independence.

Figure 3. Bifurcation diagram of the quasi-stationary axisymmetric solutions as a function of $\textit{Ma}$ and $\theta$ . The diagram is divided into three regimes, where the stable solutions are either an A $\rightarrow$ R solution (blue colour), an R $\rightarrow$ A solution (yellow colour) or a 2V solution (plum colour). The transition between these regimes is marked by fold (red) or Hopf (blue) bifurcations. Above the upper Hopf bifurcation curves, no stable quasi-stationary solutions are found (grey).

We find a complex bifurcation diagram where multiple solution branches coexist across a broad spectrum of $\textit{Ma}$ and $\theta$ values (figure 3). This diagram is segmented into three distinct regimes, each representing stable quasi-stationary axisymmetric A $\rightarrow$ R, 2V or R $\rightarrow$ A solutions (blue, plum and yellow colours, respectively). Transitions between these regimes are characterised by fold or Hopf bifurcations, indicating shifts in the stability of solution branches. By calculating the oscillation amplitude $A$ of the root-mean-square velocity of transient oscillatory solutions near selected Hopf bifurcations, we observe that these ultimately reach a stable limit cycle, $A \sim \sqrt {{Ma} - {Ma}_{{cr}}}$ . This allows us to classify these Hopf bifurcations as supercritical. Under the assumptions and flow conditions of this subsection, we have not observed any subcritical Hopf bifurcations.

It is important to note the absence of stable quasi-stationary solutions beyond the upper Hopf bifurcation curves. In this region, the flow field is dominated by nonlinear effects, which can lead to seemingly chaotic behaviour, as previously observed in both experimental and numerical studies (Machrafi et al. Reference Machrafi, Rednikov, Colinet and Dauby2010; Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ). Additionally, the upper Hopf bifurcation curves are not smooth, contrary to the other bifurcation curves. This irregularity arises from the coexistence of multiple solution branches that lose stability at different $\textit{Ma}$ values, resulting in a complex and non-smooth bifurcation structure. Interestingly, the bifurcation diagram reveals a network of bistable regions. Within these areas, the flow field can adopt solutions from different regimes based on initial conditions, demonstrating significant hysteresis.

The calculated set of solutions provides a comprehensive overview of the stability of the quasi-stationary axisymmetric solutions within binary droplets with negative $\textit{Ma}$ . However, the azimuthal stability of these solutions has not yet been addressed in this section, which is the focus of the following sections.

3.1. Regimes with vortices from rim to apex or with two vortices

We first focus on the 2V regime (figure 4 a). All of these base solutions exhibit an instability for $m = 1$ . If instabilities are triggered, the flow field will be dominated by a single vortex. The upper vortex merges with its counterpart in the rim region of the droplet, resulting in a single vortex (see bottom right of figure 4 a). Consequently, the 2V flow field is merely an artefact of the imposed axisymmetry during the computation of the base solutions.

Figure 4. Azimuthal bifurcation diagram for the stability of the quasi-stationary axisymmetric solutions in the 2V (a) and R $\rightarrow$ A (b) solution regimes. In (a), the bifurcations that limit the stability of the A $\rightarrow$ R and R $\rightarrow$ A solution regimes are shown with lower opacity. Similarly, in (b), the bifurcations for the A $\rightarrow$ R and 2V solutions are depicted with lower opacity (see all curves in figure 3). In the 2V regime, all base solutions are unstable for $m=1$ , leading to a single vortex, as depicted at the bottom right of (a). The 2V regime is therefore interpreted as an artefact of the imposed axisymmetry in § 2.3, where the upper vortex merges with its counterpart in the rim region of the droplet, as depicted at the bottom left of (a). In the R $\rightarrow$ A regime, multiple bifurcations are observed corresponding to the instability of modes $m=1$ to $m=8$ . The adjacent plots show an isometric view of an azimuthally stable solution (bottom right), and solutions subject to the linear effects of $m=2$ (top right) and $m=5$ (top left) azimuthally unstable perturbations. Exclusively in the R $\rightarrow$ A regime, the eigenvalues and eigenfunction had a non-zero imaginary part, indicating rotational motion, as depicted by the green arrows in the two upper plots.

In contrast, the R $\rightarrow$ A vortex regime (see figure 4 b) presents a variety of bifurcations, along with a stable axisymmetric regime (e.g. at $\theta =140^\circ$ and $-{Ma}=100$ , as shown in the bottom right of figure 4 b). For lower $\theta$ within this regime, an $m=1$ bifurcation at a relatively low $-{Ma}_{{cr}} \sim 100$ breaks the axisymmetry. Interestingly, the axisymmetric stability can be regained by sufficiently increasing $|{Ma}|$ . However, if $|{Ma}|$ continues to rise, the base solution becomes unstable for $m=2$ , and subsequently for $m=3$ up to $m=8$ . As $|{Ma}|$ exceeds the critical value for each bifurcation, the axisymmetry can be disrupted in a corresponding mode $m$ , as depicted in the adjacent plots of figure 4(b). These plots are derived by expanding the $\xi$ field into a sum of the base solution and a small amplitude perturbation in the direction of the corresponding eigenvector that triggers the instability of mode $m$ . While in the long-term limit the resulting $\xi$ profile will be governed by nonlinear coupling of additional excited modes, we only represent the linear effects of each mode to offer a clear visualisation of the azimuthal instability. We highlight the presence of a non-zero imaginary part in the eigenvalue and eigenfunction of the velocity field for the R $\rightarrow$ A regime, indicating rotational motion in the perturbed flow field, as depicted by the green arrows in the top plots. This result is consistent with the reports of Babor & Kuhlmann (Reference Babor and Kuhlmann2023) for thermal Marangoni flow in sessile droplets.

3.2. Regime with vortex from apex to rim

In the regime of the A $\rightarrow$ R solution, we observe multiple bifurcations for various $m$ values (see figure 5). For larger $\theta$ , greater than $45^\circ$ , the axisymmetric base solution remains stable for $|{Ma}|$ values below a critical $m=1$ bifurcation (see adjacent plot at the bottom right of figure 5), while it becomes unstable for $m=1$ above this bifurcation (see adjacent plot at the top right of figure 5). Additionally, there exists a narrow region where certain combinations of $\textit{Ma}$ and $\theta$ lead to an $m=2$ unstable solution. This region is observed for $\theta$ values slightly above $90^\circ$ and $\textit{Ma} \sim -10^4$ .

Figure 5. Azimuthal bifurcation diagram assessing the stability of the quasi-stationary axisymmetric solutions in the A $\rightarrow$ R regime. The bifurcations that limit the stability of the 2V and R $\rightarrow$ A solution regimes are shown with lower opacity (see all curves in figure 3). The bifurcation curves depict the range from $m=1$ to $m=30$ . The adjacent plots show a top view of the droplet in a stable axisymmetric regime (bottom right) or subject to the linear effects of $m=1$ (top right), 10 (bottom left) and 20 (top left) unstable perturbations on the azimuthal instability.

For lower $\theta$ , we observe a multitude of bifurcations, as illustrated in figure 5, ranging from $m=1$ to $30$ (instabilities with $m\gt 30$ are not indicated in figure 5). The bottom left and top left adjacent plots provide representations for $m=10$ and $20$ , respectively. Interestingly, the critical Marangoni number $|{Ma}_{\textit{cr}}|$ for the onset of all instabilities is particularly low for smaller $\theta$ . This observation contradicts the explanation given for Marangoni instabilities in a two-dimensional box in § 1. In a two-dimensional box, the instabilities are driven by a positive feedback loop, where a perturbation locally reduces surface tension at the interface, inducing a Marangoni flow. This, in turn, generates a return flow in the bulk, drawing more fluid with low surface tension towards the perturbed region. However, in the case of flat droplets (i.e. small $\theta$ ), the bulk domain becomes so shallow that vertical diffusion averages out any considerable compositional gradients in height direction, consistent with the assumptions of lubrication theory. Thus, the vertical return flow is absent, breaking the positive feedback loop. Contrary to expectations that Marangoni instabilities would be suppressed for low $\theta$ values, we observe the opposite behaviour.

To deeply understand this intriguing occurrence, we introduce a lubrication model in the following section. This model reduces the number of parameters to one, namely just a Marangoni number. In the limit of small $\theta$ , the contact angle can be absorbed in the non-dimensionalisation.

4.1. Lubrication model and its control parameter

In the limit of small $\theta$ , the vertical gradients are negligible compared with the radial or azimuthal ones. Rather than solving the full axisymmetric equations, we can then simplify the problem by applying lubrication theory. This allows us to simply consider a one-dimensional grid using a polar coordinate system. Alternatively, a two-dimensional circular mesh could be used, taking a Cartesian coordinate system, albeit with an increased computational cost. We choose the first option in this section.

In the lubrication approximation of an evaporating droplet (Eggers & Pismen Reference Eggers and Pismen2010; Diddens et al. Reference Diddens, Kuerten, van der Geld and Wijshoff2017a ), the governing equations for the height profile $h$ , pressure $p$ and vertically averaged mass fraction $y_A$ are given by, assuming a lateral flow profile parabolic in the height direction:

(4.1) \begin{align}&\qquad\qquad\quad\, \partial _t h = -\boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{Q} - \dfrac {j}{\rho }, \end{align}

(4.2) \begin{align}&\qquad\qquad\quad p + \boldsymbol{\nabla } \boldsymbol{\cdot } (\sigma \boldsymbol{\nabla } h) = 0, \end{align}

(4.3) \begin{align}& \partial _t (h y_A) + \boldsymbol{\nabla }\boldsymbol{\cdot } (\boldsymbol{Q} y_A) = \boldsymbol{\nabla } \boldsymbol{\cdot } (D h \boldsymbol{\nabla } y_A) - \dfrac {j_A}{\rho }, \end{align}

(4.4) \begin{align}&\qquad\qquad \boldsymbol{Q} = -\dfrac {h^3}{3 \mu } \boldsymbol{\nabla } p + \dfrac {h^2}{2 \mu } \boldsymbol{\nabla } \sigma , \end{align}

where $j$ is the evaporation rate, $\rho$ the mass density, $D$ the diffusion coeficient and $\mu$ dynamic viscosity. Note that the transient effects are more pronounced for flat droplets as compared with droplets with higher contact angle. However, we calculate the quasi-stationary solutions of the lubrication model for simplicity. We assume that the evaporation of one component is balanced by the condensation of the other in the quasi-stationary limit, such that the total volume remains constant. Thereby, we consider the evaporation rate $j_A$ to be given by the theoretically calculated limit for a flat pure droplet (Popov Reference Popov2005):

(4.5) \begin{equation} j_A = \dfrac {2 D_A^{\textit{vap}} (c_A^{{eq}} - c_A^\infty )}{{\unicode{x03C0}} R \sqrt {1 -\tilde{r}^2}} \, , \end{equation}

where $\tilde{r}$ is the dimensionless radial coordinate. As in § 2.1, we assume that the droplet maintains its initial shape, i.e. we consider the limit of zero capillary number ${\textit{Ca}}=0$ and $\partial _t h = 0$ . In the limit of the lubrication theory, the equilibrium shape corresponds to a parabolic height profile. We emphasise that these approximations may limit the applicability of the results in real systems, but they capture the essential physics that explain the onset of the Marangoni instability. By reformulating $\xi = y \theta / {Ev}$ and applying the mentioned assumptions, while using the non-dimensional scales introduced in § 2.1 – but now modifying the length scale to the base radius $R$ instead of $V^{1/3}$ – we can reduce the system of (4.1)–(4.4) to

(4.6) \begin{align}&\qquad\qquad\quad \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } \tilde {\boldsymbol{Q}} = 0, \end{align}

(4.7) \begin{align}& {\tilde {h} \partial _{\tilde {t}} \xi + } \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } (\tilde {\boldsymbol{Q}} \xi ) = \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } (\tilde {h} \tilde {\boldsymbol{\nabla }} \xi ) - \tilde {j}, \end{align}

(4.8) \begin{align}&\quad \tilde {\boldsymbol{Q}} = - \dfrac {\tilde {h}^3}{3} \tilde {\boldsymbol{\nabla }} \tilde {p} + {Ma}_R \dfrac {\tilde {h}^2}{2} \tilde {\boldsymbol{\nabla }} \xi . \end{align}

Here, $\tilde {h} = 0.5 (1-\tilde{r}^2)$ is the dimensionless height profile, $\tilde {j} = 2/({\unicode{x03C0}} \sqrt {1-\tilde{r}^2})$ is the evaporation rate in the flat droplet limit and $\textit{Ma}_R$ is the modified Marangoni number, defined as

(4.9) \begin{align} {Ma}_R = \dfrac {R \partial _{y_A} \sigma }{D_0 \mu _0} {Ev}. \end{align}

We impose a homogeneous bulk correction for $\xi$ , i.e. $\langle \xi \rangle = 0$ , to allow for quasi-stationary solutions with a fixed compositional average. While the system of (4.6)–(4.8) can be solved analytically (see Appendix A), introducing a small amplitude $\epsilon$ of an azimuthal perturbation to the base solution makes the augmented system too complex for analytical solutions. Therefore, we employ numerical methods to determine the azimuthal stability of the quasi-stationary axisymmetric solutions and compare these results with those obtained using the full model described in § 2.1. As $\theta \rightarrow 0$ , the results from the full model converge to the lubrication model, validating the latter (see figure 6 b). Remarkably, the number of unstable modes $m$ is found to be consistent in both models. Interestingly, this number is larger for low $\theta$ values than for large $\theta$ values (see figure 6 a), i.e. a cascade of azimuthal instabilities sets in for decreasing $\theta$ at unexpectedly small negative Marangoni numbers. By expanding the $\xi$ field into a sum of the base solution and a small amplitude perturbation in the direction of the corresponding eigenvector that triggers the instability of mode $m$ , we can visualise the linear effects of the azimuthal instability for the lubrication model.

Figure 6. (a) Onset of azimuthal instabilities for different $m$ as a function of $\textit{Ma}_R$ in the limit of small $\theta$ . The adjacent plots show a top view of the droplet subject to the linear effects $m=5$ (top left), $m=20$ (bottom left) and $m=60$ (bottom right) unstable perturbations on the azimuthal instability. (b) Comparison of the azimuthal stability of the quasi-stationary axisymmetric solutions between the full model as $\theta \rightarrow 0$ (solid) and the lubrication model (dashed).

4.2. Driving mechanism of the azimuthal instabilities

In the absence of bulk gradients, the driving mechanism of the Marangoni instabilities for low $\theta$ values is different from that described in the two-dimensional box of § 1, yet it bears resemblance.

Figure 7. Side (a) and top (b) views of an evaporating droplet prior to the onset of the Marangoni instability. In this figure, all variables and operators are dimensionless, but tildes were dropped for brevity. The sign of the pressure gradients indicate whether the flow is dominated by Marangoni flow (positive pressure gradient) or pressure-driven flow (negative pressure gradient). Should a disturbance locally reduce the surface tension at the contact line, it triggers an azimuthal Marangoni flow (yellow arrows) that transports fluid away from the perturbed region. This initiates pressure-driven return flow (green arrow) that attracts more fluid with lower surface tension from the centre of the droplet (greater-height region) towards the contact line (lower-height region). Ultimately, an azimuthal Marangoni instability is triggered (c).

For droplets with small $\theta$ and a negative Marangoni number ( $\textit{Ma}_R$ ), the larger evaporation rate at the contact line leads to an enhanced concentration of fluid with greater surface tension. Should a disturbance locally reduce the surface tension at contact line, it triggers an azimuthal Marangoni flow that transports fluid away from the perturbed region. Continuity of the fluid initiates a pressure-driven return flow, which attracts more fluid with lower surface tension towards the contact line (see figure 7). This creates a positive feedback loop, which is the driving mechanism of the azimuthal instabilities observed in droplets with low contact angle $\theta$ .

Naturally, the same reasoning still holds for droplets with larger contact angle $\theta$ . However, in this latter case, the presence of bulk gradients becomes more pronounced, primarily driving the Marangoni instabilities. A simple observation of the scaling of the different terms in (4.8) reveals that the flow induced by $\textit{Ma}_R$ scales with $\tilde{h}^2$ , while the flow driven by pressure scales with $\tilde{h}^3$ . As a consequence, in areas where the height profile is minimal, the Marangoni flow prevails, but in areas with a significant height profile, pressure-driven flow takes precedence. Thus, in droplets with a low $\theta$ , the flow near the contact line is driven predominantly by Marangoni flow, whereas pressure-driven circulation dominates the centre (as illustrated by the arrows in figure 7).

Therefore, the geometry plays a fundamental role in the driving mechanism of the onset of azimuthal instabilities. In the next subsection, we propose a geometry where such instabilities can be obtained for positive- $\textit{Ma}_R$ droplets (e.g. glycerol–water mixtures).

4.3. Instabilities in positive Marangoni mixtures

The insights gained from prior analyses elucidate that, in droplets with small $\theta$ , if the fluid with larger surface tension has the lowest height profile, azimuthal instabilities can be triggered. For positive- $\textit{Ma}_R$ mixtures, the enhanced evaporation at the contact line leads to higher surface tension in the centre, which is the zone of larger height. Thereby, positive- $\textit{Ma}_R$ mixtures consistently exhibit stable axisymmetric solutions. Any disturbance reducing the surface tension at the centre is promptly neutralised by a Marangoni flow in the opposite direction.

Figure 8. Side (a) and top (b) views of a droplet placed on a shallow pit. Due to the inverted height profile, $\textit{Ma}_R$ (yellow arrows) is predominant at the centre over the pressure-driven flow (green arrows), resulting in the onset of azimuthal Marangoni instabilities, depicted in (b). After enough evaporation time (6 s), blobs emerge close to the contact line in the realistic case of an $R=$ 0.5 mm, $\theta =10^\circ$ , positive- $\textit{Ma}_R$ droplet in an $h_R= 0.5\,\unicode{x03BC} \textrm {m}$ pit, with initial $y_{\textit{glycerol}}= 25\,\%$ (c).

Altering the scenario by placing the droplet within a shallow pit with a concave liquid–gas interface reverses the height profile, making the centre the lowest point (see figure 8 a). This set-up assumes that the depression is shallow enough for lubrication theory to hold and that the fluid remains pinned to the pit’s sidewalls. Under these conditions, a perturbation reducing surface tension at the centre instigates a Marangoni flow that shifts fluid azimuthally away from the disturbed area, thereby initiating an azimuthal Marangoni instability akin to the mechanism described in § 4.2, albeit in a reversed direction (see figure 8 b). This behaviour is unexpected when following the explanation provided in § 1 for Marangoni instabilities in a two-dimensional box. In that scenario, positive- $\textit{Ma}$ mixtures cannot initiate Marangoni instabilities.

To simulate the evaporation dynamics of a realistic glycerol–water droplet with a radius of $R= 0.5\,\textrm {mm}$ , $\theta = 10^{\circ }$ , in a pit of depth $h_R= 0.5 \, \unicode{x03BC} \textrm {m}$ , starting with a $y_{\textit{glycerol}}= 25\,\%$ , we apply (4.1)–(4.4) in a two-dimensional circular domain. As evaporation progresses, blobs begin to form near the rim of the well due to the azimuthal Marangoni instability (see figure 8 c). We emphasise that the pit must be shallow enough for the onset of instabilities, otherwise bulk gradients would stabilise the flow. While this finding is remarkable, a detailed investigation of the azimuthal Marangoni dynamics of droplets in such a well is beyond the scope of the work.

4.4. Capillary effects in the Marangoni instability

Until now, we assumed that the droplet retains a perfectly spherical-cap shape (parabolic shape in the lubrication limit) during evaporation, neglecting capillary effects on the onset of instabilities. In this subsection, we explore non-zero capillary numbers to examine the influence of capillary effects, in particular shape deformations induced by surface-tension gradients and Marangoni flow, on the onset of azimuthal instabilities. Utilising the scaling introduced in § 4.1 and defining the capillary number as

(4.10) \begin{equation} \textit{Ca} = \dfrac {D_0 \mu _0}{\theta R \sigma _0}, \end{equation}

we can obtain the non-dimensionalised form of (4.1)–(4.4), namely,

(4.11) \begin{align}&\qquad\qquad\quad \partial _{\tilde {t}} \tilde {h} + \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } \tilde {\boldsymbol{Q}} = 0, \end{align}

(4.12) \begin{align}& \tilde {p} + \boldsymbol{\nabla } \boldsymbol{\cdot } \big (\theta ^2 (\textit{Ca}^{-1} + {Ma}_R \xi ) \boldsymbol{\nabla } \tilde {h} \big ) = 0, \end{align}

(4.13) \begin{align}&\,\, \partial _{\tilde {t}}(\xi \tilde {h}) + \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } (\tilde {\boldsymbol{Q}} \xi ) = \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } (\tilde {h} \tilde {\boldsymbol{\nabla }} \xi ) - \tilde {j}. \end{align}

The definition of $\tilde {\boldsymbol{Q}}$ remains the same as in (4.8). Again, we neglect the ‘coffee-stain effect’ in equation (4.11), in the sense that we assume that the volume of the droplet does not change. By solving the system of equations (4.11)– (4.13) for various capillary numbers and contact angles of $10^\circ$ (figure 9 a) and $20^\circ$ (figure 9 b), we can identify the critical $\textit{Ma}_R$ for the onset of azimuthal instabilities across different wavenumbers $m$ (see figure 9) and understand how the results of figure 6 are affected by shape deformations. We ensure that the capillary number $\textit{Ca}$ is set such that the effective dimensionless surface tension, given in (4.12) by $\sigma = \theta ^2 (\textit{Ca}^{-1} + {Ma}_R \xi$ ), remains positive. Beyond this limit, we discard the solutions as they are physically unrealistic, as depicted by the grey areas of figures 9(a) and 9(b), limited by the black dashed lines.

Figure 9. Critical Marangoni number ( ${\textit{Ma}}_R$ ) for the onset of azimuthal instabilities as a function of the capillary number $\textit{Ca}$ for wavenumbers $m=$ 1, 2, 3 and 4 and for contact angles of $10^\circ$ (a) and $20^\circ$ (b). The grey area, limited by the black dashed line, indicates where the surface tension is negative, rendering physically unrealistic solutions, therefore discarded from the analysis. The pink and yellow regions highlight the areas where the droplet maintains a nearly spherical-cap shape and where it transitions to a pancake shape, respectively.

The results indicate that, within the physical constraints of the capillary number ( $\textit{Ca}$ ), the critical Marangoni number ( $\textit{Ma}_R$ ) for the onset of azimuthal instabilities is significantly influenced by shape deformations. For sufficiently large Marangoni numbers, the droplet adopts a pancake shape, as observed in both experimental and numerical studies (Diddens et al. Reference Diddens, Kuerten, van der Geld and Wijshoff2017a ; Pahlavan et al. Reference Pahlavan, Yang, Bain and Stone2021; Yang et al. Reference Yang, Pahlavan, Stone and Bain2023). In figure 9(a,b), we highlight the regions in the phase diagram where the droplet maintains a nearly spherical-cap shape (pink) and where it transitions to a pancake shape (yellow). The critical transition between these shapes was determined by identifying the set of parameters ( $\textit{Ma}_R$ , $\textit{Ca}$ ) at which the vertically averaged pressure at the highest point of the droplet along the axis of symmetry becomes zero.

An increase in the modulus of the critical Marangoni number is observed across all azimuthal wavenumbers $m$ , with a more pronounced effect for $m=1$ and for $\theta =10^\circ$ , where $|{Ma}_R|$ increases from approximately $200$ to $10^5$ as $\textit{Ca}$ rises from 0 to approximately $10^{-3}$ . Interestingly, the effects of shape deformations are less pronounced for droplets with larger contact angles, see figure 9(b). In evaporating droplets with positive Marangoni number and with unpinned contact line, the Marangoni flow can drive fluid towards the apex strongly enough to contract the droplet (‘Marangoni contraction’, Karpitschka, Liebig & Riegler Reference Karpitschka, Liebig and Riegler2017). Similarly, a pinned droplet with negative Marangoni number can form a pancake shape (Pahlavan et al. Reference Pahlavan, Yang, Bain and Stone2021; Yang et al. Reference Yang, Pahlavan, Stone and Bain2023), where strong Marangoni flow drives fluid from the apex of the droplet towards the contact line, causing deformation and an increase in the height profile at the contact line. The microscopic contact angle also becomes higher for pancake droplets. This agrees with the observations of our model: instabilities for flat droplets start at the rim and the higher the contact angle, the higher $|{Ma}_R|$ is required for instabilities.

6.1. Slowly evaporating axisymmetric water–ethanol droplet

We consider a slowly evaporating axisymmetric water–ethanol droplet. The mass fractions of water, ethanol and air in the farfield of the gas phase are adjusted at each step to maintain constant ${Ev}=0.001$ and ${Ev}_{{tot}}=0.005$ . These values are chosen to ensure that $\textit{Ma}$ remains within the limit in which stable axisymmetric quasi-stationary solutions exist and that the thermal Marangoni number $\textit{Ma}_T$ remains at least one order of magnitude lower than the solutal Marangoni number $\textit{Ma}$ , while still allowing for the droplet to evaporate.

The droplet has an initial contact angle of $\theta = 120^\circ$ and an initial volume of $V_0 = 0.1\,\unicode{x03BC} \textrm {l}$ . The contact line is pinned. The substrate has a thickness of 0.1 mm, it is made of borosilicate, and it has air beneath it. Thermal effects are incorporated into the transient model, as well as the effects of local changes in liquid composition on the evaporation rate, Stefan’s flow and natural convection in the gas phase. At each time step, the solutal Marangoni number $\textit{Ma}$ and the contact angle are calculated. These are then used as input parameters for the quasi-stationary model of § 2.1. As figure 3 shows, the quasi-stationary solutions exhibit hysteresis. To directly compare the results of the full transient simulations with the quasi-stationary solutions, we use the results of the transient model as an initial condition for the quasi-stationary calculations. The results are shown in figure 11. Supplementary movies 1 and 2 (available at https://doi.org/10.1017/jfm.2025.10519) show the full transient evolution of the droplet and its comparison with the quasi-stationary model (movie 1 for composition, movie 2 for velocity magnitude the temporal evolution of the temperature field is shown in supplementary movie 3).

Despite the absence of transient effects, thermal effects, the effects of local changes in liquid composition on the evaporation rate, Stefan’s flow and natural convection in the gas phase, the quasi-stationary solutions are in good agreement with the full transient simulations. In the transient model, the lower temperature at the apex of the droplet induces a thermal Marangoni flow that drives fluid from the apex to the rim. It also decreases the vapour pressure of ethanol at the apex, leading to lower evaporation at the apex – opposite to what is seen in the quasi-stationary model. Nevertheless, the flow is driven predominantly by solutal Marangoni effects, as demonstrated by the strong agreement of flow direction and compositional profile between the two models. In figure 11(a) and 11(b), the flow is in the R $\rightarrow$ A direction, in figure 11(c), it is in the 2V direction, and in figure 11(d), it is in the A $\rightarrow$ R direction, showing the presence of all three quasi-stationary solution regimes.

Figure 11. Slowly evaporating water–ethanol droplet with an initial contact angle of $\theta = 120^\circ$ and an initial volume of $V_0 = 0.1\,\unicode{x03BC} \textrm {l}$ . The contact line is pinned to the substrate. The far-field mass fractions of ethanol and water vapours are adjusted such that the evaporation number ${Ev} = 0.001$ and the total evaporation number ${Ev}_{{tot}}=0.005$ . Different stages of the droplet’s evolution are shown at $t = 100\,\textrm {s}$ (a), $t = 7000\,\textrm {s}$ (b), $t = 11000\,\textrm {s}$ (c) and $t = 15000\,\textrm {s}$ (d). Each panel shows the mass fraction of ethanol within the droplet and the mass fraction of ethanol in the vapour phase. Results are presented for both the full transient simulations on the left of each panel (i) and the corresponding quasi-stationary solutions on the right (ii).