3.1. Unveiling instabilities in diffuse interfaces

We identify an instability that arises solely in diffuse interfaces and is absent in the sharp-interface limit, driven by Marangoni stresses. We term this the diffuse-interface Marangoni instability.

As shown in figure 2(a), short-wavelength instabilities emerge when the interface has finite thickness, even under conditions where sharp-interface theory predicts stability (equal viscosity and diffusivity in both phases) (Sternling & Scriven Reference Sternling and Scriven1959; Schwarzenberger et al. Reference Schwarzenberger, Köllner, Linde, Boeck, Odenbach and Eckert2014). For $ \textit{Ma} = 10^4$ , $\textit{Ca} = 0.01$ , $\alpha = 1.7$ and $Sc=1000$ , an unstable band appears at $\delta =0.2$ and narrows as $\delta$ decreases, disappearing entirely as $\delta$ approaches $0.02$ – an estimated sharp-interface limit. This limit is calculated by using ‘pyoomph’ (Diddens & Rocha Reference Diddens and Rocha2024), following the methodology described in Bauer, Frink & Kreckel (Reference Bauer, Frink and Kreckel2002), Heil & Hazel (Reference Heil and Hazel2006). Disabling Marangoni stresses by omitting the $-\text{i}k\tilde {c} W^b$ term in (2.3c ) suppresses the instability, yielding a single $\delta$ -independent growth curve (blue triangles in figure 2 a), which agrees well with the sharp-interface limit. This confirms that the observed instability stems from solutal Marangoni effects rather than numerical artefacts.

Figure 2.

(a) Growth rate $\omega _i$ versus wavenumber $k$ for different $\delta$ . The grey-shaded area indicates $\omega _i \gt 0$ , i.e. instability. Stability is restored when Marangoni force is off. Parameters are $ \textit{Ma}=10^4$ , $\textit{Ca}=0.01$ , $\alpha =1.7$ . The black solid line with square markers shows the corresponding sharp-interface result. (b–e) Instability intervals (enclosed regions) under varying $\textit{Sc}$ , $\textit{Ca}$ , $\textit{Pe}_\phi$ , and $ \textit{Ma}$ , respectively.

To probe the instability mechanism, we vary key parameters individually: Marangoni number $ \textit{Ma} \in [10^3, 10^5]$ , capillary number $\textit{Ca} \in [10^{-3}, 10^{-1}]$ , Schmidt number $Sc \in [10^2, 10^6]$ and mobility index $\alpha \in [1.3, 1.7]$ . These span $Re \in [10^{-2}, 10^2]$ and $\textit{Pe}_\phi \in [0.3, 1.6\times 10^5]$ , with $\textit{Re} = \textit{Ma}/\textit{Sc}$ and $\textit{Pe}_\phi = \textit{UL}/\delta ^{\alpha }$ . Default parameters are $ \textit{Ma} = 10^4$ , $\textit{Ca} = 0.01$ , $\alpha = 1.5$ and $Sc = 1000$ , chosen based on typical physical values to ensure realistic, weakly deformed interfaces, with $\alpha$ set according to established recommendations (Jacqmin Reference Jacqmin1999; Magaletti et al. Reference Magaletti, Picano, Chinappi, Marino and Casciola2013).

Figure 2(b–e) shows the neutral stability curves $k_{\textit{NS}}$ versus interfacial thickness $\delta$ . Here, $\textit{Sc}$ has little influence, consistent with prior studies (Reichenbach & Linde Reference Reichenbach and Linde1981; Li et al. Reference Li, Wan, Zhang and Tan2025). Accordingly, inertia is negligible for the instability threshold within the $Re$ studied range ( $Re = 0.01$ –100). In contrast, $ \textit{Ma}$ , $\textit{Ca}$ and $\textit{Pe}_\phi \sim \delta ^{-\alpha }$ significantly alter the instability. Lower $\textit{Ca}$ or higher $\alpha$ suppresses short-wavelength modes but leaves the critical thickness $\delta _{\textit{cr}}$ unchanged. In contrast, reducing $ \textit{Ma}$ narrows the unstable wavenumber band and increases $\delta _{\textit{cr}}$ , shrinking the unstable regime. These results indicate that while $\delta$ , $ \textit{Ma}$ , $\textit{Ca}$ and $\textit{Pe}_\phi$ govern instability onset, $\delta _{\textit{cr}}$ is uniquely determined by $ \textit{Ma}$ .

3.2. Scaling laws for critical diffuse-interface thickness

In the laminar regime, the critical thickness $\delta _{\textit{cr}}$ depends uniquely on $ \textit{Ma}$ , not $\textit{Ca}$ or $\textit{Pe}_\phi$ , indicating that $\delta _{\textit{cr}}$ is governed primarily by solutal Marangoni advection–diffusion dynamics near the interface ((2.2b ) and (2.3b )). The circulation (figure 1 b) transports solute toward the interface, steepening concentration gradients, while diffusion tends to smooth them. This competition features the Marangoni interfacial flows (Manikantan & Squires Reference Manikantan and Squires2020). Since advection in the normal direction transports fluids of different concentrations across the layers, $\delta L$ is taken as the characteristic length scale. Balancing diffusion (time scale $t_D\sim (\delta L)^2/D$ ) and advection ( $t_A\sim \delta L /U$ ) yields the scaling

(3.1) \begin{equation} \delta _{\textit{cr}}\sim Ma^{-1}. \end{equation}

As $ \textit{Ma}$ increases, stronger advection requires a thinner interface to maintain stability. This $ \textit{Ma}$ – $\delta _{\textit{cr}}$ scaling is confirmed as $\delta _{\textit{cr}} \rightarrow 0$ in figure 3(a), and remains robust under changes in $\textit{Ca}$ or $\textit{Pe}_\phi$ , as shown in figures 2(c) and 2(d). This relation provides a physical criterion for selecting $\delta$ in diffuse-interface simulations with Marangoni effects, ensuring convergence to the sharp-interface limit (Jacqmin Reference Jacqmin1999). While thin-interface limit criteria exist for passive two-phase flows (Magaletti et al. Reference Magaletti, Picano, Chinappi, Marino and Casciola2013; Demont et al. Reference Demont, Stoter and Van Brummelen2023), the present finding extends them to Marangoni-driven instabilities, consistent with thermocapillary flow observations by Verschueren et al. (Reference Verschueren, Van De Vosse and Meijer2001).

Figure 3.

(a) Scaling of the critical interface thickness $\delta _{\textit{cr}}$ with $ \textit{Ma}$ , showing $\delta _{\textit{cr}} = 590/\textit{Ma}$ for $ \textit{Ma} \gtrsim 4000$ , independent of $\textit{Ca}$ and $\alpha$ . Deviations from the scaling emerge with increasing $\alpha$ or decreasing $\textit{Ca}$ . (b) Energy budget showing the relative contributions of $\langle\mathcal{E}_i\rangle = \langle \mathcal{E}_\eta \rangle$ (viscous dissipation), $\langle\mathcal{E}_\gamma^{cap}\rangle$ (capillary forces), $\langle\mathcal{E}_\gamma^{tilt}\rangle$ (interfacial tilt), compared to the solutal Marangoni effect $\langle\mathcal{E}_\gamma^{mag}\rangle$ , as $\delta _{\textit{cr}}$ increases. (c,d) Base-state phase field $\phi ^b_{\textit{s}}$ and solute perturbation $c^\prime$ for the stable case ( $\delta = 0.05$ ) and unstable case ( $\delta = 0.8$ ). Velocity vectors are overlaid on $\phi ^b_{\textit{s}}$ ; streamlines are shown in $c^\prime$ . The temporal evolutions from DNS are given in supplementary movies S1 and S2, respectively, available at https://doi.org/10.1017/jfm.2025.10895. (e) Cross-sectional profiles of $u$ (at $x = \pi$ ) and $c^\prime$ (at $x = \pi /2$ ) along the $y$ -direction for the stable (solid lines) and unstable (dashed lines) cases.

To explain the deviation from the scaling law (3.1) at larger interfacial thickness (i.e. for $ \textit{Ma} \lesssim 4000$ ), we perform an energy-based analysis to quantify the competing mechanisms. Multiplying the linearised perturbation equations ((2.3c ) and (2.3d )) by the velocity perturbations yields an energy budget:

(3.2) \begin{equation} \begin{aligned} \frac {\partial E_{\textit{kin}}}{\partial t} = &\underbrace {\frac {{1}}{Re} \left [-\left (\boldsymbol{\nabla }u^{\prime }\right )^2-\left (\boldsymbol{\nabla }v^{\prime }\right )^2\right ] }_{\mathcal{E}_\eta :\, \mathrm{Viscous}\,\mathrm{dissipation}} + \underbrace {\frac {{1}}{Re}\frac {1}{\textit{Ca}} \left [ \frac {\varGamma ^b}{\left |\boldsymbol{\nabla }\phi _{\textit{s}}^b\right |} \frac {\partial ^2 \phi _{\textit{s}}^{\prime }}{\partial x^2} v^{\prime } \right ] W^b}_{\mathcal{E}^{\textit{cap}}_\gamma :\, \mathrm{Capillary}\,\mathrm{effect}} \\ &+ \underbrace {\frac {{1}}{Re}\frac {1}{\textit{Ca}}\left [\frac {1}{\left |\boldsymbol{\nabla }\phi _{\textit{s}}^b\right |} \frac {\partial \phi _{\textit{s}}^{\prime }}{\partial x} \frac {\mathrm{d} \varGamma ^b}{\mathrm{d} y} u^{\prime }\right ] W^b }_{\mathcal{E}^{\textit{tilt}}_\gamma :\, \mathrm{Tilt}-\mathrm{induced}\,\mathrm{Marangoni}}+ \underbrace { \frac {{1}}{Re}\left [-\frac {\partial c^{\prime }}{\partial x} u^{\prime }\right ] W^b }_{\mathcal{E}^{\textit{mag}}_\gamma :\, \mathrm{Solutal}\,\mathrm{Marangoni}\,\mathrm{effect}}. \end{aligned} \end{equation}

Here, the kinetic energy density is defined as $E_{kin} = ( u^{\prime^2} + v^{\prime^2} )/2$ , since the base flow $\boldsymbol{U}^b$ is zero and total velocity reduces to the perturbation field ( $\boldsymbol{u}=\boldsymbol{u}^\prime$ ). The prefactor $ \textit{Ma}/Sc$ is replaced by $Re\lesssim 4$ .

The viscous dissipation $\mathcal{E}_\eta$ reflects energy loss due to shear. The solutal Marangoni term $\mathcal{E}_\gamma ^{\textit{mag}}$ originates from tangential solute gradients $-\partial c^\prime /\partial x$ . The capillary term $\mathcal{E}_\gamma ^{\textit{cap}}$ arises from interface curvature and base-state tension $\varGamma ^b$ . The tilt-induced Marangoni term $\mathcal{E}_\gamma ^{\textit{tilt}}$ stems from interfacial tilt $\partial \phi^\prime_s / \partial x$ , which projects normal interfacial tension gradients $\mathrm{d} \varGamma ^b/\mathrm{d} y$ across the diffuse layer (Berry Reference Berry1971; Marchand et al. Reference Marchand, Weijs, Snoeijer and Andreotti2011) onto the tangential direction – an effect inherently tied to the finite thickness of diffuse interfaces and absent in sharp-interface models.

Along the $ \textit{Ma}$ – $\delta _{\textit{cr}}$ paths (figure 3 a), integrating over the domain $[2\pi /k \times 4L]$ gives the power balance $\langle \mathcal{E}_\eta \rangle +\langle \mathcal{E}_\gamma ^{\textit{cap}}\rangle +\langle \mathcal{E}_\gamma ^{\textit{tilt}}\rangle +\langle \mathcal{E}_\gamma ^{\textit{mag}}\rangle =0$ ( $\omega _i=0$ ), where the angled brackets $\langle \boldsymbol{\cdot }\rangle$ denote a spatial average over the entire computational domain. Divergence terms vanish under periodic, homogeneous boundary conditions (Zhang et al. Reference Zhang, Lashgari, Zaki and Brandt2013). Figure 3(b) shows the normalised contributions of each term by $\langle \mathcal{E}_\gamma ^{\textit{mag}}\rangle$ as $\delta _{\textit{cr}}$ increases. For small $\delta _{\textit{cr}}$ ( $\lesssim 0.1$ ), the balance is primarily between viscous dissipation $\mathcal{E}_\eta$ and solutal Marangoni forcing $\mathcal{E}_\gamma ^{\textit{mag}}$ , consistent with the scaling law. As $\delta _{\textit{cr}}$ increases, the influence of capillary $\mathcal{E}_\gamma ^{\textit{cap}}$ starts at a certain $\delta _{\textit{cr}}$ , beyond which deviations from the ideal $\delta _{\textit{cr}} \sim \textit{Ma}^{-1}$ scaling become significant. Thus, interface deformation affects the limits of the scaling law’s validity.

Figures 3(c) and 3(d) show the base-state phase field $\phi ^b_{\textit{s}}$ (blue), velocity field $\boldsymbol{u}=\boldsymbol{u}^\prime$ (arrows) normalised by its global maximum absolute value $|u|_{max}$ , and perturbation concentration field $c^\prime$ (red) normalised by its global maximum absolute value $|c^\prime|_{max}$ for two representative cases: a stable system ( $\delta = 0.05$ ) and an unstable one ( $\delta = 0.8$ ), corresponding to the solid and open star symbols in figure 3(a), respectively. In the stable case (figure 3 c), vortices are symmetrically distributed about the interface, as evidenced by the velocity component $u$ profiles at $x=\pi$ (figure 3 e, black solid line). This balanced Marangoni recirculation gives rise to an anti-symmetric $c^\prime$ field, confirmed by the profile at $x=\pi /2$ (figure 3 e, red solid line). In contrast, the unstable case (figure 3 d) exhibits anti-symmetry breaking, with flow field shifting downward into source phase (figure 3 e, black dashed line). The zero-crossing point of $c^\prime$ shifts upward from the interface centreline ( $y=2$ ), indicating solute enrichment in the receiver phase and depletion in the source phase (figure 3 e, red dashed line). We further performed DNS (Li et al. Reference Li, Wan, Zhang and Tan2025) for both cases to reveal the evolution of the different fields, with the results provided in the supplementary movies. These symmetry-breaking shifts in both velocity and concentration profiles motivate the next section, in which we investigate their physical origin and impact on the instability.

3.3. Mechanism of instability in diffuse interfaces: anti-flux versus pro-flux shifts

In classical sharp-interface systems with equal viscosity and diffusivity between two fluids, interfacial Marangoni instability is absent (Sternling & Scriven Reference Sternling and Scriven1959). As illustrated in figure 4(a), symmetric convection brings high ( $C^H$ ) and low ( $C^L$ ) concentrations toward the interfacial stagnation point. Equal diffusivity and viscosity ensure that net interfacial concentration remains unchanged, yielding no amplification of the Marangoni flow $u_m$ . In contrast, a diffuse interface introduces local asymmetry. Figure 4(b) shows that inflow of low-concentration fluid reduces interfacial concentration $C^-$ on the receiver phase $\phi _{\textit{r}}$ side, while inflow of high-concentration fluid increases $C^+$ on the source $\phi _{\textit{s}}$ side. Within the interfacial layer of thickness $\delta$ , the resulting imbalance cannot be instantly homogenised, leading to a tangential concentration gradient $\boldsymbol{\nabla}_{\kern-2pt s} C (y)$ and an associated velocity profile $u(y)$ .

Figure 4.

(a) Schematic of interfacial solute and velocity fields in the sharp-interface limit, where symmetric advection yields no net concentration change. (b) Profiles of solute concentration $C(y)$ and tangential velocity $u(y)$ across the diffuse interface. (c) Mechanism of velocity offset $\ell _\delta ^u$ : asymmetric solute advection enhances Marangoni flow on the source side, displacing the velocity peak away from the centreline. (d) Mechanism of solute offset $\ell _\delta ^{c^\prime }$ : coupled vertical flow and interfacial diffusion shift the position of zero net solute change. (e,f) Cross-sectional profiles of $u$ and $c^\prime$ for the unstable case ( $\delta = 0.8$ ) at three representative times ( $t_1=10^4$ , $t_2=5 \times 10^4$ and $t_3=16 \times 10^4$ ) from DNS, showing the offsets $\ell _\delta ^u$ and $\ell _\delta ^{c^\prime }$ . Here, $u$ and $c^\prime$ are normalised by their global maximum absolute values $|u|_{max}$ and $|c^\prime|_{max}$ , respectively. (g,h) Normalised offsets $\ell _\delta ^u/\delta _{cr}$ and $\ell _\delta ^{c^\prime }/\delta _{cr}$ versus the modified Marangoni number $ \textit{Ma}_\delta = Ma \, \delta$ , showing collapse across varying $\textit{Ca}$ and $\alpha$ . Instability sets in at $ \textit{Ma}_{\textit{cr}}^{\delta }\approx 590$ . (i) Phase diagram in the $\ell _\delta ^u$ – $\ell _\delta ^{c^\prime }$ plane. Stable (solid) and unstable (open) regimes are separated by the scaling $\ell _\delta ^u \sim (\ell _\delta ^{c^\prime })$ .

As shown in figure 4(c), the local concentration decrease on the receiver side ( $\phi _{\textit{r}}$ ) reduces the interfacial tension gradient $(\boldsymbol{\nabla}_{\kern-2pt s} \gamma )^-$ , thereby weakening the Marangoni flow (Anderson et al. Reference Anderson, McFadden and Wheeler1998). Conversely, the increased concentration on the source side ( $\phi _{\textit{s}}$ ) enhances $(\boldsymbol{\nabla}_{\kern-2pt s} \gamma )^+$ , strengthening the flow. This asymmetry in surface tension gradients induces a shear layer within the interface, with suppressed velocity near $\phi _{\textit{r}}$ and enhanced velocity near $\phi _{\textit{s}}$ . The resulting shear promotes an interfacial instability reminiscent of, yet physically distinct from, the classical Kelvin–Helmholtz instability (Chandrasekhar Reference Chandrasekhar2013), since the present mechanism arises from Marangoni-driven shear within a diffuse interface rather than inertia-dominated interfacial waves. This asymmetry in surface tension gradients induces a shear layer within the interface, with suppressed velocity near $\phi _{\textit{r}}$ and enhanced velocity near $\phi _{\textit{s}}$ . We quantify the shift of the peak velocity ( $u_m$ )  position $y_{u_m}$ away from the interface centreline ( $y=2L$ ) by the anti-flux offset, $\ell _{\delta }^{u}=2L-y_{u_m}$ , which characterises the destabilising effect of advection-driven shear. Diffusion counters this destabilisation by smoothing concentration gradients ( $C^-\uparrow$ and $C^+\downarrow$ ), as depicted in figure 4(d). The position $y_{C^0}$ where the net concentration change vanishes ( $C^0$ ) shifts along the solute flux into the receiver side. This stabilising effect is captured by the pro-flux offset $\ell _{\delta }^{c^\prime }=y_{C^0}-2L$ . Figures 4(e) and 4(f) illustrate the evolution of the velocity and concentration profiles within the diffuse interface at three representative instants, as obtained from DNS of the unstable case ( $\delta =0.8$ ), thereby confirming the mechanisms discussed.

To quantify the competition between advection and diffusion across the interface, we introduce a modified Marangoni number $ \textit{Ma}_\delta = \textit{Ma}\, \delta$ , using the interfacial thickness $\hat {\delta }$ as the characteristic length scale. As shown in figures 4(g) and 4(h), both the velocity offsets $\ell _{\delta }^{u}$ and the concentration offset $\ell _{\delta }^{c^\prime }$ increase monotonically with $ \textit{Ma}_\delta$ , following distinct power-law trends that collapse across different $\textit{Ca}\in [10^{-3}, 10^{-1}]$ and $\alpha \in [1.3, 1.7]$ , confirming their universality. A critical transition emerges at $ \textit{Ma}_\delta \approx 590$ : below this threshold, the system remain stable, with $\ell _{\delta }^{u}$ growing sublinearly and $\ell _{\delta }^{c^\prime }$ superlinearly; beyond it, the system becomes unstable, and both offsets grow with fixed exponents, approximately 1.13 for $\ell _{\delta }^{u}$ and 0.75 for $\ell _{\delta }^{c^\prime }$ .

To further capture this interplay, figure 4(i) plots $\ell _{\delta }^{u}$ against $\ell _{\delta }^{c^\prime }$ , revealing a clear separation between stable (bluish solid) and unstable (reddish open) regimes by the power-law relation $\ell _{\delta }^{u} \sim (\ell _{\delta }^{c^\prime })^m$ . The critical exponent $m=1$ at $ \textit{Ma}_\delta =590$ (black solid line with star symbols) marks a crossover from diffusion-dominated stabilisation ( $m\lt 1$ ) to advective-driven destabilisation ( $m\gt 1$ ), highlighting the offset competition as the underlying mechanism governing interfacial instability.