2. Saffman model for a dusty Kolmogorov flow

We consider an Eulerian model for a dilute suspension of heavy particles with two-way coupling introduced by Saffman (Reference Saffman1962) long ago. The model considers a dilute mono-disperse suspension of small, heavy, spherical particles with density $\rho _p$ and radius $a$ transported by a Newtonian fluid with density $\rho _f$ and viscosity $\mu$. Particle size is assumed to be much smaller than any scale in the flow such that the particle Reynolds number is negligible. The particle volume fraction $\phi _v = N_p v_p / V$, defined in terms of the volume of each particle $v_p=4{\rm \pi} a^3/3$ and the number of particles $N_p$ contained in the total volume $V$, is assumed to be negligible while the mass fraction $\phi =\phi _v \rho _p/\rho _f$ can be of order one since it is assumed $\rho _p \gg \rho _f$ (as in a dilute suspensions of water droplets in air).

Within the model, the fluid density field remains constant because of the assumption of vanishing $\phi _v$ and it is transported by the incompressible velocity field of the fluid phase ${\boldsymbol u}(\boldsymbol x,t)$. The solid phase is described by a number density field $\theta (\boldsymbol x,t) = n(\boldsymbol x,t) /(N_p/V)$, where $n(\boldsymbol x,t)$ is the local number of particles per unit volume. The normalization gives $\langle \theta \rangle = 1$, where the brackets $\langle [\cdot ] \rangle$ denote the average over the volume $V$. The number density field $\theta$ is transported by a compressible particle velocity field ${\boldsymbol v}(\boldsymbol x,t)$.

For small volume fractions ($\phi _v < 10^{-3}$) the dynamics of the particle-laden flow can be described by a two-way coupling, which takes into account the interactions between individual particles and the surrounding flow, but neglects the interactions between particles (collisions and friction) and the particle–fluid–particle interactions (fluid streamlines compressed between particles) (Elghobashi Reference Elghobashi1994). In the two-way coupling regime, the exchange of momentum between the two phases can no longer be neglected (Balachandar & Eaton Reference Balachandar and Eaton2010). For small heavy particles, such an exchange is mainly mediated by the viscous drag force which is proportional to the difference between particle and fluid velocities.

The accurate modelling of the coupling between the particles and the flow is a challenging task. In Lagrangian–Eulerian approaches based on the point-particle method, it requires taking into account the local perturbation to the fluid due to the presence of the particle (Horwitz & Mani Reference Horwitz and Mani2016). The Eulerian model proposed by Saffman is based on a simplified assumption, namely that the coupling is obtained by imposing the local conservation of the total momentum of the fluid and particle phases. This leads to the following equations (Saffman Reference Saffman1962):

(2.1)\begin{gather} \partial_t {\boldsymbol u} + {\boldsymbol u}\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol u} ={-} {\boldsymbol \nabla}p + \nu \nabla^2 {\boldsymbol u} + {\boldsymbol f} + \frac{\phi}{\tau}\theta ({\boldsymbol v} - {\boldsymbol u}), \end{gather}

(2.2)\begin{gather}\partial_t {\boldsymbol v} + {\boldsymbol v}\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol v} ={-} \frac{1}{\tau} ({\boldsymbol v} - {\boldsymbol u}), \end{gather}

(2.3)\begin{gather}\partial_t \theta + {\boldsymbol \nabla}\boldsymbol{\cdot} \left({\boldsymbol v}\theta\right) = 0, \end{gather}

where $\tau = (2/9) a^2 \rho _p/(\rho _f \nu )$ is the relaxation time of the particles, $\nu =\mu /\rho _f$ is the kinematic viscosity, p is the pressure and ${\boldsymbol f}$ is an external forcing.

In the limit of very tiny particles, i.e. small $\tau$, the Saffman model reduces to the Navier–Stokes equation for an incompressible flow with an increased density, and thus a smaller viscosity (Saffman Reference Saffman1962). Indeed, when $\tau \to 0$ from (2.2) one has ${\boldsymbol v}={\boldsymbol u}$. For small $\tau$ one can expand ${\boldsymbol v}={\boldsymbol u}+\tau \tilde {\boldsymbol v}+O(\tau ^2)$ and (2.2) gives, at leading order, $\tilde {\boldsymbol v}=-(\partial _t {\boldsymbol u}+{\boldsymbol u}\boldsymbol {\cdot }{\boldsymbol \nabla }{\boldsymbol u})+O(\tau )$. Substituting now ${\boldsymbol v}={\boldsymbol u} -\tau (\partial _t {\boldsymbol u}+{\boldsymbol u}\boldsymbol {\cdot }{\boldsymbol \nabla }{\boldsymbol u})$ in (2.3) one obtains that the particle density field remains constant at leading order. Finally, using $\theta =1+O(\tau )$ and $({\boldsymbol v}-{\boldsymbol u})/\tau = -(\partial _t {\boldsymbol u}+{\boldsymbol u}\boldsymbol {\cdot }{\boldsymbol \nabla }{\boldsymbol u}) + O(\tau )$ in (2.1) gives

(2.4)\begin{equation} \partial_t {\boldsymbol u} + {\boldsymbol u}\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol u} ={-} {\boldsymbol \nabla}p + \frac{\nu}{1+\phi} \nabla^2 {\boldsymbol u} + \frac{{\boldsymbol f}}{1+\phi}, \end{equation}

i.e. the Navier–Stokes equation for an incompressible velocity field with forcing and viscosity rescaled by the factor $(1+\phi )$.

Remarkably, we show that the same result is also recovered in the limit of large $\phi$. Indeed, from (2.1) one can write

(2.5)\begin{equation} {\boldsymbol u} = {\boldsymbol v} + \frac{\tau}{\phi} \frac{1}{\theta} \left(-\partial_t {\boldsymbol u} + {\boldsymbol u} \boldsymbol{\cdot} {\boldsymbol \nabla u} - {\boldsymbol \nabla} p + \nu \nabla^2 {\boldsymbol u} + {\boldsymbol f} \right), \end{equation}

showing that the difference between ${\boldsymbol u}$ and ${\boldsymbol v}$ is of order $1/\phi \ll 1$. Substituting ${\boldsymbol v}={\boldsymbol u} + O(1/\phi )$ in (2.3) implies $\theta =1+O(1/\phi )$ which, together with (2.5) in (2.2), yields (2.4) multiplied by $(1+\phi )$, i.e. again a Navier–Stokes equation with rescaled forcing and viscosity. We remark that the limit of large $\phi$ is physically questionable since it could violate the assumption of negligible volume fraction. Nonetheless, it is mathematically well defined and we will use it the following to discuss our results.

We consider the case of a monochromatic periodic forcing ${\boldsymbol f} = F \cos (K y) \hat {\boldsymbol x}$ which produces the Kolmogorov laminar fixed point ${\boldsymbol u}({\boldsymbol x})={\boldsymbol v}({\boldsymbol x})={\boldsymbol U}(y) \equiv U \cos (K y)$ with $U=F/(\nu K^2)$ and $\theta ({\boldsymbol x})=1$. We remark that, in general, the Kolmogorov flow is a stationary solution also for $\theta ({\boldsymbol x})=g(y)$ with $g$ an arbitrary function. Nonetheless, the solution with uniform density $\theta$ is physically the more relevant as it survives in the presence of an arbitrarily small diffusivity.

The non-dimensional parameters of the model are the Reynolds number $Re=U/(\nu K)$, defined in terms of the amplitude of the laminar flow $U$ and on the only characteristic length of the flow $K^{-1}$, the Stokes number $St=\tau \nu K^2$, defined as the ratio between the particle relaxation time $\tau$ and the viscous time $\tau _\nu =1/(\nu K^2)$ and the mass fraction $\phi$. In the following, we will study the linear stability of the laminar fixed point as a function of $Re$, $St$ and $\phi$.

We conclude this section with a comment about the limitations of the Saffman model. Beside the assumption of small volume fraction, in the case of turbulent flows at high Reynolds numbers the validity of the model (2.1)–(2.3) is in general limited to small Stokes numbers $St < 1$. This is due to the phenomenon of caustics (Wilkinson & Mehlig Reference Wilkinson and Mehlig2005) which would imply a multi-valued particle velocity field breaking the validity of the continuum description. Nonetheless, for the specific case of the linear stability of a laminar parallel flow considered here, in the laminar fixed point the particle velocity field is equal to the fluid velocity field, and therefore the model is well defined for arbitrary value of $St$.

3. Linear stability analysis

We study the linear stability of an infinitesimal perturbation of the basic Kolmogorov flow. To this aim we expand (2.1)–(2.3) around the laminar fixed point ${\boldsymbol u}({\boldsymbol x},t)={\boldsymbol U}(y)+{\boldsymbol u}'({\boldsymbol x},t)$, ${\boldsymbol v}({\boldsymbol x},t)={\boldsymbol U}(y)+{\boldsymbol v}'({\boldsymbol x},t)$, $\theta ({\boldsymbol x},t)=1+\theta '({\boldsymbol x},t)$, and obtain the linearized equations for the perturbations

(3.1)\begin{gather} \partial_t {\boldsymbol u}' + {\boldsymbol U}\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol u}' + {\boldsymbol u}'\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol U}={-} {\boldsymbol \nabla}p' + \nu \nabla^2 {\boldsymbol u}' + \frac{\phi}{\tau}({\boldsymbol v}' - {\boldsymbol u}'), \end{gather}

(3.2)\begin{gather}\partial_t {\boldsymbol v}' + {\boldsymbol U}\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol v}' + {\boldsymbol v}'\boldsymbol{\cdot}{\boldsymbol \nabla}{\boldsymbol U} ={-} \frac{1}{\tau}({\boldsymbol v}' - {\boldsymbol u}'), \end{gather}

(3.3)\begin{gather}\partial_t \theta' + {\boldsymbol U}\boldsymbol{\cdot}{\boldsymbol \nabla}\theta'+ {\boldsymbol \nabla}\boldsymbol{\cdot} {\boldsymbol v}'=0. \end{gather}

We observe that, at this order, the density field becomes a passive scalar since it does not enter (3.1)–(3.2). Therefore, the evolution of $\theta '$ can be neglected.

A remarkable simplification of the linear stability analysis can be achieved by invoking the Squire theorem for parallel flows (Squire Reference Squire1933), which states that it suffices to consider two-dimensional perturbations, since three-dimensional perturbations are more stable. From the original formulation, the theorem has been extended to various systems, including magnetohydrodynamic equations (Hughes & Tobias Reference Hughes and Tobias2001), stratified flows (Balmforth & Young Reference Balmforth and Young2002) and viscoelastic flows (Bistagnino et al. Reference Bistagnino, Boffetta, Celani, Mazzino, Puliafito and Vergassola2007). In Appendix A, we report the derivation of the Squire theorem for the dusty fluid model (3.1)–(3.2).

In the following we will therefore consider the two-dimensional version of the linearized equation. It is convenient to rewrite the fluid velocity fluctuation in terms of a streamfunction ${\boldsymbol u}'=(\partial _y \varPsi,-\partial _x \varPsi )$ and the compressible particle velocity in terms of a particle streamfunction $\varPsi _p$ and potential $\varPhi _p$ as ${\boldsymbol v}'=(\partial _y \varPsi _p+\partial _x \varPhi _p, -\partial _x \varPsi _p+\partial _y \varPhi _p)$. In terms of these fields the linear equations (3.1)–(3.2) read

(3.4)\begin{gather} \partial_t \nabla^2 \varPsi +U \cos(K y) (K^2+\nabla^2)\partial_x \varPsi -\nu \nabla^4 \varPsi+ \frac{\phi}{\tau} \nabla^2 (\varPsi -\varPsi_p) =0, \end{gather}

(3.5)\begin{gather}\partial_t \nabla^2 \varPsi_p \!+\!U \cos(K y) \left[(K^2\!+\!\nabla^2)\partial_x \varPsi_p -K^2 \partial_y \varPhi_p \right]\!-\!U k \sin(K y) \nabla^2 \varPhi_p + \frac{\nabla^2 (\varPsi_p - \varPsi)}{\tau} =0, \end{gather}

(3.6)\begin{gather}\partial_t \nabla^2 \varPhi_p +U \cos(K y) \partial_x \nabla^2\varPhi_p- 2 U K \sin(K y) \left(\partial_x \partial_y \varPhi_p - \partial_x^2 \varPsi_p \right)+ \frac{1}{\tau} \nabla^2 \varPhi_p =0. \end{gather}

For a Newtonian fluid ($\phi =0$), the laminar solution is known to be linearly stable to perturbations at wavenumbers larger than $K$, and to become unstable to large-scale transverse perturbations (i.e. in the direction $x$ transverse to the direction of modulation $z$) above the critical value $Re_c=\sqrt {2}$ (Meshalkin & Sinai Reference Meshalkin and Sinai1961; Sivashinsky & Yakhot Reference Sivashinsky and Yakhot1985). As discussed in § 2, in the limit of small inertia ($\tau \ll 1$) or large mass fraction ($\phi \gg 1$) the Saffman model recovers the Navier–Stokes equation with a rescaled viscosity $\nu /(1+\phi )$. Therefore, in these limits we expect the critical Reynolds number to become $Re_c = \sqrt {2}/(1+\phi )$, i.e. the presence of tiny particles, or a large mass fraction of particles, makes the flow more unstable.

4. Multiple-scale analysis

The general dependence of the critical Reynolds number on the parameters $\tau$ and $\phi$ can be obtained by a standard multiple-scale analysis (Bensoussan et al. Reference Bensoussan, Lions and Papanicolaou2011) of the linearized equations (3.4)–(3.6). The main idea of the multiple-scale method is to search for a perturbation which varies on spatial scales much larger than those of the base flow. For this purpose, beside the small-scale variables $x$, $y$ and $t$, the multiple-scale method introduces the large-scale spatial variables $X=\varepsilon x$, $Y=\varepsilon y$ and a corresponding slow time $T=\varepsilon ^2 t$, where the small parameter $\varepsilon$ is the ratio between the characteristic scales of the basic flow and the perturbation. The relative powers of $\varepsilon$ in the space and time variables reflect the diffusive dynamics expected at large scales. The two sets of variables are then assumed to be independent, so that by averaging over the small scales it is possible to obtain an effective diffusion-like equation for the large scales, which defines an eddy viscosity. A change of sign of the eddy viscosity corresponds to a change of the stability of the perturbation. In particular, the system becomes unstable when the eddy viscosity becomes negative (Sivashinsky & Yakhot Reference Sivashinsky and Yakhot1985; Dubrulle & Frisch Reference Dubrulle and Frisch1991).

The choice of the multiple-scale method to study the stability of the dusty Kolmogorov flow is motivated by the fact that, in the Newtonian case (at $\phi =0$), the most unstable perturbation is indeed at large scale, and the multiple-scale prediction for the critical Reynolds number is correct. For simplicity of the calculation, and in analogy with the Newtonian case, we also assume that the most unstable perturbation is transverse, i.e. depends on the large-scale variable $X$ only and not on $Y$. The validity of these assumptions for the dusty gas at $\phi > 0$ will be checked by extensive numerical simulations of the linear systems in § 5. In particular, we anticipate that, while the transverse nature of the most unstable perturbation was always confirmed, in a certain parameter region the scale separation appears to be violated, when this happens the multiple-scale approach does not provide the correct prediction.

Before proceeding, it is convenient to rewrite (3.4) in terms of a co-streamfunction defined as $\chi =\varPsi +\phi \varPsi _p$. Such a choice allows us to remove the apparent singularity of the Stokes drag in the limit $\tau \to 0$ in (3.4). Indeed, by combining (3.4) and (3.5) we obtain

(4.1)$$\begin{gather} \partial_t \nabla^2 \chi + U \cos(K y) (K^2+\nabla^2)\partial_x \chi- \phi U K \left(K \cos(K y) \partial_y \varPhi_p +\sin(K y) \nabla^2 \varPhi_p \right)\nonumber\\ - \nu \nabla^4 (\chi-\phi \varPsi_p) = 0 , \end{gather}$$

which removes the explicit dependence of (3.4) on $\tau$. The linear systems for the perturbative analysis is hence formed by the set of equations (4.1), (3.5) and (3.6).

Following the multiple-scale method, we assume a perturbative expansion of the fields

(4.2)\begin{equation} \left. \begin{aligned} \chi(X,y,T) & = \chi_0(X,y,T)+\varepsilon \chi_1(X,y,T)+\varepsilon^2 \chi_2(X,y,T), \\ \varPsi_p(X,y,T) & = \varPsi_{p,0}(X,y,T)+\varepsilon \varPsi_{p,1}(X,y,T)+\varepsilon^2 \varPsi_{p,2}(X,y,T),\\ \varPhi_p(X,y,T) & = \varPhi_{p,0}(X,y,T)+\varepsilon \varPhi_{p,1}(X,y,T)+\varepsilon^2 \varPhi_{p,2}(X,y,T) . \end{aligned} \right\} \end{equation}

The derivative operators are transformed as $\partial _x \to \varepsilon \partial _X$, $\partial _t \to \varepsilon ^2 \partial T$. Notice that the base flow does not depend on $x$ and $t$ and therefore the same holds for the perturbation. By inserting the expansions (4.2) and the fast/slow variable decomposition into (4.1), (3.5) and (3.6) we obtain, at order $\varepsilon ^0$, that the zero-order fields do not depend on the fast variable, i.e. $\chi _0(X,y,T)=a_0(X,T)$, $\varPsi _{p,0}(X,y,T)=b_0(X,T)$ and $\varPhi _{p,0}(X,y,T)=c_0(X,T)$. At the order $\varepsilon ^2$, the absence of secular terms requires $c_0=0$.

The solvability condition is obtained by integrating (3.5)–(4.1) over one period of the fast variable $y$. The first non-trivial condition is obtained at order $\varepsilon ^3$ and gives a relation among the large-scale fields

(4.3)\begin{equation} a_0(X,T)=(1+\phi) b_0(X,T) . \end{equation}

At order $\varepsilon ^4$, we finally get the diffusion equation for the slow field $a_0$

(4.4)\begin{equation} \frac{2}{Re} (1+\phi) \partial_T \partial_X^2 a_0(X,T) + \nu Re \left(1 - \frac{2}{ Re^2} + 2 \phi + \phi^2 - \phi St \right) \partial_X^4 a_0(X,T) = 0, \end{equation}

which defines the eddy viscosity

(4.5)\begin{equation} \nu_e = \nu \frac{Re^2}{2(1+\phi)} \left(\frac{2}{Re^2} -(1+\phi)^2+ \phi St \right). \end{equation}

The critical Reynolds number is finally obtained by the condition $\nu _e=0$ at which the eddy viscosity becomes negative, indicating that the basic flow is linearly unstable (Sivashinsky & Yakhot Reference Sivashinsky and Yakhot1985; Dubrulle & Frisch Reference Dubrulle and Frisch1991)

(4.6)\begin{equation} Re_c = \sqrt{\frac{2}{(1+\phi)^2 - \phi St}}. \end{equation}

For $\phi =0$ the critical Reynolds number predicted by (4.6) recovers correctly the Newtonian value $Re_c=\sqrt {2}$. Interestingly, the same value is recovered also on the neutral curve $St=2+\phi$. For $St \to 0$ we obtain the Saffman limit $R_c=\sqrt {2}/(1+\phi )$. For $St<2$, (4.6) predicts a monotonic decrease of $Re_c$ as a function of $\phi$ ($Re_c(\phi ) \le Re_c(0)$), indicating that tiny particles always destabilize the flow. On the contrary, for $St>2$, the critical $Re$ depends on $\phi$ in a non-monotonic way: for $\phi < \phi _{max}=(St-2)/2$, $Re_c$ increases monotonically above the Newtonian value $\sqrt {2}$, it reaches a maximum at $\phi _{max}$ after which it monotonically decreases and for $\phi >St-2$ it goes below $\sqrt {2}$. Therefore, increasing the mass fraction particles first stabilize the flow up to a maximum then the stabilizing effect decreases and, finally, for large enough mass fraction particles make the flow more unstable than the in the Newtonian case. We remark that, according to (4.6), the dusty Kolmogorov flow should always be stable (i.e. $Re_c \to \infty$) for $St \ge (1+\phi )^2/\phi \ge 4$. We will see that this is actually an overestimation of the stability due to the fact that the main assumption of the multiple-scale analysis (instability to large-scale perturbations) does not hold in a certain region of the parameter space ($\phi$, $St$).

5. Numerical analysis

To check the validity of the analytical result (4.6) obtained with the multiple-scale analysis, we performed an extensive numerical study of the linearized equations in two dimensions (3.4)–(3.6) by means of a pseudo-spectral method in a square domain of size $L =2{\rm \pi}$ with periodic boundary conditions. For each set of values of the parameters $\phi$ and $St$ in the range $0 \le \phi \le 6$ and $0 \le St \le 6$ we have studied the stability of the system on varying the $Re$ number. The latter has been varied by changing the amplitude of the forcing $F$ while keeping fixed the viscosity $\nu =10^{-3}$ and the scale of the base flow $1/K$. Simulations were done at two different resolutions, with $128^2$ and $256^2$ grid points and forcing wavenumber $K=32$ and $K=64$, respectively, to check finite size effects. A random initial perturbation has been imposed on each Fourier mode $(k_x,k_y)$ in the range $0 \le |{\boldsymbol k}| \le K$. The stability of each mode and its growth rate is determined by the temporal evolution of its amplitude after a short transient. The critical Reynolds number was determined by means of the bisection method based on the total kinetic energy of the fluid.

In figure 1(a) we plot the critical Reynolds number as a function of the Stokes number for different mass fraction values $\phi$. At small $St$, $Re_c$ is smaller than that of the single-phase fluid ($Re_c=\sqrt {2}$ for $\phi =0$) and the numerical results are in agreement with the theoretical prediction (4.6). In particular, in the limit $St \ll 1$, the critical Reynolds number recovers the Saffman limit $Re_c = \sqrt {2}/(1+\phi )$ (not shown). The critical Reynolds number increases monotonically with $St$ at fixed $\phi$, eventually becoming larger than $\sqrt {2}$, meaning that large particle inertia has a stabilizing effect on the flow. At a qualitative level, the physical mechanisms of the stabilizing/destabilizing effect of the particles have been already discussed by Saffman (Reference Saffman1962). Particles with small $St$ follow the flow almost like tracers, so that their effect is simply to increase the density of the suspension. Therefore, the dusty gas behaves as a Newtonian flow with a reduced kinematic viscosity (see (2.4)) which makes the flow more unstable. Conversely, particles with large inertia do not follow the perturbation of the flow, but they ‘carry on with the velocity of the base flow’ (Saffman Reference Saffman1962). The disturbance has therefore to flow around the particles, dissipating its energy because of the viscous drag. Our numerical results show that the stabilizing effect at large $St$ is weaker than the prediction (4.6). In particular, the multiple-scale method predicts unconditioned stability (i.e. $Re_c =\infty$) for $St \ge (1+\phi )^2/\phi$, while in the numerical simulations $Re_c$ remains finite.

Figure 1. Critical Reynolds number (a) as a function of the Stokes number $St$ for different values of $\phi$ and (b) as a function of the mass fraction $\phi$ for different values of $St$. The values of the parameters $\phi$ and $St$ are reported in the legend. Solid curves denote the multiscale prediction (4.6); symbols the numerical results; dashed lines display the Newtonian value $\sqrt {2}$; dash-dotted line in panel (b) shows the Saffman limit $\sqrt {2}/(1+\phi )$.

The behaviour of $Re_c$ as a function of the mass fraction $\phi$ for fixed values of $St$, shown in figure 1(b), gives further insights into the stability of the system. In agreement with (4.6), we find that $Re_c$ is monotonically decreasing for $St \le 2$ (i.e. tiny particles always destabilize the flow). Conversely, at $St>2$ the particles at low concentration stabilize the flow while at sufficiently large concentrations $\phi \ge St - 2$ they have a destabilizing effect. It is interesting to note that a similar non-monotonic behaviour as a function of the mass loading has been observed also for the skin-friction coefficient in Lagrangian–Eulerian simulations of inertial particles in a vertical channel flow (Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2018). The physical mechanism of the destabilizing effect at large $\phi$ is similar to that of the case of small $St$. The strong drag exerted by the large mass fraction forces the fluid to follow closely the particle velocity (see (2.5)). As a consequence, the dusty gas behaves almost as a single-phase fluid with a larger density and therefore a smaller kinematic viscosity, which reduces its stability. From figure 1(b) it is evident that the agreement between the multiple-scale result and numerical simulations is very good for any $\phi$ up to $St=3$. For $St \ge 4$, the multiple-scale result (4.6) overestimates the $Re_c$ for an intermediate interval of values of $\phi$ around $\phi _{max} = (St-2)/2$. Nonetheless, also in these cases ($St=4$ and $St=5$) the multiple-scale prediction works well for small and large values of $\phi$.

In order to understand why the multiple-scale analysis fails in predicting the correct $Re_c$ at large $St$, we computed numerically the growth rate $\sigma$ as a function of the wavenumber $k_x$ of the perturbation (i.e. the dispersion relation) for different values of the parameters $\phi$ and $St$. In figure 2(a) we show the dispersion relation computed at the critical point $Re=Re_c$ for $St=4$ and three values of the mass fraction. For small and large mass fractions ($\phi =0.1$ and $\phi =5$) we observe that the growth rate $\sigma$ is a monotonically decreasing function of $k_x$ and the unstable mode is the smallest available wavenumber $k_c = k_{min} \equiv 2{\rm \pi} /L$. In these cases, the hypothesis of large scale separation is justified and indeed the predictions of the multiple-scale analysis are in agreement with the numerical results. Conversely, for $\phi =1$ the curve $\sigma (k_x)$ is non-monotonic and the unstable mode appears to be at $k_c \simeq 0.3 K$, therefore the instability is no more triggered by large-scale perturbations and multiple-scale analysis fails to predict the instability. Similar behaviours have been observed also for $St=5$ and $St=6$ (not shown). To systematically investigate the region of parameters for which the multiple-scale analysis is not expected to work we numerically studied the dependence of the unstable (transverse) mode $k_c$ on $\phi$ and $St$ at $Re=Re_c$, shown in figure 2(b). For $St \ge 4$ and intermediate values of $\phi$ we found $k_c \simeq 0.3 K$, while for $St \le 3$ (not shown) we always found $k_c=k_{min}$ in agreement with the multiple-scale assumption. By comparing figures 2(b) and 1(b) we clearly observe the correspondence between the theoretical–numerical agreement in figure 1(b) and the fact that $k_c \ll K$.

Figure 2. Panel (a) shows the non-dimensional growth rate $\sigma /(\nu K^2)$ as a function of the $x$-wavenumber $k_x$ at the onset of the instability $Re=Re_c$ for $St=4$ and different values of the mass fraction $\phi$ as labelled. Panel (b) shows the first unstable mode $k_c$ normalized by $K$ as a function of $\phi$ for different values of $St$ as labelled.