2.1. Equilibrium statistical mechanics

Interfaces separate coexisting phases. In equilibrium, all thermodynamic properties of interfaces follow from the intensive interfacial free energy $f_I$ , which is precisely the interfacial tension $\sigma _{\alpha, \beta }$ between coexisting bulk phases $\alpha$ and $\beta$ . For pedagogical reasons, we illustrate this by considering the ferromagnetic, spin-half, Ising model, which is used to model magnets that are anisotropic (in spin space). This Ising model, which can be mapped onto a simple lattice-gas model for a fluid that can have two coexisting bulk phases (see below), has the Hamiltonian

(2.1) \begin{equation} \mathcal {H} = - J \sum _{\langle i,j\rangle } S_iS_j - h \sum _i S_i . \end{equation}

The spins $S_i = \pm 1$ occupy the $N \sim L^d$ sites $i$ of a $d$ -dimensional hypercubic lattice, with side $L$ , in a region $\Omega$ ; the exchange interaction $J$ couples spins on nearest-neighbour pairs of sites $\langle i,j\rangle$ ; a ferromagnetic interaction $(J \gt 0)$ favours spin alignment, i.e. spins of the same sign on nearest-neighbour sites; the external magnetic field $h$ favours $S_i = sgn (h)$ and it distinguishes between the two phases: up-spin $\uparrow$ (with $S_i = +1, \forall i$ , at zero temperature $T=0$ ) and down-spin $\downarrow$ (with $S_i = -1, \forall i$ , at $T=0$ ). If we are interested only in bulk statistical physics, it is convenient to use periodic boundary conditions (PBCs) in all $d$ directions.

The transformation $n_i = (1+S_i)/2$ , yields the lattice gas variables $n_i = 0 \, \textrm {or} \, 1$ and defines the lattice-gas model, whose low- and high-density phases are the analogues of the $\downarrow$ and $\uparrow$ phases in the ferromagnetic Ising model and in which the chemical potential is the counterpart of $h$ in the Ising model (see, e.g., Griffiths Reference Griffiths1972; Goldenfeld Reference Goldenfeld2018). (A similar transformation relates the Ising-model Hamiltonian to a lattice model for a binary alloy.) Just as a positive (negative) external field $h$ favours the up-spin $\uparrow$ (down-spin $\downarrow$ ) phase, a large and positive (negative) values of the chemical potential favour the high-density (low-density) phase in the lattice-gas model (see table 1).

Table 1.

Correspondences between Ising ferromagnetic systems, with lattice-gas or binary-mixture counterparts, and their phases and phase transitions; here, $\uparrow$ and $\downarrow$ represent the up-spin and down-spin phases of an Ising ferromagnet (see the text). Rows 2–4 mention equilibrium phases and transitions; the last row mentions non-equilibrium active fluids which can exhibit active phase-separation (see § 7.6).

The intensive bulk free energy $f_B$ of the Ising model, which is a function of the temperature $T$ , $h$ and $J$ , is defined as follows:

(2.2) \begin{eqnarray} f_B(T,h,J) &=& \lim _{N\to \infty } \frac {1}{N}\left [F(T,h,J,\Omega )\right ]; \nonumber \\ F(T,h,J,\Omega )] &=& -k_BT\ln [\mathcal {Z}] ; \nonumber \\ {\mathcal {Z}} &=& \sum _{\{S_i = \pm 1\}} \exp [-\mathcal {H}/(k_BT)]; \end{eqnarray}

where, $k_B,\,\, F$ and $\mathcal Z$ are, respectively, the Boltzmann constant, the total free energy and the partition function and the sum is over all the spin states. Note that $F$ depends on $\Omega$ and, therefore, on BCs; by contrast, $f_B$ does not depend on $\Omega$ , the BCs, and boundary couplings and fields because we have taken the thermodynamic limit $N \to \infty$ . If the thermodynamic limit exists (as it does for the Ising model we consider), then: (i) $f_B$ is a convex-up function of its arguments, which guarantees continuity of $f_B$ with respect to $T$ , $h$ and $J$ ; (ii) $\partial f_B/\partial T, \partial f_B/\partial h$ and $\partial f/\partial J$ exist almost everywhere ( $f_B$ has slope discontinuities at first-order phase boundaries); and (iii) $\partial f_B/\partial T,\, \partial f_B/\partial h$ and $\partial f/\partial J$ are, respectively, monotone, non-increasing functions of $T,\, h$ and $J$ . Bulk thermodynamic functions follow from derivatives of $f_B$ ; e.g. the magnetisation per site

(2.3) \begin{equation} M = \sum _i \langle S_i \rangle /N = -\partial f_B/ \partial h , \end{equation}

where the angular brackets denote the thermal average. When the second derivatives of $f_B$ exist, they lead to stability conditions such as the positivity of the specific heat at constant $h$ , to wit, $C_h \equiv -k_BT\partial ^2f_B/\partial T^2 \geqslant 0$ .

This magnetisation $M \in [-1,1]$ is the order parameter, which is positive (negative) in the up-spin (down-spin) phase; its lattice-gas counterpart is the density $\rho = \sum _i \langle n_i \rangle /N \in [0,1]$ . For dimension $d \gt 1$ , the Ising-model phase diagram in the $h/J-k_BT/J$ plane (henceforth, we set $k_B=1$ and $J=1$ ) shows a first-order phase boundary, $h=0$ for $T \lt T_c$ , where $T_c$ is the critical or Curie temperature at which this model exhibits a continuous phase transition from the ferromagnetic to the paramagnetic phase (figure 2a ). Along the first-order boundary, the $\uparrow$ and $\downarrow$ phases coexist. In the $T{-}M$ counterpart (figure 2 b) of the $h{-}T$ phase diagram, two-phase coexistence occurs everywhere below the coexistence curve in the peach-shaded region: the system undergoes macroscopic phase separation, forming distinct $\uparrow$ and $\downarrow$ phases, rather than a fully mixed state; in equilibrium, these coexisting phases are separated by an interface (see below). If we use the Ising-model–lattice-gas correspondence (table 1), we see that its binary-mixture counterpart is the coexistence of immiscible fluids, below the critical temperature above which the two fluids mix.

We do not concentrate on the critical properties of the Ising model here; but we do mention critical-point behaviours occasionally. Therefore, we note that, in the vicinity of the critical point at $h=0$ and $T=T_c$ , thermodynamic functions display power-law forms, with universal exponents and scaling functions (see, e.g., Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Kardar Reference Kardar2007; Goldenfeld Reference Goldenfeld2018). For example, at $h=0$ and $(T-T_c)/T_C \ll 1$ , $M \sim [(T_c-T)/T_c]^\beta$ , for $T \lesssim T_c$ ; and the specific heat $C_h \sim [(T_c-T)/T_c]^{-\alpha }$ diverges; $\beta$ and $\alpha$ are universal critical exponents that depend on $d$ and the number of components $n$ of the order parameter (for the Ising model we have a scalar order parameter and $n=1$ ). If we were to obtain order-parameter correlation functions and, therefrom, the correlation length $\xi$ , then we would find the divergence $\xi \sim [(T_c-T)/T_c]^{-\nu }$ . At such a critical point, only two exponents are independent; the others can be related by scaling laws that are linear relations between these exponents (see, e.g., Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Kardar Reference Kardar2007; Goldenfeld Reference Goldenfeld2018). For the $d=2$ Ising model, $\alpha = 0$ because $C_h$ has a logarithmic divergence at $h=0$ and $T=T_c$ , $\beta = 1/8$ , and $\nu = 1$ .

Figure 2.

Schematic phase diagrams for the $d$ -dimensional Ising model (2.1) in the magnetic-field $h$ and temperature $T$ plane (a) and the $T$ and magnetisation $M$ plane (b). There is a first-order phase boundary (brown line), $h=0$ for $T \lt T_c$ , that ends in a critical point at $h=0,\, T=T_c$ (blue point); along the first-order boundary, the $\uparrow$ and $\downarrow$ phases coexist. In the $T{-}M$ phase diagram, two-phase coexistence occurs in the peach-shaded region everywhere below the coexistence curve (black), atop which is the critical point.

To define the intensive interfacial free energy $f_I$ (per unit area of the interface), we must distinguish between BCs, for the Ising Hamiltonian (2.1), that yield a $(d-1)-$ dimensional interface between coexisting bulk phases and those that do not. We illustrate two such BCs, $++$ and $-+$ , via the schematic diagrams for $d=2$ in figure 3. We use periodic BCs in $(d-1)$ directions; in the remaining direction, denoted by $x$ , we have two $(d-1)$ -dimensional surfaces at $x=-L/2$ and $x=+L/2$ ; the $++$ BC does not yield an interface; by contrast, the $-+$ BC yields an interface with $N_I \sim L^{(d-1)}$ the number of sites in the $T=0$ interface. (The precise value of $x$ at the interface location depends on whether we work in the fixed- $h$ or fixed- $M$ ensemble; in the former, the interface can lie anywhere between $x=-L/2$ and $x=+L/2$ , whereas, in the latter, the position of the interface is such that the proportion of the $\uparrow$ and $\downarrow$ regions yields the fixed magnetisation $M$ .) Here,

(2.4) \begin{equation} f_I(T,h,J) = \lim _{N_I \to \infty ; N\to \infty }\frac {1}{N_I}\left [F^{-+} - F^{++}\right ], \end{equation}

where $F^{++}$ and $F^{-+}$ are, respectively, the total free energies of the Ising model with $++$ and $+-$ BCs (figure 3). (Note that the free-energy contributions from the two surfaces cancel when we subtract $F^{++}$ from $F^{-+}$ in (2.4).) We have mentioned above that bulk statistical mechanics and phase diagrams follows from $f_B$ and its derivatives; similarly, interfacial statistical mechanics and interfacial phase diagrams follow from $f_I$ and its derivatives. If the order parameter has more than one component, e.g. if $n=2$ as in the XY model or a superfluid, then $f_I$ has to be replaced by the helicity modulus $\Upsilon$ , because the interface is not sharp (see, e.g., Fisher, Barber & Jasnow Reference Fisher, Barber and Jasnow1973); in such cases, $-+$ BCs lead to a correction to the total free energy that scales as $L^{(d-2)}$ . We restrict ourselves to scalar order parameters that have $n=1$ . Furthermore, $f_I$ gives us the strict definition of the interfacial (or surface) tension (see, e.g., Rottman & Wortis Reference Rottman and Wortis1984). Both $f_I$ and surface tension have the same physical units: $f_I$ is typically given in energy per unit area, whereas the surface tension is specified as a force per unit length, both of which are dimensionally equivalent. In particular, if $T \to T_c$ from below, the coexisting phases come together at the critical point (the consolute point in a binary-fluid mixture) and the interfacial free energy vanishes as $f_I \sim [(T_c-T)/T_c]^{\mu }$ , with the exponent $\mu = 2 \nu$ .

Figure 3.

Schematic diagrams for (a) $++$ and (b) $-+$ BCs for the $d$ -dimensional Ising model (2.1); $d=2$ in this illustration. In $d-1$ directions ( $y$ here) we use periodic BCs.

It is important to distinguish between surface and interfacial free energies. The former is associated with a geometrically imposed surface (say on the (1,0,0, …) face of the hypercubic lattice we consider). The intensive surface free energy is (see, e.g., Fisher & Caginalp Reference Fisher and Caginalp1977; Caginalp & Fisher Reference Caginalp and Fisher1979; Rottman & Wortis Reference Rottman and Wortis1984)

(2.5) \begin{equation} f_s(T,h,J) = \lim _{N_s \to \infty ; N\to \infty } \frac {1}{2N_s}\left [F^{++} - Nf_B\right ], \end{equation}

where $N_s \sim L^{(d-1)}$ is the number of sites in the geometrically imposed surface (or surfaces); for the $++$ BC in figure 3, there are two surfaces at $x=-L/2$ and $x=+L/2$ , so we have a factor of $2$ in (2.5). Geometrically imposed surfaces do not fluctuate, but interfaces do, because of thermal or other (e.g. turbulent) fluctuations.

It is well known that the $d$ -dimensional Ising model can be solved exactly (i.e. $f_B$ can be obtained analytically) only if (i) $d=1$ or (i) $d=2$ and $h=0$ (see, e.g., Onsager Reference Onsager1944; Huang Reference Huang2008; Thompson Reference Thompson2015; Baxter Reference Baxter2016). In other dimensions, we must either use approximations, such as mean-field theory, or numerical simulations, such as the Monte Carlo method (see, e.g., Plischke & Bergersen Reference Plischke and Bergersen1994; Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Gould et al. Reference Gould, Tobochnik, Meredith, Koonin, McKay and Christian1996; Kardar Reference Kardar2007; Goldenfeld Reference Goldenfeld2018). For $d=2$ , a variety of elegant results can be obtained for $f_I$ ; these are of relevance to equilibrium crystal shapes (see Rottman & Wortis Reference Rottman and Wortis1984).

We give a brief introduction to mean-field theory because it leads directly to our discussion of the CH system. We consider a generalisation of the Hamiltonian (2.1) in which the magnetic field $h$ is replaced by a site-dependent magnetic field $h_i$ , so the second term becomes $\sum _i h_i S_i$ . In its most rudimentary form, mean-field theory (known as Curie–Weiss theory for our Ising-model example) uses

(2.6) \begin{eqnarray} S_iS_j &=& (S_i - \langle S_i \rangle ) (S_j - \langle S_j \rangle ) + S_i \langle S_j \rangle + S_i \langle S_j\rangle - \langle S_i \rangle \langle S_j \rangle ;\nonumber \\ S_iS_j&\simeq & + S_i \langle S_j \rangle + S_i \langle S_j\rangle - \langle S_i \rangle \langle S_j \rangle . \end{eqnarray}

The second equation follows from the first one if we neglect $(S_i - \langle S_i \rangle ) (S_j - \langle S_j \rangle )$ because it is quadratic in deviations from the equilibrium values of the site magnetisations $M_i \equiv \langle S_i \rangle$ . We now substitute the second row of (2.6) into the Hamiltonian (2.1); this yields a mean-field Hamiltonian $H_{MF}$ in which a spin at site $i$ experiences an effective magnetic field $h^{eff}_i = h_i + J\sum _{j\in [nni]} M_j$ , where $[nni]$ indicates all the nearest-neighbour sites of $i$ . The single-site magnetisation for $\mathcal {H}_{MF}$ follows simply and we get the Curie–Weiss self-consistency equations

(2.7) \begin{equation} M_i = \tanh \left [\frac {(h_i + J\sum _{j\in [nni]} M_j)}{(k_BT)}\right ], \end{equation}

which can have many solutions. The variational formulation of this mean-field theory (see, e.g., Falk Reference Falk1970; Girardeau & Mazo Reference Girardeau and Mazo1973; Plischke & Bergersen Reference Plischke and Bergersen1994) implies that we must select the solution that leads to the global minimum of the variational free energy

(2.8) \begin{align} \mathcal {F}_{CW}(\{\mathcal {M}_i\}) &= -\sum _ih_i\mathcal {M}_i - J\sum _{\lt i,j\gt } (\mathcal {M}_i\mathcal {M}_j) \nonumber \\ &\quad + \frac {k_BT}{2}\sum _i\left [(1+\mathcal {M}_i)\ln \frac {(1+\mathcal {M}_i)}{2} + (1-\mathcal {M}_i)\ln \frac {(1-\mathcal {M}_i)}{2}\right ], \end{align}

where $\mathcal {M}_i$ are the variational parameters and the subscript $CW$ stands for Curie–Weiss. Note that $\mathcal {F}(\{\mathcal {M}_i\})$ is not the equilibrium free energy, so it might not be a convex function of its arguments. If we minimise $\mathcal {F}(\{\mathcal {M}_i\})$ by setting $\partial \mathcal {F}(\{\mathcal {M}_i\})/\partial \mathcal {M}_i = 0$ , we recover the self-consistency equations (2.7). At the global minimum of $\mathcal {F}(\{\mathcal {M}_i\})$ , the equilibrium value of $\mathcal {M}_i$ is $\mathcal {M}_{i,eq}$ ; when we substitute this value of $\mathcal {M}_i$ in $\mathcal {F}(\{\mathcal {M}_i\})$ , we obtain the mean-field expression for the equilibrium free energy, which satisfies all the convexity properties mentioned earlier.

If the variational parameters $\mathcal {M}_i$ vary slowly in space, say over a length scale $\ell$ , with $a/\ell \ll 1$ , where $a$ is the lattice spacing of the hypercubic lattice, then we can make the following continuum approximation:

(2.9) \begin{eqnarray} \mathcal {M}_i &\to & \phi (\mathbf {r}) ; \,\,\,\, h_i\to h(\mathbf {r}); \,\,\,\, \sum _i \to \frac {1}{a^d}\int {\textrm{d}}\mathbf {r} ; \nonumber \\ &-& J\sum _{\lt i,j\gt }(\mathcal {M}_i\mathcal {M}_j) \to \left [-\frac {qJ}{2a^d}\int {\textrm{d}}\mathbf {r}[\phi (\mathbf {r})]^2+\frac {J}{2a^{(d-2)}}\int {\textrm{d}}\mathbf {r}[\nabla \phi (\mathbf {r})]^2 \right ] \nonumber \\ &+& O((a/\ell )^4), \end{eqnarray}

where $q = 2d$ is the nearest-neighbour coordination number for the $d$ -dimensional hypercubic lattice and $\phi$ is a scalar order parameter. We now use (2.9) to obtain the continuum limit of the variational free energy functional (2.8),

(2.10) \begin{eqnarray} a^d\mathcal {F}_{CW}(\phi, \nabla \phi ) &=& \int {\textrm{d}}\mathbf {r} \left [-h(\mathbf {r}) \phi (\mathbf {r}) - \frac {qJ}{2} [\phi (\mathbf {r})]^2 + \frac {Ja^2}{2} [\nabla \phi (\mathbf {r})]^2 + \frac {k_BT}{2} V_{CW}(\phi ) \right ] , \nonumber \\ V_{CW}(\phi ) &\equiv & \left [(1+\phi )\ln \frac {(1+\phi )}{2} + (1-\phi )\ln \frac {(1-\phi )}{2}\right ], \end{eqnarray}

where, in the Curie–Weiss approximation, the first three terms in (2.10) yield the energy contribution and the fourth term the entropy contribution (see, e.g., Puri Reference Puri2004); furthermore, $\phi \in [-1,1]$ .

If we expand $V_{CW}(\phi )$ to $O(\phi ^4)$ and set $a=1$ , we obtain the Landau–Ginzburg (LG) variational free-energy functional (see, e.g., Kadanoff et al. Reference Kadanoff, Götze, Hamblen, Hecht, Lewis, Palciauskas, Rayl, Swift, Aspnes and Kane1967; Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Kardar Reference Kardar2007; Landau & Lifshitz Reference Landau and Lifshitz2013; Goldenfeld Reference Goldenfeld2018),

(2.11) \begin{eqnarray} \mathcal {F}_{LG}(\phi, \nabla \phi ) &=& \int {\textrm{d}}\mathbf {r} \left [ g(\phi ) + \frac {\sigma }{2}[\nabla \phi ]^2 \right ] , \nonumber \\ g(\phi ) &\equiv & \left [-h \phi + \frac {k_B}{2} (T - T_c^{CW}) \phi ^2 + \frac {k_BT}{12} \phi ^4 +O(\phi ^6)\right ] , \end{eqnarray}

where we suppress the $\mathbf {r}$ argument of $h$ and $\phi$ , the Curie–Weiss critical temperature $T_c^{CW} \equiv qJ$ , and $\sigma \equiv J$ is the bare surface-tension cost for large gradients in the scalar-order-parameter field $\phi$ . If we consider uniform ordering, then minimisation of $g$ via $\partial g(\phi )/\partial \phi = 0$ yields the mean-field-theory phases, transitions, and exponents for this Ising model. We present a schematic plot of $g(\phi )$ in figure 4(a) for $T \lt T_c^{CW}$ and $h=0$ . Note that $g$ displays two, equally deep quadratic minima at

(2.12) \begin{equation} \phi _{b} \equiv \sqrt {\frac {3(T_c^{CW} - T)}{T}} \,\,\,\, {\textrm {and}} \,\,\,\, -\phi _{b}, \end{equation}

so the $\uparrow$ and $\downarrow$ phases coexist along the first-order phase boundary $h=0,\, 0\leqslant T \lt T_c^{\textit{CW}}$ , whose counterpart in the $T-\phi$ phase is the coexistence curve shown in figure 4(b). If $T \to T_c^{CW}$ from below, with $h=0$ , these minima merge at $T = T_c^{CW}$ , the coefficient of the quadratic term vanishes and elementary steps yield the mean-field order-parameter exponent $\beta _{MF}=1/2$ and interfacial-free-energy exponent $\mu _{MF} = 2\nu _{MF} = 1$ . The interfacial free energy vanishes for $T \gt T_c$ , for there are no coexisting phases and, therefore, no interface; and $g(\phi )$ has only one quadratic minimum if $T \gt T_c^{CW}$ and $h=0$ .

Figure 4.

Schematic plots of (a) the variational free energy $g(\phi )$ (see (2.11)) versus the order parameter $\phi$ , for magnetic field $h=0$ and temperature $T \lt T_c$ , with two equally deep minima at $\phi = \pm \phi _b$ , the LG equilibrium values; (b) the LG phase diagram in the $T-\phi$ plane showing the coexistence curve (also called the binodal), below which the phases coexist, and the spinodal curve, which is the locus of points at which $d^2g(\phi )/{\textrm{d}}\phi ^2 = 0$ .

To study non-uniform order-parameter configurations (see, e.g., Pandit & Wortis Reference Pandit and Wortis1982; Bray Reference Bray2002; Puri Reference Puri2004; Puri & Wadhawan Reference Puri and Wadhawan2009), we must minimise $\mathcal {F}_{CW}$ or $\mathcal {F}_{LG}$ . We illustrate this for $\mathcal {F}_{LG}$ below by considering BCs (cf. figure 3) that yield an interface for $T \lt T_c^{CW}$ and $h=0$ :

(2.13) \begin{eqnarray} \delta \mathcal {F}_{LG}/\delta \phi & =& 0 ; \nonumber \\ &\Rightarrow & \left [ -\sigma \nabla ^2\phi + k_B (T - T_c^{CW}) \phi + \frac {k_BT}{3} \phi ^3 \right ] = h(\mathbf {r}). \end{eqnarray}

If we assume that $\phi$ varies only in the $x$ direction and that $h(\mathbf {r}) = 0$ , we obtain the following ordinary differential equation, which must be solved with the BCs given in the second row:

(2.14) \begin{eqnarray} \sigma d^2\phi /{\textrm{d}}x^2 &=& k_B (T - T_c^{CW}) \phi + \frac {k_BT}{3} \phi ^3 ; \nonumber \\ {\textrm {BCs}}\,\,\,\, \phi (x) &=& \pm \phi _b \,\,\,\, {\textrm {at}} \,\,\,\, x = \pm \infty . \end{eqnarray}

This has the following well-known solution (see, e.g., Pandit & Wortis Reference Pandit and Wortis1982; Puri Reference Puri2004; Puri & Wadhawan Reference Puri and Wadhawan2009):

(2.15) \begin{eqnarray} \phi (x) &=& \phi _b \tanh \left [ (x-x_0)/\xi _{MF} \right ]; \nonumber \\ \xi _{MF} &=& \sqrt {\frac {\sigma }{2k_B(T_c^{CW} - T)}}; \end{eqnarray}

here, $x_0$ is the point at which this interfacial (or kink) profile goes through zero and $\xi$ is the width of the interface. If we portray the solutions of (2.14) in the phase space $[{\textrm{d}}\phi /{\textrm{d}}x,\,\phi ]$ , the bulk phases correspond to hyperbolic fixed points at ${\textrm{d}}\phi /{\textrm{d}}x=0,\,\phi =\pm \phi _{b}$ and the interfacial profile is a heteroclinic trajectory that goes from one of these fixed points to the other (see Pandit & Wortis Reference Pandit and Wortis1982). The mean-field free energy $f_{I,MF}$ (surface tension) for the kink solution follows by subtracting $\mathcal {F}_{LG}$ without an interface from its value with an interface,

(2.16) \begin{eqnarray} L^{(d-1)}f_{I,MF}(\phi, \nabla \phi ) &=& \int _{-\infty }^{\infty } {\textrm{d}}x \left [ \Delta g(\phi ) + \frac {\sigma }{2}[{\textrm{d}}\phi /{\textrm{d}}x]^2 \right ] , \nonumber \\ \Delta g(\phi ) &\equiv & \left [\frac {k_B}{2} (T - T_c^{CW}) [\phi ^2 - \phi _b^2] + \frac {k_BT}{12} [\phi ^4 - \phi _b^4]\right ] . \end{eqnarray}

If we multiply the first equation in (2.13) by ${\textrm{d}}\phi /{\textrm{d}}x$ and integrate, we obtain (see, e.g., Puri & Wadhawan Reference Puri and Wadhawan2009) for details)

(2.17) \begin{eqnarray} \frac {k_B}{2}\phi ^2 + \frac {k_B}{12}\phi ^4 + \frac {\sigma }{2}[{\textrm{d}}\phi /{\textrm{d}}x]^2 &=& \frac {k_B}{2}(T_c^{CW} - T)\phi _b^2 + \frac {k_B}{12}\phi _b^4 ,\nonumber \\ {\textrm {whence}} \,\,\,\, f_{I,MF} = \int _{-\infty }^{\infty } {\textrm{d}}x [{\textrm{d}}\phi /{\textrm{d}}x]^2. \end{eqnarray}

Note that as $T$ approaches $T_c^{CW}$ from below, $\phi _b \to 0$ , so $f_{I,MF}$ vanishes as $ \sim (T_c^{CW} - T)^{\mu _{MF}}$ , with $\mu _{MF} = 2\nu _{MF} = 1$ , where $\nu _{MF}$ is the exponent that characterises the divergence of $\xi _{MF}$ (see (2.15)). To go beyond this mean-field treatment, we must include the thermal fluctuations that are neglected by mean-field theory and which modify the mean-field values of critical exponents for $d \leqslant 4$ (see, e.g., Plischke & Bergersen Reference Plischke and Bergersen1994; Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Amit & Martin-Mayor Reference Amit and Martin-Mayor2005; Kardar Reference Kardar2007; Goldenfeld Reference Goldenfeld2018; Ma Reference Ma2018; Zinn-Justin Reference Zinn-Justin2021).

In most of this paper, we concentrate on physical parameter ranges in which systems are far from critical points. Therefore, we continue to use the LG framework outlined above and PDEs such as (2.13) with no noise terms that account for thermal fluctuations. If turbulence is present in the flows we consider, then turbulent fluctuations overpower thermal fluctuations; the latter are normally not required in studies of turbulent flows. For studies that suggest the inclusion of thermal noise in dissipation-range fluid turbulence see Betchov (Reference Betchov1977) and Bandak et al. (Reference Bandak, Goldenfeld, Mailybaev and Eyink2022).

We can also use the above inhomogeneous mean-field theory for the following:

1. (i) to study the wetting transition on an attractive substrate, which we illustrate by the schematic diagram in figure 5 (see, e.g., Cahn (Reference Cahn1977), Pandit & Wortis (Reference Pandit and Wortis1982) and Pandit, Schick & Wortis (Reference Pandit, Schick and Wortis1982) for a detailed treatment of this problem);

2. (ii) to study curvature corrections to the surface tension (see Fisher & Wortis Reference Fisher and Wortis1984 for details) at a curved interface; the dominant contribution at a curved interface yields Laplace’s result (see, e.g., Fisher & Wortis Reference Fisher and Wortis1984; Rowlinson & Widom Reference Rowlinson and Widom2013) for the pressure difference $\Delta P$ between the interior (say $A$ ) and exterior (say $B$ ) phases that are separated by a curved interface with surface tension $\sigma _{AB}$ and radius of curvature $R_{AB}$ :

   (2.18) \begin{equation} \Delta P = \frac {2 \sigma _{AB}}{R_{AB}}. \end{equation}

Figure 5.

Schematic diagram of a droplet on a solid substrate for (a) partial wetting, (b) complete drying and (c) complete wetting. The white, green and blue regions indicate the three phases A, B and C (the last of which is a solid substrate); the surface tensions between the coexisting phases are $\sigma _{AB}$ , $\sigma _{AC}$ and $\sigma _{BC}$ ; the contact angle $\theta$ is given by the Young–Dupré equation $\sigma _{AC}=\sigma _{BC}+\sigma _{AB} \cos (\theta )$ ; in (a) $0 \lt \theta \lt \pi$ , (b) $\theta = \pi$ and (c) $\theta = 0$ .

2.2. Kinetic Ising models

We turn now to a brief description of some time-dependent phenomena in multiphase systems. For pedagogical reasons, we have used the Ising model to illustrate various aspects of the equilibrium statistical mechanics of systems with interfaces. We adopt a similar strategy here by beginning with the kinetic Ising model (see, e.g., Kawasaki Reference Kawasaki1970; Puri Reference Puri2004; Puri & Wadhawan Reference Puri and Wadhawan2009), employing a mean-field approximation, taking a continuum limit and obtaining, therefrom, TDGL descriptions for the spatiotemporal evolution of the order parameter in two cases, viz., when it is (i) not conserved and (ii) conserved. Case (ii) yields the CH equation. Our discussion follows that of Puri & Wadhawan (Reference Puri and Wadhawan2009), which the reader should consult for details.

The Ising-model Hamiltonian (2.1) contains only one component (say the $z$ component in spin space), so it does not have intrinsic dynamics. We assume instead that the spins in this model are in contact with a heat bath (provided, e.g. by coupling to phonons in a crystalline solid) that leads to stochastic flipping of these spins. To study time-dependent phenomena we use, therefore, master equations for the conditional probability $P(\{S_i^0\},0;\{S_i\},t)$ that the spin at site $i$ is in the state $S_i$ at time $t$ , given that its $t=0$ state was $S_i^0$ ; henceforth, for notational convenience, we suppress the arguments $\{S_i^0\},0$ . (The braces indicate that we consider all the $N$ spins in the $d-$ dimensional hypercubic lattice). We give two such master equations:

1. (i) for the Glauber model (see Glauber (Reference Glauber1963), indicated by the subscript $G$ in (2.19)), in which the total magnetisation $\sum _i S_i$ is not conserved because we consider single-spin flips;

2. (ii) for the Kawasaki model (see Kawasaki Reference Kawasaki1966, Reference Kawasaki1970), indicated by the subscript $K$ in (2.28), in which the total magnetisation is conserved because we consider exchanges of nearest-neighbour spins.

2.2.1. Non-conserved order parameter: TDGL theory

The master equation for the Glauber model is

(2.19) \begin{align} \frac {{\textrm{d}}P_G(\{S_i\},t)}{dt} &= - \sum _{j=1}^N w_G(S_1, \ldots, S_j, \ldots, S_N|S_1, \ldots, -S_j, \ldots, S_N)P_G(\{S_i\},t)\nonumber \\ &\quad + \sum _{j=1}^N w_G(S_1, \ldots, -S_j, \ldots, S_N|S_1, \ldots, S_j, \ldots, S_N)P_G(\{S^{\prime}_i\},t) ; \end{align}

on the right-hand side the first and second terms arise, respectively, from the loss and gain of probability. The loss is associated with the spin flip $S_j \to -S_j$ ; the gain is because of the flip $S^{\prime}_j \to -S^{\prime}_j$ , in a state $\{S^{\prime}_i\}$ with $S^{\prime}_i = S_i\,\,\,\, \forall \,\, i\neq j$ and $S^{\prime}_j = -S_j$ . The master equation (2.19) has the long-time ( $t \to \infty$ ) equilibrium solution

(2.20) \begin{equation} P_{eq}(\{S_i\}) = \frac {1}{\mathcal {Z}}\exp \left [\frac {-\mathcal {H}}{k_BT}\right ] \end{equation}

if the transition matrix $w_G(\{S_i\}|\{S^{\prime}_i\})$ satisfies the condition of detailed balance (see, e.g., Van Kampen 1981; Martinelli 1999; Puri & Wadhawan Reference Puri and Wadhawan2009)

(2.21) \begin{equation} w_G(\{S_i\}|\{S^{\prime}_i\}) P_{eq}(\{S_i\}) = w_G(\{S^{\prime}_i\}|\{S_i\}) P_{eq}(\{S^{\prime}_i\}). \end{equation}

This condition does not specify $w_G(\{S_i\}|\{S^{\prime}_i\})$ completely; one convenient choice is the Suzuki–Kubo form (see, e.g., Suzuki & Kubo Reference Suzuki and Kubo1968; Puri & Wadhawan Reference Puri and Wadhawan2009)

(2.22) \begin{equation} w_G(\{S_i\}|\{S^{\prime}_i\}) = \frac {1}{2\tau _G}\left [1-\tanh \left (\frac {\Delta \mathcal {H}}{2k_BT}\right )\right ], \end{equation}

where $\tau _G$ is the relaxation time scale and $\Delta \mathcal {H}$ is the change in the energy because of the spin flip. If we substitute (2.22) in the master equation (2.19) and then calculate the mean value of the magnetisation at site $i$ ,

(2.23) \begin{equation} \langle S_i \rangle = \sum _{\{S_j\}} S_i P_G(\{S_j\},t), \end{equation}

standard calculations yield (see, e.g., Puri & Wadhawan Reference Puri and Wadhawan2009)

(2.24) \begin{equation} \tau _G\frac {{\textrm{d}}}{{\textrm{d}}t}\langle S_i \rangle = -\langle S_i \rangle +\left \langle \tanh \left [\frac {(h_i + J\sum _{j\in [nni]} S_j)}{(k_BT)}\right ]\right \rangle . \end{equation}

This equation can be solved analytically only if $d=1$ and $h_i=0\,\,\,\, \forall i$ (see, e.g., Glauber Reference Glauber1963; Siggia Reference Siggia1977; Pandit, Forgacs & Rujan Reference Pandit, Forgacs and Rujan1981). For $d \gt 1$ , we must either use approximations, such as mean-field theory, or finite-size calculations, numerical simulations or renormalisation groups (see, e.g., Pandit et al. Reference Pandit, Forgacs and Rujan1981; Münkel Reference Münkel1993; Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Kardar Reference Kardar2007; Täuber Reference Täuber2013; Goldenfeld Reference Goldenfeld2018; Ma Reference Ma2018). We follow Puri & Wadhawan (Reference Puri and Wadhawan2009) and use following the mean-field approach, which builds on the Curie–Weiss approximation (2.7) for the equilibrium magnetisation $M_i$ . A direct expansion of the second term on the right-hand side of (2.29) yields moments, of all orders, of the spins $S_i$ . In the mean-field dynamical model we use the approximation

(2.25) \begin{equation} \tau _G\frac {\textrm{d}}{{\textrm{d}}t}\langle S_i \rangle \simeq -\langle S_i \rangle + \tanh \left [\frac {(h_i + J\sum _{j\in [nni]} \langle S_j\rangle )}{(k_BT)}\right ] , \end{equation}

whose steady-state solution yields the Curie–Weiss self-consistency equation for the mean-field magnetisation given in (2.7). If we now take the continuum limit (2.9), we get

(2.26) \begin{equation} \tau _G\frac {\partial \phi }{\partial t} = -\frac {\delta \mathcal {F}_{CW}}{\delta \phi } , \end{equation}

and, if we make a small- $\phi$ expansion (valid, strictly speaking, near the Curie–Weiss critical point), we obtain the TDGL equation

(2.27) \begin{equation} \frac {\partial \phi }{\partial t} = -D\frac {\delta \mathcal {F}_{LG}}{\delta \phi } , \end{equation}

where we measure time in units of $\tau _G$ (so we set $\tau _G =1$ ). The TDGL equation (2.27) is often postulated phenomenologically and its right-hand side includes a Gaussian white-noise term $\theta (\mathbf {r},t)$ with zero-mean and a variance $2Dk_BT$ that satisfies the fluctuation-dissipation theorem; this is referred to as Model A in the critical-dynamics literature (see, e.g., Hohenberg & Halperin Reference Hohenberg and Halperin1977; Puri & Wadhawan Reference Puri and Wadhawan2009; Täuber Reference Täuber2013; Ma Reference Ma2018). This noise term plays an important role in the calculation of critical exponents, in general, and the dynamic critical exponent $z$ , in particular, which characterises the critical-point divergence of the correlation time $\tau _c \sim \xi ^z$ . At the simplest mean-field level, Model A yields $z^A_{MF} = 2$ .

2.2.2. Conserved order parameter: CH theory

The master equation for the Kawasaki model is (see Kawasaki Reference Kawasaki1966, Reference Kawasaki1970; Puri & Wadhawan Reference Puri and Wadhawan2009)

(2.28) \begin{align} & \frac {{\textrm{d}}P_K(\{S_i\},t)}{\textrm{d}t} \nonumber \\& \quad = - \sum _{j=1}^N\sum _{k\in [nnj]} w_K(S_1, \ldots, S_j,S_k, \ldots, S_N|S_1, \ldots, S_k,S_j, \ldots, S_N)P_K(\{S_i\},t)\nonumber \\& \qquad + \sum _{j=1}^N\sum _{k\in [nnj]} w_K(S_1, \ldots, S_k,S_j, \ldots, S_N|S_1, \ldots, S_j,S_k, \ldots, S_N)P_K(\{S^{\prime}_i\},t) , \end{align}

where the subscript $K$ stands for Kawasaki, we set $h=0$ in the Hamiltonian (2.1), and we work in the ensemble with fixed total magnetisation $M_{tot}=\sum _i S_i$ . Note that (2.28) accounts for nearest-neighbour spin exchanges $S_j \leftrightarrow S_k$ that conserve $M_{tot}$ . If we proceed as we did in the Glauber case, with the Suzuki–Kubo form (2.22) for the transition matrix $w_K$ , we obtain (see, e.g., Puri & Wadhawan (Reference Puri and Wadhawan2009) for details)

(2.29) \begin{align} \tau _K\frac {\textrm{d}}{{\textrm{d}}t}\langle S_i \rangle &= -2d\langle S_i \rangle + \sum _{j\in [nni]} \langle S_j\rangle \nonumber \\ &\quad+\sum _{k\in [nni]} \left \langle (1-S_iS_k)\tanh \Bigg [\frac {J\big (\sum _{j\in [nni;\neq k]}S_j-\sum _{j\in [nni;\neq i]}S_j\big )}{(k_BT)}\Bigg ]\right \rangle . \end{align}

The dynamical analogue of the Curie–Weiss approximation, followed by a continuum approximation, and an expansion in powers of $\phi$ , yields the CH equation (see Cahn & Hilliard Reference Cahn and Hilliard1958, Reference Cahn and Hilliard1959; Cahn Reference Cahn1961; Puri & Wadhawan Reference Puri and Wadhawan2009)

(2.30) \begin{equation} \frac {\partial \phi }{\partial t} = D \nabla ^2 \left [\frac {\delta \mathcal {F}_{LG}}{\delta \phi } \right ] , \end{equation}

where we set $\tau _K =1$ ; note that $\int {\textrm{d}}\mathbf {r} \phi$ is conserved by this CH equation (2.30), which is often postulated phenomenologically (see, e.g., Hohenberg & Halperin Reference Hohenberg and Halperin1977; Puri & Wadhawan Reference Puri and Wadhawan2009; Täuber Reference Täuber2013; Ma Reference Ma2018) and which has been used extensively in studies of phase separation, nucleation, and spinodal decomposition. If the right-hand side of the CH equation (2.30) includes a $\phi$ -conserving Gaussian white-noise term $\theta (\mathbf {r},t)$ with zero-mean and with a variance that satisfies the fluctuation-dissipation theorem, it is known as Model B or the Cahn–Hilliard–Cook equation. This noise term plays an important role in the calculation of critical exponents. For Model B, at the simplest level, (2.30) yields $z^B_{MF} = 4$ , which should be contrasted with the Model-A result $z^A_{MF} = 2$ .

Before we proceed to our discussion of the CHNS PDEs, we summarise below the types of time-dependent phenomena that can be studied by using the TDGL and CH equations.

1. (i) Space and time dependent correlation functions, e.g.

   (2.31) \begin{equation} C(\mathbf {r},t) \equiv \langle \phi (\mathbf {R},t_0) \phi (\mathbf {R}+\mathbf {r},t_0+t)\rangle _{t_0,\mathbf {R}} - \langle \phi (\mathbf {R},t_0)\rangle _{t_0,\mathbf {R}}\rangle ^2_{t_0,\mathbf {R}}, \end{equation}
   which can be used to study critical dynamics (see, e.g., Hohenberg & Halperin Reference Hohenberg and Halperin1977; Täuber Reference Täuber2013) and the kinetics of phase separation (see, e.g., Bray Reference Bray2002; Puri Reference Puri2004; Puri & Wadhawan Reference Puri and Wadhawan2009).

2. (ii) Consider the schematic plot of $g(\phi )$ , which appears in (2.11) for $\mathcal {F}_{LG}$ , in figure 4(b) (we set $h=0$ ) and the associated schematic phase diagram shown in figure 4(a). The early stages of phase separation can occur via nucleation, or droplet-type fluctuations, and spinodal decomposition, or small-amplitude long-wavelength fluctuations (see, e.g., Oki et al. (Reference Oki, Sagane and Eguchi1977), Gunton et al. (Reference Gunton, San Miguel and Sahni1983) and, in particular, figure 3 in Gunton et al. (Reference Gunton, San Miguel and Sahni1983), which has been reproduced from Oki et al. (Reference Oki, Sagane and Eguchi1977)). At the level of LG theory, the loci of points along which $g''(\phi ) = 0$ is the spinodal curve that provides a clear boundary between the regions with nucleation-dominated ( $g''(\phi ) \gt 0$ ) growth and spinodal decomposition ( $g''(\phi ) \lt 0$ ) in figure 4(a). If we go beyond mean-field or LG theory and consider the early stages of phase separation, the sharp spinodal curve is replaced by a crossover regime across which the phase-separating system moves from nucleation-dominated to spinodal-decomposition-initiated growth (see, e.g., Oki et al. Reference Oki, Sagane and Eguchi1977; Gunton et al. Reference Gunton, San Miguel and Sahni1983).

3. (iii) Late stages of phase separation. Here, we must distinguish Lifshitz–Slyozov and Lifshitz–Allen–Cahn scaling laws, for conserved and non-conserved $\phi$ , respectively; the former (latter) leads to domain growth in time $t$ that is characterised by the power-law growth of the length $\mathcal {L}(t) \sim t^{1/3}$ $\mathcal {L}(t) \sim t^{1/2}$ as discussed, e.g., in Lifshitz & Slyozov (Reference Lifshitz and Slyozov1961), Lifshitz (Reference Lifshitz1962), Fogedby & Mouritsen (Reference Fogedby and Mouritsen1988), Bray (Reference Bray2002), Puri (Reference Puri2004) and Puri & Wadhawan (Reference Puri and Wadhawan2009). Hydrodynamical effects can modify the domain-growth power-law exponent (see, e.g., Bray Reference Bray2002; Puri Reference Puri2004; Puri & Wadhawan Reference Puri and Wadhawan2009).

3.1. Binary-fluid CHNS

1. (i) For the binary-fluid case (see, e.g., Guo et al. Reference Guo, Yu, Lin, Wise and Lowengrub2021; Borcia et al. Reference Borcia, Borcia, Bestehorn, Sharma and Amiroudine2022) we use the following free-energy functional for our CHNS equations:

   (3.1) \begin{eqnarray} \mathcal F[\phi, \nabla \phi ] = \int _{\Omega } {\textrm{d}}\Omega \left [\mathcal V(\phi ) + \frac {3}{4} \sigma \epsilon |\nabla \phi |^2\right ], \end{eqnarray}
   where $\mathcal V(\phi )$ is the potential, $\sigma$ is the bare interfacial tension and $\epsilon$ is the width of the interface. The following two forms of $\mathcal V(\phi )$ are used in the CHNS literature:

   1. (a) Curie–Weiss-type potential,

      (3.2) \begin{eqnarray} \mathcal V_{CW}(\phi ) = \frac {a}{2} \left [(1+\phi )\ln (1+\phi ) + (1-\phi )\ln (1-\phi )\right ] - \frac {b}{2} \phi ^2; \end{eqnarray}
      this potential has the virtue that $\phi \in [-1,1]$ . (In some studies $\mathcal {V}_{CW}$ is called the Flory–Huggins logarithmic potential (see, e.g., Giorgini & Temam Reference Giorgini and Temam2020), as in the theory of polymer solutions (see, e.g., Rubinstein & Colby Reference Rubinstein and Colby2003).)

   2. (b) LG-type $\phi ^4$ potential,

      (3.3) \begin{eqnarray} \mathcal V_{LG}(\phi ) = \frac {3}{16} \frac {\sigma }{\epsilon }(\phi ^2-1)^2; \end{eqnarray}
      this parametrisation has the virtue that the global minima of $\mathcal V_{LG}(\phi )$ are at $\phi _b = \pm 1$ . However, now $\phi \in [-\infty, \infty ]$ . We use $\mathcal V_{LG}(\phi )$ in most of our CHNS studies.

2. (ii) If we allow the shear viscosity $\eta$ and the density $\rho$ to depend on $\phi$ , so that the two coexisting phases have different viscosities ( $\eta _1$ and $\eta _2$ ) and densities ( $\rho _1$ and $\rho _2$ ), far away from interfaces, the incompressible CHNS equations can be written as follows:

   (3.4) \begin{eqnarray} \rho (\phi ) (\partial _t {\boldsymbol {u}} + (\boldsymbol {u} \cdot \nabla ) \boldsymbol {u}) &=& -\nabla P + \nabla \cdot \left [\eta (\phi ) (\nabla \boldsymbol {u} + \nabla \boldsymbol {u}^{T})\right ] \nonumber \\ &-& \phi \nabla \mu + \rho (\phi ) \boldsymbol {g} - \alpha \boldsymbol {u} + \boldsymbol {f}^{ext}; \nonumber \\ \nabla \cdot \boldsymbol {u} &=& 0; \nonumber \\ \partial _t \phi + (\boldsymbol {u} \cdot \nabla ) \phi &=& M \nabla ^2\mu ; \nonumber \\ \mu &=& \frac {\delta \mathcal {F}}{\delta \phi }; \nonumber \\ \eta (\phi ) &=& \eta _1\left (\frac {1+\phi }{2}\right ) + \eta _2\left (\frac {1-\phi }{2}\right ); \nonumber \\ \rho (\phi ) &=& \rho _1\left (\frac {1+\phi }{2}\right ) + \rho _2\left (\frac {1-\phi }{2}\right ). \end{eqnarray}
   Here, $\boldsymbol {g}$ , $\boldsymbol {f}^{ext}$ and $\alpha$ are, respectively, the acceleration because of gravity, an external forcing, and the coefficient of friction (which is required typically in two-dimensional (2-D) settings, e.g. to account for bottom friction).

3. (iii) Boussinesq approximation. If the density difference between the two phases is not too large, in the NS part of the above equations we can use the Boussinesq approximation (see, e.g., Celani et al. Reference Celani, Mazzino, Muratore-Ginanneschi and Vozella2009; Boffetta et al. Reference Boffetta, Mazzino, Musacchio and Vozella2010; Lee & Kim Reference Lee and Kim2013a ,Reference Lee and Kim b ; Shah, Saeed & Yuan Reference Shah, Saeed and Yuan2017; Khan & Shah Reference Khan and Shah2019; Zanella et al. Reference Zanella, Tegze, Tellier and Henry2020; Forbes, Turner & Walters Reference Forbes, Turner and Walters2022; Huang & Yang Reference Huang and Yang2022). Furthermore, this approximation holds only when the vorticity is not too strong; otherwise, non-Boussinesq effects, such as interfacial instabilities, can emerge even at low Atwood numbers (see, e.g., Ramadugu, Perlekar & Govindarajan Reference Ramadugu, Perlekar and Govindarajan2022). In this Boussinesq approximation the NS equations can be written as follows:

   (3.5) \begin{eqnarray} \partial _t {\boldsymbol {u}} + (\boldsymbol {u} \cdot \nabla ) \boldsymbol {u} &=& - \frac {1}{\rho _0}\nabla P + \nu \nabla ^2 \boldsymbol {u}- \frac {1}{\rho _0}(\phi \nabla \mu ) + \frac {[\rho (\phi ) - {\rho }_0]}{{\rho }_0} \boldsymbol {g} - \alpha \boldsymbol {u},\nonumber \\ \nabla \cdot \boldsymbol {u} &=& 0, \end{eqnarray}
   where the mean density is $\rho _0 = (\rho _1+\rho _2)/2$ . This approximation neglects density differences except in the term with gravity. We can write
   (3.6) \begin{equation} \frac {(\rho (\phi ) - \rho _0)}{\rho _0} = -\frac {\Delta \rho }{\rho _0} \phi \equiv -\mathcal A \phi , \end{equation}
   with the Atwood number $\mathcal A \ll 1$ .

In table 2 we give the important dimensionless parameters that govern the states of the binary-fluid CHNS system. The larger the Reynolds number ${Re}$ , the more the inertia-induced turbulence, which can lead, e.g., to coarsening arrest (§ 6.3). The Bond number ${Bo}$ plays a crucial role in gravity-driven flows that include the evolution of antibubbles (§ 6.4) or a bubble passing through a liquid–liquid interface (§ 7.4); the Ohnesorge number ${Oh}$ is important in the coalescence of droplets or liquid lenses (§ 7.5); and the Weber ( ${We}$ ) and Capillary ( ${Ca}$ ) numbers affect interface-induced low- ${Re}\,$ turbulence (§ 6.5). The Cahn number governs the thickness of the interface. Other dimensionless numbers must be introduced as we increase the complexity of the multiphase system; e.g., in three-phase flows, density, surface-tension and viscosity ratios are relevant (§ 7.2), and the activity parameter is important in active CHNS systems (§§ 7.8 and 7.9).

Table 2.

Important dimensionless numbers for the binary-fluid CHNS system.

4.1. Pseudospectral method

The pseudospectral method, a widely used numerical technique for solving hydrodynamical PDEs, was pioneered over $50$ years ago by Patterson Jr & Orszag (Reference Patterson and Orszag1971) and has been used, inter alia, to study fluid turbulence (see, e.g., McWilliams Reference McWilliams1984; Canuto et al. Reference Canuto, Hussaini and Quarteroni1988; Pandit, Perlekar & Ray Reference Pandit, Perlekar and Ray2009; Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2010; San & Staples Reference San and Staples2012; Buaria & Sreenivasan Reference Buaria and Sreenivasan2022, Reference Buaria and Sreenivasan2023), magnetohydrodynamic (MHD) turbulence (see, e.g., Basu et al. Reference Basu, Sain, Dhar and Pandit1998; Müller & Biskamp Reference Müller and Biskamp2000; Verma Reference Verma2004; Gómez, Mininni & Dmitruk Reference Gómez, Mininni and Dmitruk2005; Sahoo et al. Reference Sahoo, Perlekar and Pandit2011, Reference Sahoo, Padhan, Basu and Pandit2020; Dritschel & Tobias Reference Dritschel and Tobias2012; Yadav, Miura & Pandit Reference Yadav, Miura and Pandit2022), CHNS turbulence (see, e.g., Scarbolo & Soldati Reference Scarbolo and Soldati2013; Perlekar et al. Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014; Pal et al. Reference Pal, Perlekar, Gupta and Pandit2016; Pandit et al. Reference Pandit2017; Vela-Martín & Avila Reference Vela-Martín and Avila2021; Pal et al. Reference Pal, Ramadugu, Perlekar and Pandit2022; Vela-Martín & Avila Reference Vela-Martín and Avila2022), the coalescence of droplets and liquid lenses (see, e.g., Padhan & Pandit Reference Padhan and Pandit2023b ); Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019), elastic turbulence in polymer solutions (see, e.g., Berti et al. Reference Berti, Bistagnino, Boffetta, Celani and Musacchio2008; Gupta & Pandit Reference Gupta and Pandit2017; Plan, Vincenzi & Gibbon Reference Plan, Gupta, Vincenzi and Gibbon2017; Singh et al. Reference Singh, Perlekar, Mitra and Rosti2024) and active turbulence and active droplets (see, e.g., Gao & Li Reference Gao and Li2017; Rana & Perlekar Reference Rana and Perlekar2020; Mukherjee et al. Reference Mukherjee, Singh, James and Ray2021, Reference Mukherjee, Singh, James and Ray2023; Li & Koch Reference Li and Koch2022; Gibbon et al. Reference Gibbon, Kiran, Padhan and Pandit2023; Kiran et al. Reference Kiran, Gupta, Verma and Pandit2023; Padhan & Pandit Reference Padhan and Pandit2023a ; Puggioni Leonardo & Stefano Reference Leonardo, Guido and Stefano2023; Backofen et al. Reference Backofen, Altawil, Salvalaglio and Voigt2024).

We employ the Fourier pseudospectral method, which is known for its accuracy and efficiency in comparison with other numerical methods (see, e.g., Orszag & Pao Reference Orszag and Pao1975; Canuto et al. Reference Canuto, Hussaini and Quarteroni1988); in its most common form, this method uses plane waves as basis functions. We solve the coupled CHNS equations in a periodic domain $\mathcal {D}\equiv l^d$ , with $l$ the length of the side of (typically) $d$ -dimensional hypercubic domain. To illustrate this method, we simplify (3.4) by making the following two key assumptions: all fluids have identical densities (i.e. $\rho _1 = \rho _2 = \rho _0 = 1$ ) and identical viscosities (i.e. $\eta _1 = \eta _2 = \eta$ ). The equations are given below for two dimensions and three dimensions.

4.1.1. Equations in two dimensions

In two dimensions, we solve the CHNS equations in its vorticity–stream function form (see, e.g., Boffetta & Ecke Reference Boffetta and Ecke2012; Pal et al. Reference Pal, Perlekar, Gupta and Pandit2016; Padhan & Pandit Reference Padhan and Pandit2023a ),

(4.1) \begin{eqnarray} \partial _t \omega + (\boldsymbol {u} \cdot \nabla ) \omega &=& \nu \nabla ^2 \omega -\alpha \omega - [\nabla \times (\phi \nabla \mu )]\cdot \hat {e_z} + f^\omega ;\nonumber \\ \nabla \cdot \boldsymbol {u} &=& 0; \quad \omega = (\nabla \times \boldsymbol {u})\cdot \hat e_z;\nonumber \\ \partial _t \phi + (\boldsymbol {u} \cdot \nabla ) \phi &=& M \nabla ^2 \mu ;\,\,\,\,\,\,\,\, \mu = -\frac {3}{2}\sigma \epsilon \nabla ^2 {\phi } + \frac {3}{4}\frac {\sigma }{\epsilon } ({\phi ^3} - {\phi }). \end{eqnarray}

4.1.2. Equations in three dimensions

In three dimensions, we use

(4.2) \begin{eqnarray} \partial _t {\boldsymbol {u}} + (\boldsymbol {u} \cdot \nabla ) \boldsymbol {u} &=& -\nabla P + \nu \nabla ^2\boldsymbol {u} - \phi \nabla \mu - \alpha \boldsymbol {u} + \boldsymbol {f}^{ext}; \nonumber \\ \nabla \cdot \boldsymbol {u} &=& 0; \nonumber \\ \partial _t \phi + (\boldsymbol {u} \cdot \nabla ) \phi &=& M \nabla ^2\mu ; \,\,\,\,\,\,\,\, \mu = -\frac {3}{2}\sigma \epsilon \nabla ^2 {\phi } + \frac {3}{4}\frac {\sigma }{\epsilon } ({\phi ^3} - {\phi }). \end{eqnarray}

Here $\nu = \eta / \rho _0$ is the kinematic viscosity. In most three-dimensional (3-D) applications, $\alpha = 0$ .

Our representation of the phase fields, velocity fields, vorticity fields and pressure fields utilises the (truncated) Fourier projection onto a grid of $N^d$ points, which are expressed as follows (carets denote spatial Fourier transform):

(4.3) \begin{eqnarray} \phi (\boldsymbol {x}, t) &=& \displaystyle \sum _{|k|\lt N} \hat {\phi }(\boldsymbol {k}, t) \exp (\dot {\iota } \boldsymbol {k} \cdot \boldsymbol {x}); \end{eqnarray}

(4.4) \begin{eqnarray} \boldsymbol {u}(\boldsymbol {x}, t) &=& \displaystyle \sum _{|k|\lt N} \hat {\boldsymbol {u}}(\boldsymbol {k}, t) \exp (\dot {\iota } \boldsymbol {k} \cdot \boldsymbol {x}); \end{eqnarray}

(4.5) \begin{eqnarray} \boldsymbol {\omega }(\boldsymbol {x}, t) &=& \displaystyle \sum _{|k|\lt N} \hat {\boldsymbol {\omega }}(\boldsymbol {k}, t) \exp (\dot {\iota } \boldsymbol {k} \cdot \boldsymbol {x}); \end{eqnarray}

(4.6) \begin{eqnarray} P(\boldsymbol {x}, t) &=& \displaystyle \sum _{|k|\lt N} \hat {P}(\boldsymbol {k}, t) \exp (\dot {\iota } \boldsymbol {k} \cdot \boldsymbol {x}). \end{eqnarray}

There are $N^d$ wavenumbers: $\boldsymbol {k} \equiv k_0 \boldsymbol {n} = k_0\sum _{i=1}^{d} n_i \hat {e}_i$ , where $k_0 = {2\pi }/{l}$ is the lowest wavenumber and $n_i$ are integers with values ranging from $(- ({N}/{2})+1)$ to ${N}/{2}$ . In each direction, the largest wavenumber is $k_{max} = {k_0 N}/{2}$ . (The spectral representation mentioned above aligns with representing the fields in physical space on an $N^d$ grid with uniform spacing $\Delta x = {l}/{N} = {\pi }/{k_{{max}}}$ .) We write, below, the CHNS equations in terms of the above Fourier modes.

4.1.3. Two-dimensional CHNS equations in Fourier space

We write (4.1) in Fourier space as follows:

(4.7) \begin{eqnarray} \partial _t \hat {\omega }(\boldsymbol {k}, t) &=& -\widehat {(\boldsymbol {u} \cdot \nabla \omega )} (\boldsymbol {k}, t) - \nu k^2 \hat {\omega }(\boldsymbol {k}, t) - \dot {\iota } \boldsymbol {k}\times \widehat {(\phi \nabla \mu )}(\boldsymbol {k}, t) -\alpha \hat {\omega }(\boldsymbol {k}, t); \end{eqnarray}

(4.8) \begin{eqnarray} \dot {\iota }\boldsymbol {k} \cdot \hat {\boldsymbol {u}}(\boldsymbol {k}, t) &=& 0; \,\,\,\,\,\, \hat {\omega }(\boldsymbol {k}, t) = \dot {\iota }\boldsymbol {k}\times \hat {\boldsymbol {u}}(\boldsymbol {k}, t); \end{eqnarray}

(4.9) \begin{eqnarray} \hat {\psi }(\boldsymbol {k}, t) &=& \frac {\hat {\omega }(\boldsymbol {k}, t)}{k^2}; \end{eqnarray}

(4.10) \begin{eqnarray} \partial _t \hat {\phi }(\boldsymbol {k}, t) &=& -\widehat {(\boldsymbol {u} \cdot \nabla \phi )}(\boldsymbol {k}, t) - M k^2 \hat {\mu }(\boldsymbol {k}, t); \end{eqnarray}

(4.11) \begin{eqnarray} \mu (\boldsymbol {k}, t) &=& \frac {3}{2}\sigma \epsilon k^2 \hat {\phi }(\boldsymbol {k}, t) + \frac {3}{4}\frac {\sigma }{\epsilon } [\hat {\phi ^3}(\boldsymbol {k}, t) - \hat {\phi }(\boldsymbol {k}, t)]. \end{eqnarray}

4.1.4. Three-dimensional CHNS equations in Fourier space

We write (4.2) in Fourier space as follows:

(4.12) \begin{align} \partial _t \hat {\boldsymbol {u}}(\boldsymbol {k}, t) = \mathbb {P}(\boldsymbol {k})\cdot \left [-\widehat {(\boldsymbol {u} \cdot \nabla \boldsymbol {u})} (\boldsymbol {k}, t) - \widehat {(\phi \nabla \mu )}(\boldsymbol {k}, t)\right ] - \nu k^2 \hat {\boldsymbol {u}}(\boldsymbol {k}, t) -\alpha \hat {\boldsymbol {u}}(\boldsymbol {k}, t); \end{align}

(4.13) \begin{eqnarray} \dot {\iota }\boldsymbol {k} \cdot \hat {\boldsymbol {u}}(\boldsymbol {k}, t) &=& 0; \end{eqnarray}

(4.14) \begin{eqnarray} \partial _t \hat {\phi }(\boldsymbol {k}, t) &=& -\widehat {(\boldsymbol {u} \cdot \nabla \phi )}(\boldsymbol {k}, t) - M k^2 \hat {\mu }(\boldsymbol {k}, t); \end{eqnarray}

(4.15) \begin{eqnarray} \mu (\boldsymbol {k}, t) &=& \frac {3}{2}\sigma \epsilon k^2 \hat {\phi }(\boldsymbol {k}, t) + \frac {3}{4}\frac {\sigma }{\epsilon } [\hat {\phi ^3}(\boldsymbol {k}, t) - \hat {\phi }(\boldsymbol {k}, t)]; \end{eqnarray}

where $\mathbb {P}$ is the transverse projection operator with components

(4.16) \begin{eqnarray} \mathbb {P}_{ij}(\boldsymbol {k}) = \left (\delta _{ij} - \frac {k_i k_j}{k^2}\right ). \end{eqnarray}

To perform Fourier transforms, we utilise the FFTW library (Frigo & Johnson Reference Frigo and Johnson1998), which employs the Cooley–Tukey fast-Fourier-transform algorithm (see Cooley & Tukey Reference Cooley and Tukey1965). The nonlinear terms such as $\boldsymbol {u} \cdot \nabla \boldsymbol {u},$ $\boldsymbol {u} \cdot \nabla \phi \,\,$ and $\mu \nabla \phi$ lead to convolution sums that require $N^{2d}$ operations; and the cubic nonlinear term $\phi ^3$ requires $N^{3d}$ operations in Fourier space. To bypass the high computational cost of these operations, we use the pseudospectral method in which we evaluate these nonlinear terms as follows. First, we transform the fields back to real space; and then we perform the multiplications in real space and transform them back to Fourier space again. This procedure improves greatly the computational efficiency of the pseudospectral method compared with the ordinary spectral methods. However, the Fourier pseudospectral method leads to aliasing errors while evaluating the nonlinear terms, which we remove by using the $1/2$ -dealiasing scheme (see, e.g., Patterson Jr & Orszag Reference Orszag1971; Hou & Li Reference Hou and Li2007; Padhan & Pandit Reference Padhan and Pandit2023a ), because of the cubic nonlinearity.

The computational cost of such a pseudospectral DNS (see, e.g., Moin & Mahesh Reference Moin and Mahesh1998) increases with the number of grid points, which must increase in turn to resolve small-scale structures. For statistically homogeneous and isotropic turbulence in three dimensions, the standard estimate (based on Kolmogorov 1941 phenomenology) is $N \sim Re^{9/4}$ , where $Re$ is the Reynolds number. A DNS of turbulence in the CHNS PDEs the grid spacing $\Delta x$ must be small enough to capture both the (Kolmogorv) dissipation scale and the widths of interfaces.

For the time integration of (4.7)–(4.15) we use a semi-implicit exponential time-differentiating Runge–Kutta2 (ETDRK2) approach, which treats the linear terms implicitly with their exact solutions as described in Cox & Matthews Reference Cox and Matthews2002). This method, which combines the ETD method with the RK2 method, offers several advantages over other numerical techniques, including enhanced accuracy, efficiency, stability and ease of implementation. We write the CHNS equations (4.7), (4.10) and (4.11) in the following general form:

(4.17) \begin{eqnarray} \frac {{\textrm{d}}q_1(t)}{{\textrm{d}}t} &=& \lambda _1 q_1(t) + \mathcal G(q_1, q_2); \end{eqnarray}

(4.18) \begin{eqnarray} \frac {{\textrm{d}}q_2(t)}{{\textrm{d}}t} &=& \lambda _2 q_2(t) + \mathcal H(q_1, q_2); \end{eqnarray}

here,

(4.19) \begin{eqnarray} q_1(t) &\equiv & \hat {\omega }(\boldsymbol {k}, t);\,\,\,\,\,\,\,\, \lambda _1 \equiv (-\nu k^2 - \alpha ); \nonumber \\ q_2(t) &\equiv & \hat {\phi }(\boldsymbol {k}, t); \,\,\,\,\,\,\,\, \lambda _2 \equiv \left(-\frac {3}{2} M \sigma \epsilon k^4 + \frac {3}{4} \frac {\sigma }{\epsilon } M k^2\right); \nonumber \\ \mathcal G(q_1, q_2) &\equiv & -\widehat {(\boldsymbol {u} \cdot \nabla \boldsymbol {\omega })} (\boldsymbol {k}, t) + \dot {\iota } \boldsymbol {k}\times \widehat {(\mu \nabla \phi )}(\boldsymbol {k}, t); \nonumber \\ \mathcal {H}(q_1, q_2) &\equiv & -\widehat {(\boldsymbol {u} \cdot \nabla \phi )}(\boldsymbol {k}, t) - \frac {3}{4} \frac {\sigma }{\epsilon }M k^2 \hat {\phi }(\boldsymbol {k}, t). \end{eqnarray}

In the ETDRK2 algorithm, (4.17) and (4.18) have the solution

(4.20) \begin{eqnarray} a_1^n &=& q_1^n \exp (\lambda _1 \Delta t) + \mathcal {G}(q_1^n, q_2^n) \left (\frac {\exp (\lambda _1 \Delta t) - 1}{\lambda _1}\right ); \end{eqnarray}

(4.21) \begin{eqnarray} a_2^n &=& q_2^n \exp (\lambda _2 \Delta t) + \mathcal {H}(q_1^n, q_2^n) \left (\frac {\exp (\lambda _2 \Delta t) - 1}{\lambda _2}\right ); \end{eqnarray}

(4.22) \begin{eqnarray} q_1^{n+1} &=& a_1^{n} + \left [\mathcal {G}(a_1^n, a_2^n) - \mathcal {G}(q_1^n, q_2^n)\right ] \left (\frac {\exp (\lambda _1 \Delta t) - 1 - \lambda _1 \Delta t}{\lambda _1^2 \Delta t}\right ); \end{eqnarray}

(4.23) \begin{eqnarray} q_2^{n+1} &=& a_2^{n} + \left [\mathcal {H}(a_1^n, a_2^n) - \mathcal {H}(q_1^n, q_2^n)\right ] \left (\frac {\exp (\lambda _2 \Delta t) - 1 - \lambda _2 \Delta t}{\lambda _2^2 \Delta t}\right ); \end{eqnarray}

for simplicity, we use $q_1^{n+1} \equiv q_1(t + \Delta t)$ and $q_1^n \equiv q_1(t)$ .

4.2. Volume-penalised CHNS

Fluid–structure interactions play a crucial role in multiphase flows, where solid boundaries or immersed objects influence phase separation, interfacial dynamics and transport properties (see, e.g., HUAT Reference Huat2015; Zheng et al. Reference Zheng, Wang, Triantafyllou and Karniadakis2021; Hester et al. Reference Hester, Carney, Shah, Arnheim, Patel, Di Carlo and Bertozzi2023; Pavuluri et al. Reference Pavuluri, Holtzman, Kazeem, Mohammed, Seers and Rabbani2023; Treeratanaphitak & Abukhdeir Reference Treeratanaphitak and Abukhdeir2023; Hou et al. Reference Hou, Wei, Iqbal, Yang, Wang, Ren and Yan2024). These interactions are challenging in computational fluid dynamics because of the enforcement of BCs at the fluid–solid interface (see, e.g., Mokbel, Abels & Aland Reference Mokbel, Abels and Aland2018; Pramanik, Verstappen & Onck Reference Pramanik, Verstappen and Onck2024). The numerical implementation of these BCs is especially difficult when the immersed solid is in motion. The immersed boundary method has been used widely in computational fluid dynamics to handle such interactions (see, e.g., Mittal & Iaccarino Reference Mittal and Iaccarino2005; Mittal & Seo Reference Mittal and Seo2023). The volume penalisation method (VPM) is a type of immersed boundary method that is simple and versatile in modelling objects in fluid flows; it is gaining in popularity, given its ease of implementation (see, e.g., Kolomenskiy & Schneider Reference Kolomenskiy and Schneider2009; Kadoch et al. Reference Kadoch, Kolomenskiy, Angot and Schneider2012; Morales et al. Reference Morales, Leroy, Bos and Schneider2014; Engels et al. Reference Engels, Kolomenskiy, Schneider and Sesterhenn2016, Reference Engels, Truong, Farge, Kolomenskiy and Schneider2022; Hester, Vasil & Burns Reference Hester, Vasil and Burns2021; Sinhababu & Bhattacharya Reference Sinhababu and Bhattacharya2022; Mittal & Seo Reference Mittal and Seo2023; Puggioni Leonardo & Stefano Reference Leonardo, Guido and Stefano2023). Without imposing explicit BCs, the VPM uses an extra force as a penalisation term in the classical NS equations; the modified equation is known as the Brinkman model (see, e.g., Angot, Bruneau & Fabrie Reference Angot, Bruneau and Fabrie1999). The VPM approximates solid boundaries in a fluid by applying rapid linear damping within a fictitious solid region. This approach allows for the simulation of complex, moving objects within general numerical solvers without the need for specialised algorithms or boundary-conforming grids. The VPM treats solids as porous media that are characterised by negligible permeability, resulting in the velocity of the adjacent fluid becoming zero at the fluid–solid interface. This method can be easily incorporated into a DNS solver that employs periodic BCs. The VPM has been successfully implemented in the NS equations, passive scalar equations, MHD equations, CH equations and the Toner–Tu–Swift–Hohenberg equations (see, e.g., Kadoch et al. Reference Kadoch, Kolomenskiy, Angot and Schneider2012; Morales et al. Reference Morales, Leroy, Bos and Schneider2014; Engels et al. Reference Engels, Kolomenskiy, Schneider and Sesterhenn2016, Reference Engels, Truong, Farge, Kolomenskiy and Schneider2022; Sinhababu, Bhattacharya & Ayyalasomayajula Reference Sinhababu, Bhattacharya and Ayyalasomayajula2021; Puggioni Leonardo & Stefano Reference Leonardo, Guido and Stefano2023). We extend the VPM method to the CHNS equations as follows (for illustration we use the 2-D CHNS system):

(4.24) \begin{eqnarray} \partial _t \omega + (\boldsymbol {u} \cdot \nabla ) \omega &=& \nu \nabla ^2 \omega -\alpha \omega - [\nabla \times (\phi \nabla \mu )]\cdot \hat {e_z} + f^\omega - \nabla \times \left [\frac {\chi }{\eta _p} \boldsymbol {u}\right ];\nonumber \\ \nabla \cdot \boldsymbol {u} &=& 0; \quad \omega = (\nabla \times \boldsymbol {u})\cdot \hat e_z;\nonumber \\ \partial _t \phi + [(1-\chi )\boldsymbol {u}] \cdot \nabla \phi &=& \nabla \cdot ([M(1 - \chi ) + M_p \chi ] \nabla \mu ). \end{eqnarray}

Here, $\eta _p$ and $M_p$ are the penalisation parameters (or permeabilities) associated with the velocity and $\phi$ fields, respectively. Here $\chi$ is the mask function to distinguish solid and fluid regions, which is defined as follows:

(4.25) \begin{eqnarray} \chi (\boldsymbol {x}) &=& \begin{cases} 1 &\text {for}\ \boldsymbol {x} \in \Omega _s ;\\0 &\text {for}\ \boldsymbol {x} \in \Omega _f ; \end{cases} \end{eqnarray}

where $\Omega _s$ and $\Omega _f$ are the volumes representing solid and fluid regions. The term $[M(1 - \chi ) + M_p \chi ]$ takes into account the no-flux BCs at the solid–fluid interface $\nabla \phi \cdot \hat {n}|_{\Omega _f} = 0$ ; $\hat {n}$ is the unit vector that is normal to the solid wall.

6.1. Rayleigh–Taylor instability: CHNS (two dimensions)

We illustrate how the CHNS framework can be used to simulate the gravity-induced RT instability (see, e.g., Drazin Reference Drazin2002; Charru Reference Charru2011; Livescu, Wei & Petersen Reference Livescu, Wei and Petersen2011), which occurs when a dense fluid is positioned, initially, above a less dense one (see, e.g., Tryggvason Reference Tryggvason1988; Youngs Reference Youngs1992; Young et al. Reference Young, Tufo, Dubey and Rosner2001; Kadau et al. Reference Kadau, Germann, Hadjiconstantinou, Lomdahl, Dimonte, Holian and Alder2004; Celani et al. Reference Celani, Mazzino, Muratore-Ginanneschi and Vozella2009; Livescu et al. Reference Livescu, Wei and Petersen2011; Lee & Kim Reference Lee and Kim2013; Khomenko et al. Reference Khomenko, Díaz, De, Angel and Luna2014; Gibbon et al. Reference Gibbon, Gupta, Krstulovic, Pandit, Politano, Ponty, Pouquet, Sahoo and Stawarz2016; Talat et al. Reference Talat, Mavrič, Hatić, Bajt and Šarler2018; Gonzalez-Gutiérrez & de Andrea González Reference Gonzalez-Gutiérrez and de Andrea González2019; Khan & Shah Reference Khan and Shah2019; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2021; Garoosi & Mahdi Reference Garoosi and Mahdi2022; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022; Pandya & Shkoller Reference Pandya and Shkoller2023). The CHNS framework has been shown to be useful in studying RT instabilities in immiscible fluids (see, e.g., Celani et al. Reference Celani, Mazzino, Muratore-Ginanneschi and Vozella2009; Lee & Kim Reference Lee and Kim2013; Gibbon et al. Reference Gibbon, Pal, Gupta and Pandit2016b ). We demonstrate this by an illustrative pseudospectral DNS of the 2-D CHNS equations, with the Boussinesq approximation (3.5)–(3.6), in a 2-D rectangular box $(L_x, L_y) = (2\pi, 4\pi )$ ; we incorporate impenetrable boundaries in the $y$ -direction by using the volume-penalisation method, withsix grid points on both the top and bottom boundaries for penalisation. The initial conditions we use are as follows:

(6.1) \begin{eqnarray} \phi (x, y, t = 0) &=& \tanh \left [\frac {y - L_y/2 - h_0 \cos (m x)}{\epsilon /2}\right ];\\ \omega (x, y, t = 0) &=& 0;\nonumber \end{eqnarray}

where $h_0$ is the amplitude of the perturbation we impose on top of a flat interface at $L_y/2$ . In figure 6(a,b), we present pseudocolour plots of $\phi$ and $\omega$ ; these show how the RT instability develops in time (which increases from left to right); the top phase (in red) has a higher density than the bottom phase (in blue); and the Atwood number is $\mathcal A = 0.6$ . We can quantify the temporal growth of the RT instability by computing the normalised square of the vertical velocity

(6.2) \begin{eqnarray} \frac {\langle u_y^2(t)\rangle }{U_0^2} \equiv \frac {2}{L_xL_y}\int _{0}^{L_x} {\textrm{d}}x \int _{\frac {L_y}{4}}^{\frac {3L_y}{4}} u_y^2(t) {\textrm{d}}y, \end{eqnarray}

with $U_0 = g h_0^2/\nu$ the natural velocity scale for this problem, which is shown in the red semi-log plot in figure 6(c). The blue semi-log plot in figure 6(c) shows the maximum scaled height $h(t)/h_0$ of the interface (the maximal deviation of the $\phi = 0$ contour line from the $t=0$ interface position $L_y/2$ ). The linear regions in both semi-log plots in figure 6(c) are consistent with the exponential growth of the RT instability at early times.

Figure 6.

Rayleigh–Taylor instability in the 2-D CHNS system: pseudocolour plots of (a) the $\phi$ field and (b) the corresponding vorticity field at different times (increasing from left to right). (c) Semi-log plots versus the scaled time $T/T_0$ of $[\langle u_y^2(t)\rangle ]/[U_0^2]$ (red) and the maximum scaled height $h(t)/h_0$ (blue) (see the text). The simulation box size is $(L_x, L_y) = (2\pi, 4\pi )$ with $512\times 1024$ grid points; $\nu = 0.01,\,\, \alpha = 0,\,\, \sigma = 0.1,\,\, g = 1, \mathcal A = 0.6,\,\, \rho _0 = 1,\,\, h_0 = 0.1$ and $m = 2$ . The characteristic velocity and time scales are $U_0 = g h_0^2/\nu$ and $T_0 \equiv h_0/U_0 = {\nu }/{g h_0}$ . The simulation box is periodic in the $x$ -direction; BCs in the $y$ -direction are implemented by the volume-penalisation method (see the text).

For a CHNS-based study of the RT instability in three dimensions, we refer the reader to Gibbon et al. (Reference Gibbon, Pal, Gupta and Pandit2016b ). This DNS study was designed to obtain $E_{\infty }(t)$ (5.3) by obtaining large- $m$ estimates for the $L^m$ norm of the energy $E(t)$ (5.2); and the results, based on the BKM-type criterion (5.4), were consistent with no finite-time singularity, given the resolution of the DNS runs.

6.2. Kelvin–Helmholtz instability: CHNS (two dimensions)

Figure 7.

Our DNS of the KH instability in the 2-D CHNS system: pseudocolour plots of (a) the $\phi$ field and (b) the corresponding vorticity field at different simulation times (increasing from left to right); the vorticity field is normalised by its absolute maximum value. (c) Plots versus the scaled time $t/T_0$ of $[\langle u_y^2(t) \rangle ]/[U_0^2]$ (red semi-log plot), where $T_0 \equiv h_0/U_0$ , and $\theta (t)$ (blue linear plot) (see (6.4) and (6.5)). The simulation box size is $(L_x, L_y) = (2\pi, 2\pi )$ with $1024 \times 1024$ grid points; and $U_0 = 2, h_0 = 0.1, \nu = 0.01; \alpha = 0.01, \sigma = 0.05, g = 1, \mathcal A = 0.01, \rho _0 = 1$ . The simulation box is periodic in the $x$ -direction and we use volume penalization in the $y$ -direction to incorporate solid boundaries; we incorporate impenetrable boundaries in the $y$ -direction by using the volume-penalisation method, with $6$ grid points on both the top and bottom boundaries for penalisation.

If there is a significant difference in the velocities of two fluid layers, which are separated by an interface, then this interface becomes unstable because of the well-known KH instability (see, e.g., Aref & Siggia Reference Aref and Siggia1981; Drazin Reference Drazin2002; Cushman-Roisin Reference Cushman-Roisin2005; Charru Reference Charru2011; Yilmaz et al. Reference Yilmaz, Davidson, Edis and Saygin2011; Lee & Kim Reference Lee and Kim2015; Tian & Chen Reference Tian and Chen2016; Hoshoudy & Cavus Reference Hoshoudy and Cavus2018; Budiana et al. Reference Budiana, Pranowo, Indarto and Deendarlianto2020; Zhou, Dong & Li Reference Zhou, Dong and Li2020; Delamere et al. Reference Delamere, Barnes, Ma and Johnson2021; Jia et al. Reference Jia, Xie, Deng, He, Yang and Fu2023; Kumar et al. Reference Kumar, Bandyopadhyay, Singh, Dharodi and Sen2023) that plays an important role in various marine, geophysical, solar and astrophysical processes (see, e.g., Mishin & Tomozov Reference Mishin and Tomozov2016; Gregg et al. Reference Gregg, D’Asaro, Riley and Kunze2018). The CHNS framework has been used to study KH instabilities in binary and ternary fluids (see Lee & Kim Reference Lee and Kim2015). To illustrate how this is done, we perform a pseudospectral DNS of the CHNS equations within the Boussinesq approximation (3.5)–(3.6) in a 2-D box $(L_x, L_y) = (2\pi, 2\pi )$ . Our simulation box is periodic in the $x$ -direction and we incorporate solid boundaries in the $y$ -direction by using the volume-penalisation method, where we consider six grid points on both the top and bottom sides for penalisation. We use the following stably stratified initial conditions, so that there is no RT instability:

(6.3) \begin{align} \phi (x, y, 0) &= \tanh \left [\frac {y - L_y/2 - h_0 \sin (2x)}{\epsilon /2}\right ];\nonumber \\ u_x(x, y, 0) &= U_0 \phi (x, y, 0); \nonumber\\ u_y(x, y, 0) &= 0. \end{align}

In figures 7(a) and 7(b), we portray the $\phi$ and vorticity fields, respectively, via pseudocolour plots. We also quantify the temporal growth of the KH instability by computing the normalised square of the vertical velocity

(6.4) \begin{eqnarray} \frac {\langle u_y^2(t) \rangle }{U_0^2} \equiv \frac {2}{L_xL_y} \int _{0}^{L_x}{\textrm{d}}x\int _{\frac {L_y}{4}}^{\frac {3L_y}{4}} u_y^2(t) {\textrm{d}}y, \end{eqnarray}

which we plot versus the scaled time $t/T_0$ in figure 7(c) (red semi-log plot), where $T_0 \equiv h_0/U_0$ ; the limits on the integral over $y$ are chosen to exclude the effects of the boundaries. The initial increase in ${\langle u_y^2(t) \rangle }/{U_0^2}$ signals the KH instability; this ratio decreases eventually because, in our DNS, the shear is present only in the initial condition. We also calculate the momentum thickness (Aref & Siggia Reference Aref and Siggia1981)

(6.5) \begin{eqnarray} \theta (t) = \int _{\frac {L_y}{4}}^{\frac {3L_y}{4}} {\textrm{d}}y \frac {\sqrt {\langle u_x^2(y, t)\rangle _x}}{\sqrt {\langle u_x^2(L_y/2, t)\rangle _x}}; \end{eqnarray}

a linear plot of $\theta (t)$ versus the scaled time $t/T_0$ is shown in figure 7(c) (blue line).

6.3. Phase separation in the binary-fluid CHNS

A homogeneous binary-fluid mixture spontaneously phase separates into domains of pure phases when the system is quenched from a high temperature to a low temperature, which lies below the critical temperature $T_c$ (see, e.g., Chaikin et al. Reference Chaikin, Lubensky and Witten1995; Bray Reference Bray2002; Puri & Wadhawan Reference Puri and Wadhawan2009); in equilibrium, pure phases are separated by an interface. If the order parameter is conserved, the early stages of phase separation proceed via nucleation or spinodal decomposition; thereafter, the system approaches the state with complete phase separation via the coarsening of domains of the pure phases (see § 2.2.2).

The CHNS framework has been used for studying the coarsening of phase-separating binary-fluid mixtures in both two and three dimensions (see, e.g., Berti et al. Reference Berti, Boffetta, Cencini and Vulpiani2005; Perlekar, Pal & Pandit Reference Perlekar, Pal and Pandit2017; Wang et al. Reference Wang, Altschuh, Ratke, Zhang, Selzer and Nestler2019). To study the coarsening process in a symmetric binary fluid mixture, we define the domain length scale ${L}(t)$ and integral length scale ${L}_I(t)$ in terms of the phase-field spectrum $S(k, t)$ and energy spectrum $E(k, t)$ , respectively, as follows:

(6.6) \begin{eqnarray} {L}(t) &=& 2\pi \frac {\displaystyle \sum _{k} k^{-1}S(k, t)}{\displaystyle \sum _{k} S(k, t)};\,\,\,\,\,\,\,\, {L}_I(t) = 2\pi \frac {\displaystyle \sum _{k} k^{-1}E(k, t )}{\displaystyle \sum _{k} E(k, t)}; \nonumber \\ S(k, t) &=& \displaystyle \sum _{k-1/2\lt k'\lt k+1/2} \left [\hat {\phi }(\mathbf {k}',t)\cdot \hat {\phi }(-\mathbf {k}',t)\right ]; \nonumber \\ E(k,t) &=& \frac {1}{2}\displaystyle \sum _{k-1/2\lt k'\lt k+1/2} \left [\hat { \boldsymbol {u}}(\mathbf {k}',t)\cdot \hat { \boldsymbol {u}}(-\mathbf {k}',t)\right ]. \end{eqnarray}

The length scale $L(t)$ displays self-similar power-law growth with $L(t) \sim t^{\beta }$ : (i) in the diffusion-dominated regime (governed by the CH equation) $\beta \simeq 1/3$ , the well-known Lifshitz–Slyozov exponent (see, e.g., Lifshitz & Slyozov Reference Lifshitz and Slyozov1961; Bray Reference Bray2002; Puri & Wadhawan Reference Puri and Wadhawan2009); (ii) in the viscous-hydrodynamic regime, $\beta \simeq 1$ (see Siggia Reference Siggia1979; Perlekar et al. Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014); and (iii) in the inertia-dominated regime, $\beta \simeq 2/3$ (Furukawa Reference Furukawa2000).

Figure 8.

Coarsening in a binary-fluid mixture from our DNS of the 2-D CHNS equations. (a) Pseudogreyscale plots of the phase-field $\phi (\boldsymbol {x}, t)$ at simulation times $t = 0, 0.24, 5.52 (\equiv t_{*}), 8.88$ , increasing from left to right; (b) pseudocolour plots of the corresponding vorticity fields $\omega (\boldsymbol {x}, t)$ (normalised by their absolute maximum for ease of visualisation). (c) Log–log plots versus $t/t_\nu$ , with $t_{\nu } = \nu ^3/\sigma ^2$ , of the scaled lengths $L(t)/L_0$ (red) and ${L}_I(t)/L_0$ (blue), with $L_0 = 2\pi$ the side of the simulation domain. (d) The energy spectrum (red line) $E(k,t=t_*)$ and the phase-field spectrum (blue line) $S(k, t=t_*)$ show power-law behaviour with scaling exponents $\sim k^{-3}$ and $\sim k^{-2}$ , respectively.

We carry out an illustrative pseudospectral DNS of the CHNS (two dimensions) equations to obtain the scaling exponent $\beta \simeq 2/3$ in the inertial-hydrodynamic regime. We illustrate domain coarsening by the pseudocolour plots of $\phi$ and $\omega$ in figures 8(a) and 8(b). The power-law growth of $L(t)$ and ${L}_I(t)$ , in the low-viscosity inertial regime, is shown in figure 8(c); the growth of both these lengths is consistent with the same power-law exponent $\beta \simeq 2/3$ . We show the energy spectrum and phase-field spectrum in figure 8(d); $E(k) \sim k^{-3}$ (consistent with the scaling exponent obtained in 2-D turbulence (Boffetta & Ecke Reference Boffetta and Ecke2012)) and $S(k) \sim k^{-2}$ (consistent with the results obtained by Furukawa (Reference Furukawa2000)). We solve the 2-D CHNS equations via a pseudospectral DNS in a doubly periodic box of size $2\pi \times 2\pi$ with grid points $1024 \times 1024$ . The relevant dimensionless numbers are the Reynolds number $Re = u_{rms} L_I/\nu$ and the Ohnesorge number $Oh = \nu ({\rho }/{\sigma L_I})^{1/2}$ . The simulation parameters are $\nu = 5\times 10^{-3}, \sigma = 0.4, \rho = 1, \epsilon = 0.018$ and $M = 10^{-4}$ . The dimensionless numbers, calculated at $t = t^*$ (see figure 8), are $Re = 24$ and $Oh = 0.01$ .

Fluid turbulence can mix immiscible fluids and lead to the arrest of phase separation, which is also referred to as coarsening arrest. The CHNS PDEs have been used to study such turbulence-induced coarsening arrest in both two dimensions and three dimensions (see, e.g., Perlekar et al. Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014, Reference Perlekar, Pal and Pandit2017). If parameters are chosen such that the CH equations lead to complete phase separation of the two fluids in the binary mixture, then the inclusion of turbulence, obtained by forcing the coupled CHNS equations, suppresses this phase separation. This can be characterised by using the lengths and spectra that we have defined in (6.6). In the absence of turbulence, the spectrum $S(k,t)$ displays an inverse cascade to small wavenumbers $k$ as time $t$ increases; this leads to the power-law growth of ${L}_I(t)$ as $t\to \infty$ , with the Lifshitz–Slyozov, viscous-hydrodynamic and inertial-hydrodynamic exponents mentioned above. This inverse cascade and the associated power-law growth of ${L}_I(t)$ are arrested by turbulence. In particular, it has been shown (see the pseudospectral DNS of Perlekar et al. (Reference Perlekar, Pal and Pandit2017) and the Lattice–Boltzmann study of Perlekar et al. (Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014) in two dimensions and three dimensions, respectively), which ${L}_I(t) \sim L_H$ as $t \to \infty$ , where the Hinze length scale $L_H \sim \varepsilon _{inj}^{-2/5}\sigma ^{3/5}$ , with $\varepsilon _{inj}$ the energy injection into the NS part of the CHNS system. We will show in § 7.8 that active turbulence, which occurs in the active-CHNS system, also leads to qualitatively similar coarsening arrest.

6.4. Spatiotemporal evolution of antibubbles

A shell of a low-density fluid inside a high-density fluid is known as an antibubble. It seems to have been described first by Hughes and Hughes (see, e.g., Hughes & Hughes Reference Hughes and Hughes1932; Dorbolo et al. Reference Dorbolo, Vandewalle, Reyssat and Quéré2006; Kim & Vogel Reference Kim and Vogel2006; Kalelkar Reference Kalelkar2017; Vitry et al. Reference Vitry, Dorbolo, Vermant and Scheid2019; Pal et al. Reference Pal, Ramadugu, Perlekar and Pandit2022; Zia et al. Reference Zia, Nazir, Poortinga, van and Cornelus2022). An antibubble has two surfaces, and a certain volume of fluid is trapped between these two surfaces. Clearly, an antibubble is unstable under gravity: if the fluid in its inner core is denser than that in the outer core, the antibubble rises; and the fluid in the shell forms a bulb at the top while its bottom thins until the shell collapses completely. For experimental investigations of the dynamics of an antibubble, see, for example, Dorbolo et al. (Reference Dorbolo, Vandewalle, Reyssat and Quéré2006), Kim & Vogel (Reference Kim and Vogel2006), Vitry et al. (Reference Vitry, Dorbolo, Vermant and Scheid2019) and Zia et al. (Reference Zia, Nazir, Poortinga, van and Cornelus2022); and for theoretical work consult, for example, Scheid et al. (Reference Scheid, Dorbolo, Arriaga and Rio2012), Zou et al. (Reference Zou, Ji, Yuan, Ruan and Fu2013), Sob’yanin (Reference Sob’yanin2015) and Pal et al. (Reference Pal, Ramadugu, Perlekar and Pandit2022). It is important to note that the inherent instability of antibubbles makes experimental studies challenging; in some cases surfactant molecules have to be introduced to obtain some stabilisation of an antibubble. Furthermore, antibubbles have several applications that include sonoporation (see, e.g., Kotopoulis et al. Reference Kotopoulis, Delalande, Popa, Mamaeva, Dimcevski, Gilja, Postema, Gjertsen and McCormack2014), drug delivery (see, e.g., Johansen et al. Reference Johansen, Kotopoulis, Poortinga and Postema2015 a); Kotopoulis et al. Reference Kotopoulis, Johansen, Gilja, Poortinga and Postema2015) and active leakage detection (see, e.g., Johansen et al. Reference Johansen, Kotopoulis, Poortinga and Postema2015 b).

It was recognised by Pal et al. (Reference Pal, Ramadugu, Perlekar and Pandit2022) that the CHNS system provides an ideal theoretical framework for the elucidation of the spatiotemporal evolution of antibubbles. The initial condition for $\phi$ is an annulus in two dimensions or a shell in three dimensions of the lighter fluid (shown in red) with the heavier fluid (shown in blue) both inside and outside the shell. The initial outer and inner radii of the shell are, respectively, $R_0$ and $R_1$ , and $h_0$ is the initial thickness of the antibubble shell, so it is natural to define the Bond number Bo with $L_0$ replaced by $h_0$ (see table 2). The spatiotemporal evolution of such an antibubble in the 2-D CHNS system is shown via pseudocolour plots of $\phi$ in figures 9(a) and 9(c) at representative times; the associated evolution of the vorticity $\omega$ is given, respectively, in figures 9(b) and 9(d). (These figures have been provided very kindly by N. Pal.) From figures 9(a) and 9(c) we see that gravity induces a thinning of the bottom of the antibubble while forming a slight dome on the top of it; eventually this makes the antibubble rupture. A similar rupture also occurs for antibubbles in the 3-D CHNS system (see Pal et al. Reference Pal, Ramadugu, Perlekar and Pandit2022 for details). These DNSs yield the scaled rupture time $\tau _1/\tau _g$ , with $\tau _g \equiv \sqrt {R_0/\mathcal {A}g}$ , the velocity $v_{rim}$ of the retracting rim of the collapsing antibubble, and they find that $v_{rim} \sim \sqrt {\sigma }$ , in agreement with the theoretical estimate of Sob’yanin (Reference Sob’yanin2015) and experiments of Scheid et al. (Reference Scheid, Dorbolo, Arriaga and Rio2012) and Zou et al. (Reference Zou, Ji, Yuan, Ruan and Fu2013). Furthermore, this CHNS study (Pal et al. Reference Pal, Ramadugu, Perlekar and Pandit2022) obtains the dependence of $\tau _1/\tau _g$ on the Bond number Bo and shows, by obtaining the energy ( $E(k)$ ) and concentration ( $S(k)$ ) spectra, that the rupture of the antibubble leads to some turbulence. Finally Pal et al. (Reference Pal, Ramadugu, Perlekar and Pandit2022) provide a comparison of the spatiotemporal evolution of antibubbles by using both the CHNS framework and a VOF DNS.

Figure 9.

Illustrative pseudocolour plots of (a) and (c) the field $\phi$ and the associated plots of the vorticity $\omega$ showing the spatiotemporal evolution of antibubbles for low (a) and (b) and high (c) and (d) values of $R_0$ and $R_1$ , the initial outer and inner radii of the shell of the antibubble (see the text and Pal et al. (Reference Pal, Ramadugu, Perlekar and Pandit2022) for details). We thank N. Pal for these figures.

6.5. Low-Reynolds-number interface-induced turbulence in the CHNS system

Emergent turbulence-type states have been found at low Reynolds numbers in a variety of systems, for instance, in fluids with polymer additives (see, e.g., Groisman & Steinberg Reference Groisman and Steinberg2000; Majumdar & Sood Reference Majumdar and Sood2011; Gupta & Pandit Reference Gupta and Pandit2017; Benzi & Ching Reference Benzi and Ching2018; Steinberg Reference Steinberg2021; Singh et al. Reference Singh, Perlekar, Mitra and Rosti2024) and in dense bacterial suspensions (see, e.g., Wensink et al. Reference Wensink, Dunkel, Heidenreich, Drescher, Goldstein, Löwen and Yeomans2012; Dunkel et al. Reference Dunkel, Heidenreich, Bär and Goldstein2013; Bratanov, Jenko & Frey Reference Bratanov, Jenko and Frey2015; Linkmann et al. Reference Linkmann, Boffetta, Marchetti and Eckhardt2019; Alert, Casademunt & Joanny 2021; Kiran et al. Reference Kiran, Gupta, Verma and Pandit2023; Mukherjee et al. Reference Mukherjee, Singh, James and Ray2023). In the former, this turbulence is driven by an increase in the Weissenberg number, whereas in the latter, it is obtained by increasing the activity of the bacterial suspension. Recently Padhan et al. (Reference Padhan, Vincenzi and Pandit2024 b) have demonstrated that low-Re turbulence occurs in the 2-D CHNS system if we increase the Weber number We or the Capillary number Ca (i.e. we decrease the interfacial tension; see table 2) and hence enhance interfacial fluctuations.

To obtain such interface-induced turbulence, Padhan et al. (Reference Padhan, Vincenzi and Pandit2024 b) begin with a periodic arrangement of vortices and antivortices, which is referred to as a vortex crystal or a cellular flow. Such cellular flows, imposed by a spatially periodic forcing with an amplitude $f_0$ and wavenumber $k_f$ , can be disordered via turbulence as shown for 2-D fluid turbulence by Perlekar & Pandit (Reference Perlekar and Pandit2010) and for 2-D fluids with polymer additives by Gupta & Pandit (Reference Gupta and Pandit2017) and Plan et al. (Reference Plan, Gupta, Vincenzi and Gibbon2017). For the 2-D CHNS system, we follow the discussion of Padhan et al. (Reference Padhan, Vincenzi and Pandit2024b ) and examine the dependence of the statistically steady state of the system as we increase the Capillary number Ca. The natural length, time, and velocity scales are, respectively, $k_f^{-1}$ , $\nu k_f/f_0$ and $U \equiv f_0/(\nu k_f^2)$ .

Figure 10.

(a–e) The scaled total kinetic energy $e(t)/e_0$ versus time $t$ [here, $e_0 = U^2$ and $U = f_0/(\nu k_f^2)$ ]; (f–j) the corresponding power spectra $\tilde {e}(f)$ of $e(t)/e_0$ . Pseudocolour plots, at a representative time, of (k–o) the vorticity $\omega$ overlaid with the $\phi = 0$ contour (black lines) and (p–t) the energy spectra $\mathcal {E}(k_x, k_y)$ (see Padhan, Vincenzi & Pandit Reference Padhan, Vincenzi and Pandit2024 b for details); the capillary number $Ca$ increases from $0.01$ in the bottom row to $0.2$ in the top row.

In figures 10(a) and 10(e) we plot the scaled total kinetic energy $e(t)/e_0$ versus time $t$ , with $e_0 = U^2$ ; figures 10(b) and 10(f) display the corresponding power spectra $|\tilde {e}(f)|$ of $e(t)/e_0$ . Figures 10(c) and 10(g) depict pseudocolour plots of the vorticity $\omega$ overlaid with the $\phi = 0$ contour (black lines) at a representative time. In figures 10(d) and 10(h) we present pseudocolour plots of the energy spectra $\mathcal {E}(k_x, k_y)$ at a representative time; these spectra are not averaged over wavenumber shells because the underlying crystalline state is not isotropic. Figures 10(a–e) and 10(f–j) help us to distinguish between states that are temporally periodic ( $Ca=0.17$ ) from those that are chaotic ( $Ca=0.19$ ); figures 10(k–o) and 10 ( p–t) aid the identification of the spatial order of the state. If the state is periodic in space, $\mathcal {E}(k_x, k_y)$ shows Bragg peaks in the reciprocal lattice of the vortex crystal (e.g., the strong red peaks in figure 10 h), where we use the standard terminology of crystal physics. We find a rich variety of states: STPO, spatially and temporally periodic (i.e. a spatiotemporal crystal of the type discussed in Perlekar & Pandit (Reference Perlekar and Pandit2010) and Gupta & Pandit (Reference Gupta and Pandit2017)); STPOG, is like STPO, but with a grain boundary separating two crystalline parts; STC denotes spatiotemporal chaos. Here, TPO denotes temporally periodic oscillations. For similar states in studies of turbulence-induced melting of a vortex crystals see Perlekar & Pandit (Reference Perlekar and Pandit2010) for fluids and Gupta & Pandit (Reference Gupta and Pandit2017) for fluids with polymer additives. In the latter case, the melting of the vortex crystal can be induced by elastic turbulence at low $Re$ ; this is akin to the low- $Re$ interface-fluctuation-induced melting we discuss here. In particular, the study of Gupta & Pandit (Reference Gupta and Pandit2017) has noted that the boundary between spatially and temporally periodic states and ones that exhibit STC is non-monotonic, i.e. there is re-entrance from one state into another and back. Interface-induced turbulence in the 2-D CHNS system also displays such re-entrance, as a function of $Ca$ ; Padhan et al. (Reference Padhan, Vincenzi and Pandit2024b ) have found the following re-entrant sequence of states: chaotic $\to$ temporally periodic $\to$ chaotic.

7.1. Ternary-fluid CHNS

Figure 11.

Pseudocolour plot of $F(\{c_i\})$ projected onto a Gibbs triangle for the CHNS3 model (7.1). The three vertices yield the three minima of $F(\{c_i\})$ : the top vertex is $(c_1, c_2, c_3) = (1, 0, 0)$ ; the left-hand vertex is $(c_1, c_2, c_3) = (0, 1, 0)$ ; the right-hand vertex is $(c_1, c_2, c_3) = (0, 0, 1)$ .

The CHNS3 model for a ternary-fluid mixture uses the following variational free energy in the domain $\Omega$ (see, e.g., Boyer & Lapuerta Reference Boyer and Lapuerta2006; Kim Reference Kim2007; Boyer et al. Reference Boyer, Lapuerta, Minjeaud, Piar and Quintard2010; Shek & Kusumaatmaja Reference Shek and Kusumaatmaja2022a ; Padhan & Pandit Reference Padhan and Pandit2023b ):

(7.1) \begin{eqnarray} \mathcal {F}(\{c_i, \nabla c_i\}) = \int \nolimits _{\Omega } {\textrm{d}}\Omega \left [\frac {12}{\epsilon }F(\{c_i\}) + \frac {3\epsilon }{8} \sum _{i=1}^{3}\gamma _i(\nabla c_{i})^2\right ], \end{eqnarray}

where the concentration fields $c_i (i = 1, 2, 3)$ are conserved order parameters that satisfy the constraint $\sum _{i=1}^{3} c_i = 1$ , $\epsilon$ is the thickness of the interface, the variational bulk free energy

(7.2) \begin{eqnarray} F(\{c_i\}) = \displaystyle \sum _{i=1}^{3} \gamma _i c_{i}^2(1-c_{i})^2, \end{eqnarray}

and the gradient terms give the surface-tension penalties for interfaces, with

(7.3) \begin{eqnarray} \sigma _{ij} = (\gamma _i + \gamma _j) / 2 \end{eqnarray}

the bare surface (or interfacial) tension for the interface between the phases $i$ and $j$ ; the equilibrium values of $c_i$ follow from the global minima of $F(\{c_i\})$ . In figure 11, we show a pseudocolour plot of $F(\{c_i\})$ projected onto a Gibbs triangle (see, e.g., Kim Reference Kim2007). The equilibrium chemical potential of fluid- $i$ is $\mu _i = ({\delta \mathcal F}/{\delta c_i}) + \beta (c_i)$ , with $\beta (c_i)$ the Lagrange multiplier that ensures $\displaystyle \sum _{i=1}^3 c_i = 1$ , whence we get

(7.4) \begin{align} \mu _i & = -\frac {3}{4}\epsilon \gamma _i \nabla ^2 c_i + \frac {12}{\epsilon } \left [ \vphantom{\frac{1}{2}} \gamma _i c_i(1-c_i)(1-2c_i) \right. \nonumber \\ & \left. - \frac {6\gamma _1 \gamma _2 \gamma _3 (c_1 c_2 c_3)}{\gamma _1 \gamma _2 + \gamma _1 \gamma _3 + \gamma _2 \gamma _3} \right]\,. \end{align}

The incompressible CHNS3 equations can be written as follows (see, e.g., Kim & Lowengrub Reference Kim and Lowengrub2005; Boyer et al. Reference Boyer, Lapuerta, Minjeaud, Piar and Quintard2010; Tóth, Zarifi & Kvamme Reference Tóth, Zarifi and Kvamme2016; Padhan & Pandit Reference Padhan and Pandit2023b ):

(7.5) \begin{align} \rho (\{c_i\}) (\partial _t {\boldsymbol {u}} + (\boldsymbol {u} \cdot \nabla ) \boldsymbol {u}) &= -\nabla P + \nabla \cdot \left [\eta (\{c_i\}) (\nabla \boldsymbol {u} + \nabla \boldsymbol {u}^{T})\right ] \nonumber \\ &\quad - \displaystyle \sum _{i=1}^{3}c_i \nabla \mu _i + \rho (\{c_i\}) \boldsymbol {g} - \alpha \boldsymbol {u} + \boldsymbol {f}^{ext}; \nonumber \\ \nabla \cdot \boldsymbol {u} &= 0; \nonumber \\ \partial _t{c_{i}} + (\boldsymbol {u}.{\nabla })c_{i} &= \frac {M}{\gamma _{i}}{\nabla }^2 \mu _{i}, \,\,\,\, i = 1\,\, \textrm {or}\,\, 2 ; \nonumber \\ \eta (\{c_i\}) &= \displaystyle \sum _{i=1}^{3}\eta _i c_i; \nonumber \\ \rho (\{c_i\}) &= \displaystyle \sum _{i=1}^{3}\rho _i c_i. \end{align}

Here, $\boldsymbol {g}$ , $\boldsymbol {f}^{ext}$ and $\alpha$ are, respectively, the acceleration because of gravity, an external forcing, and the coefficient of friction. $\eta _i$ and $\rho _i$ are the viscosity and density of fluid $i$ , respectively. The CHNS3 model becomes the binary-fluid CHNS model for $(c_1, c_2, c_3) = (c, 0, 0)$ .

7.1.1. Boussinesq approximation

If we use the Boussinesq approximation, the CHNS3 equations can be written as follows:

(7.6) \begin{eqnarray} \partial _t {\boldsymbol {u}} + (\boldsymbol {u} \cdot \boldsymbol {\nabla }) \boldsymbol {u} &=& - \frac {1}{\rho _0}\boldsymbol {\nabla } P + \nu \nabla ^2 {\boldsymbol {u}}- \frac {1}{\rho _0}\sum _{i=1}^{3}(c_i \boldsymbol {\nabla } \mu _i) + \frac {[\rho (\{c_i\}) - {\rho }_0]}{{\rho }_0}\boldsymbol {g} - \alpha \boldsymbol {u},\,\, \end{eqnarray}

(7.7) \begin{eqnarray} \nabla \cdot \boldsymbol {u} &=& 0. \end{eqnarray}

We write the density in the form

(7.8) \begin{eqnarray} \rho (\{c_i\}) &=& \sum \limits _{i=1}^{3} \rho _i c_i \nonumber \\ &=& \rho _1 c_1 + \rho _2 c_2 + \rho _3 (1-c_1-c_2) \nonumber \\ &=& \rho _3 + (\rho _1 - \rho _3) c_1 + (\rho _2 - \rho _3) c_2; \end{eqnarray}

we use $\rho _0 = \rho _3$ , so

(7.9) \begin{eqnarray} \frac {[\rho (\{c_i\}) - {\rho }_0]}{{\rho }_0} &=& \left(\frac {\rho _1-\rho _3}{\rho _0} \right) c_1 + \left(\frac {\rho _2-\rho _3}{\rho _0}\right) c_2 \nonumber \\ &=& \mathcal A_1 c_1 + \mathcal A_2 c_2, \end{eqnarray}

where $\mathcal A_1$ and $\mathcal A_2$ are the Atwood numbers and $\mathcal A_1, A_2 \ll 1$ .

Figure 12.

Coarsening in a ternary-fluid mixture from our DNS of the 2-D CHNS3 equations. (a) Pseudocolour plots of the phase-fields $c_2 - c_1$ at simulation times $t = 0.8, 2, 7, 13$ , increasing from left to right. (b) Pseudocolour plots of the corresponding vorticity fields $\omega (\boldsymbol {x}, t)$ (normalised by their absolute maximum for ease of visualisation). (c) Log–log plots versus $t$ of the scaled lengths $L_1(t)/L_0$ and ${L}_2(t)/L_0$ , with $L_0 = 2\pi$ the side of the simulation domain. (d) The energy spectrum $E(k,t)$ at simulation times $t = 2, 7, 13$ . Simulation parameters: viscosity $\nu = 10^{-3}$ ; grid points $256 \times 256$ ; surface tension coefficients $(\sigma _{12}, \sigma _{23}, \sigma _{13}) = (1, 1, 1)$ .

7.2. Phase separation in the ternary-fluid CHNS3

We turn now to the phase separation of ternary-fluid mixtures that occurs in a variety of settings (see, e.g., Huang, de La Cruz & Swift Reference Huang, de La Cruz and Swift1995; Singh et al. Reference Singh and Puri2015 a; Wang et al. Reference Wang, Altschuh, Ratke, Zhang, Selzer and Nestler2019; Shek & Kusumaatmaja Reference Shek and Kusumaatmaja2022b ). In figure 12, we show illustrative results for coarsening in such a mixture, from our DNS of the 2-D CHNS3 equations, which we visualise via pseudocolour plots of the difference of the phase-fields $(c_2 - c_1)$ (figure 12 a) and the corresponding vorticity fields $\omega (\boldsymbol {x}, t)$ , normalised by their absolute maximum (figure 12 b), at simulation times $t = 0.8, 2, 7, 13$ , which increase from left to right. We define the following scaled lengths (cf. (6.6) for the two-fluid CHNS):

(7.10) \begin{equation} L_1(t) = 2\pi \tfrac {\displaystyle \sum _k S_1(k, t)}{\displaystyle \sum _k k S_1(k, t)};\qquad L_2(t) = 2\pi \tfrac {\displaystyle \sum _k S_2(k, t)}{\displaystyle \sum _k k S_2(k, t)}. \end{equation}

The scaled lengths $L_1(t)/L_0$ and ${L}_2(t)/L_0$ , with $L_0 = 2\pi$ , increase with time, in a manner that is consistent with $\sim t^{2/3}$ , the power-law growth for the inertia-dominated regime (see the log–log plot in figure 12 c). This power-law growth has also been seen in the molecular-dynamics simulations of Singh & Puri (Reference Singh and Puri2015a ). We expect that such ternary-fluid phase separation will also be suppressed by turbulence in the CHNS3 equations, just as turbulence leads to coarsening arrest in the binary-fluid case (see, e.g., Perlekar et al. Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014, Reference Perlekar, Pal and Pandit2017), but we must account for additional non-dimensional control parameters such as the ratios of the surface tensions between the three fluid phases that coexist in equilibrium. To the best of our knowledge, a complete study of such turbulence-induced coarsening arrest in three-phase fluid mixtures has not been carried out so far.

7.3. Spatiotemporal evolution of droplets in turbulent flows

We follow and generalise, to the case of compound droplets, the investigation of the spatiotemporal evolution of droplets in binary-fluid turbulent flows by Pal et al. (Reference Pal, Perlekar, Gupta and Pandit2016), who studied the advection of a droplet, initially circular with a diameter $d_0$ , by a turbulent binary-fluid flow for which they used the CHNS system in two dimensions. The droplet was active, insofar as it affected the flow, and was, in turn, deformed by the flow. Pal et al. (Reference Pal, Perlekar, Gupta and Pandit2016) obtained the acceleration components of the droplet centre of mass and showed that their probability distribution function (PDF) had wide, non-Gaussian tails. They uncovered multifractal fluctuations in the time series of the scaled perimeter (see below) of the droplet. Finally they showed that the droplet fluctuations led to an enhancement of the energy spectrum $E(k)$ for large wave numbers $k$ , and thence to dissipation reduction, as in fluid turbulence with polymer additives (see, e.g., Perlekar, Mitra & Pandit Reference Perlekar, Mitra and Pandit2006; Perlekar, Mitra & Pandit Reference Perlekar, Mitra and Pandit2010; Gupta, Perlekar & Pandit Reference Gupta, Perlekar and Pandit2015). We refer the reader to Pal et al. (Reference Pal, Perlekar, Gupta and Pandit2016) for the details of their DNS. Our investigations of a compound droplet are motivated by studies of compound droplets in external electric fields (see Santra et al. Reference Santra, Panigrahi, Das and Chakraborty2020), the examination of the breakup of double-emulsion droplets in linear flows (see Stone & Leal Reference Stone and Leal1990), the numerical and theoretical investigations of the dynamics of a compound vesicle (a lipid bilayer membrane enclosing a fluid with a suspended particle) in a shear flow (see Veerapaneni et al. Reference Veerapaneni, Young, Vlahovska and Bławzdziewicz2011) and a compound vesicle in shear flow (see Sinha & Thaokar Reference Sinha and Thaokar2019). Experiments have also suggested that a compound droplet can be used as a model for a white-blood cell (see Levant & Steinberg Reference Levant and Steinberg2014). Furthermore, the deformation and breakup of a compound droplet has been studied in a 3-D oscillatory shear flow (see Liu et al. Reference Liu, Lu, Li, Yu and Sahu2021) and in a channel flow (see Lanjewar & Ramji Reference Lanjewar and Ramji2024).

Figure 13.

The pseudocolour plots of (a) the phase field $c_2$ for a simple droplet, and the phase field $c_2 - c_1$ for the compound droplets with radius ratios (b) $R_{in}/R_{out} = 0.5$ and (c) $R_{in}/R_{out} = 0.9$ .

We now give illustrative results for the turbulent advection of a simple droplet, initially circular with a radius $R_0$ , and for a compound droplet, initially concentric circles with inner and outer radii $R_{in}$ and $R_{out}$ ; for the former we employ the binary-fluid CHNS in two dimensions; and for the latter we use the 2-D ternary-fluid CHNS; in both cases we impose periodic BCs, obtain statistically steady turbulence via a Kolmogorov-type forcing with amplitude $F_{0}$ , as in Pal et al. (Reference Pal, Perlekar, Gupta and Pandit2016), and a forcing wavenumber $k_f$ that yields an energy spectrum in which the forward-cascade regime is dominant. The pseudocolour plots in figure 13(a) show the phase field $\phi$ , for a simple droplet, at different representative times, which increase from left to right. For the compound droplet we present pseudocolour plots of $(c_2 - c_1)$ with the radius ratios (b) $R_{in}/R_{out} = 0.5$ (figure 13 b,c) $R_{in}/R_{out} = 0.9$ (figure 13 c).

To characterise droplet fluctuations, we use the time dependence of the deformation parameter defined in Pal et al. (Reference Pal, Perlekar, Gupta and Pandit2016),

(7.11) \begin{equation} \Gamma (t)=\frac {{\mathcal {S}}(t)}{{\mathcal {S}}_{0}(t)}-1, \end{equation}

with ${\mathcal {S}}(t)$ , the perimeter of the droplet, and ${\mathcal {S}}_0(t)$ , the perimeter of an undeformed droplet. In the binary-fluid case we track the perimeter via the $\phi =0$ contour; in the ternary-fluid case we use the contours $c_1 = 0.5$ and $c_2 = 0.5$ for the perimeters of the inner and outer droplets, respectively. The strength of the non-dimensionalised forcing is given by the Grashof number $Gr=L^{4}F_{0}/\nu ^{2}$ ; and the forcing-scale Weber number $We\equiv \rho k_f^{-3}F_0/\sigma$ measures the non-dimensionalised inverse surface tension; these parameters are chosen such that the droplets are not torn asunder during our DNSs. In figure 14(a) we plot $\Gamma (t)$ versus time (scaled by the forcing time scale $T = \nu k_f/F_0$ ) for a simple droplet (red line) and compound droplets with radius ratios $R_{in}/R_{out} = 0.5$ (green line) and $R_{in}/R_{out} = 0.9$ (blue line). In figures 14(b) and 14(c) we plot, respectively, the PDF of $\Gamma$ and the multifractal spectrum $D(h)$ of the time series $\Gamma (t)$ . From figures 14(a) and 14(c) we clearly see that the temporal fluctuations $\Gamma$ , its PDF and its multifractal spectrum $D(h)$ are similar for simple droplet interface and outer droplet interface in the compound droplet with a radius ratio $R_{in}/R_{out} = 0.5$ . By contrast, all these measures are reduced for the perimeter of the outer droplet interface in a compound droplet with a radius ratio $R_{in}/R_{out} = 0.9$ because the fluctuations of the outer interface are constrained by the presence of the inner droplet.

The energy spectra $E(k)$ for a single-phase turbulent fluid, a binary fluid with a simple droplet and a ternary fluid with a compound droplet (with radius ratios $R_{in}/R_{out} = 0.5$ and $R_{in}/R_{out} = 0.9$ ) are given in figure 14(d). These spectra show that droplet fluctuations, of both single and compound droplets, enhance $E(k)$ at large $k$ , relative to its counterpart for single-fluid turbulence. This enhancement is reminiscent of a similar enhancement of $E(k)$ in fluid turbulence by polymer additives, where it is associated with dissipation reduction can be understood as a $k$ -dependent correction to the viscosity (see Perlekar et al. Reference Perlekar, Mitra and Pandit2010; Gupta et al. Reference Gupta, Perlekar and Pandit2015). Such a scale-dependent correction to the viscosity also occurs for a single droplet in a turbulent flow as has been discussed by Pal et al. (Reference Pal, Perlekar, Gupta and Pandit2016).

Figure 14.

(a) Plot versus time of the deformation parameter $\Gamma (t)$ for a simple droplet (red line) and compound droplets with radius ratios $R_{in}/R_{out} = 0.5$ (green line) and $R_{in}/R_{out} = 0.9$ (blue line); the horizontal axis is scaled by the forcing time scale $T = \nu k_f/f_0$ . The corresponding (b) PDFs of $\Gamma$ (semi-log plots) and (c) the multifractal spectra $D(h)$ of $\Gamma$ . (d) Log–log plots of the energy spectrum $E(k)$ for a single-phase turbulent fluid, a binary fluid with a simple droplet, and a ternary fluid with a compound droplet (with radius ratios $R_{in}/R_{out} = 0.5$ and $R_{in}/R_{out} = 0.9$ ); the inset shows a pseudocolour plot of the vorticity in the presence of a compound droplet.

7.4. Bubble passing through an interface between two fluids

Figure 15.

A bubble passing through a fluid–fluid interface. (a) The isosurface plot of the $(c_1, c_2)$ fields. (b) The isocontour plot of the $z$ component of the vorticity field. (c) The kinetic energy time series of the droplet (blue line), where $e(t) = \langle |\boldsymbol {u}(\boldsymbol {x})|^2\rangle _{\boldsymbol {x}\in (c_1 \geqslant 0.5)}$ . The red line shows line shows the kinetic energy time series of the fluid–fluid interface defined as $e(t) = \langle |\boldsymbol {u}(\boldsymbol {x})|^2\rangle _{\boldsymbol {x} \in (0.1 \leqslant c_2\leqslant 0.9)}$ . We use the characteristic velocity and time scales as $U_0 = g\epsilon ^2/\nu$ and $T_0 = {\nu }/{g\epsilon }$ .

How do bubbles or droplets pass through an interface between two fluids? This problem has attracted considerable attention (see, e.g., Manga & Stone Reference Manga and Stone1995; Dietrich et al. Reference Dietrich, Poncin, Pheulpin and Li2008; Bonhomme et al. Reference Bonhomme, Magnaudet, Duval and Piar2012; Li et al. Reference Li, Al-Otaibi, Vakarelski and Thoroddsen2014; Natsui et al. Reference Natsui, Takai, Kumagai, Kikuchi and Suzuki2014; Singh & Bart Reference Singh and Bart2015; Feng et al. Reference Feng, Muradoglu, Kim, Ault and Stone2016; Prosperetti Reference Prosperetti2017; Emery, Raghupathi & Kandlikar Reference Emery, Raghupathi and Kandlikar2018; Emery & Kandlikar Reference Emery and Kandlikar2019; Kumar, Rohilla & Das Reference Kumar, Rohilla and Das2019; Choi & Park Reference Choi and Park2021; Chowdhury, Mahapatra & Sen Reference Chowdhury, Mahapatra and Sen2022; Rabbani & Ray Reference Rabbani and Ray2024) in the fluid dynamics, chemical engineering and microfluidics literature. Modern experiments that use high-speed cameras to track the passage of bubbles through fluid–fluid interfaces have led to theoretical and numerical investigations of this problem.

In figure 15 we present our results from an illustrative DNS of a bubble passing through a fluid–fluid interface; for this we employ the three-component 3-D CHNS3 equations (7.7)–(7.9) and choose parameters such that the bubble does go through the interface and does not get trapped there. For such studies it is natural to use the characteristic velocity and time scales $U_0 = g\epsilon ^2/\nu$ and $T_0 = {\nu }/{g\epsilon }$ , an elongated simulation domain ( $(L_x, L_y, L_z) = (4\pi, \pi, \pi )$ , with $512 \times 128 \times 128$ grid points) with periodic BCs in the directions normal to gravity, and volume penalisation in the direction of gravity to incorporate solid boundaries. We use six grid points at both top and bottom boundaries in our volume-penalisation scheme and the following simulation parameters: $\nu = 3.5\times 10^{-3}$ , $g=1$ , $\mathcal A_1 = 0.132$ , $\mathcal A_2 = 0.132$ , $\rho _1 = 1.132, \rho _2 = 0.868, \rho _3 = 1$ (with $\rho _3$ as the reference density), and $\sigma _{12} = \sigma _{13}= 0.5$ and $\sigma _{23} = 0.01$ . Figures 15(a) and 15(b) show, respectively, isosurface plots of the $c_1$ and $c_2$ fields and isocontour plots of the $z$ component of the vorticity field, which illustrate how a $c_1$ droplet of the first fluid (in red), with concentration field $c_1$ , passes through the $c_2-c_3$ interface (light brown) between the second and third fluids. As it passes through the interface, this droplet entrains some of the heavy fluid, which forms a slender neck that collapses eventually to yield droplets of the heavy fluid that fall back onto the interface (see, e.g., Singh & Bart Reference Singh and Bart2015 and Emery et al. Reference Emery, Raghupathi and Kandlikar2018). In figure 15(c) we show that we can track the passage of this bubble through the interface by monitoring the temporal evolution of the droplet’s kinetic energy [ $e_d(t) \equiv \langle |\boldsymbol {u}(\boldsymbol {x})|^2\rangle _{\boldsymbol {x}\in (c_1 \geqslant 0.5)}$ (blue line)] and the energy of the $c_2-c_3$ interface ( $e_I(t) = \langle |\boldsymbol {u}(\boldsymbol {x})|^2\rangle _{\boldsymbol {x} \in (0.1 \leqslant c_2\leqslant 0.9)}$ (red line)); both these quantities display maxima when the bubble passes through the interface.

7.5. The coalescence of liquid lenses and droplets

The coalescence of liquid droplets and lenses is a problem of fundamental importance in fluid mechanics and statistical mechanics (see, e.g., Paulsen et al. Reference Paulsen, Carmigniani, Kannan, Burton and Nagel2014; Pal Reference Pal2016; Heinen et al. Reference Heinen, Hoffmann, Diewald, Seckler, Langenbach and Vrabec2022; Padhan & Pandit Reference Padhan and Pandit2023b ; Scheel et al. Reference Scheel, Xie, Sega and Harting2023). Our DNS for the CHNS3 (Padhan & Pandit Reference Padhan and Pandit2023b ) model can be used to examine the development of liquid-lens mergers in phase-separated ternary-fluid systems as we summarise below. The coexistence of three immiscible fluids leads to three distinct interfaces with three interfacial tensions: $\sigma _{ij}$ is the surface-tension coefficient for the $ij$ interface, where the integers $i$ and $j$ ( $= 1, 2$ or $3$ ) label the coexisting phases. We prepare neutrally buoyant, symmetrical or asymmetrical, lenses in two dimensions by starting our DNS with the following configuration for a single circular droplet of fluid 1, with radius $R_0$ and centre at $(\pi, \pi )$ , placed at the interface between fluids 2 and 3 (see figure 16 a):

(7.12) \begin{eqnarray} c_1(x,y,0) &=& \frac {1}{2}\left [1-\tanh \left (\frac {\sqrt {(x-\pi )^2+(y-\pi )^2}-R_{0}}{2 \sqrt {2} \epsilon }\right )\right ];\nonumber \\ c_2(x, y, 0) &=& \frac {1}{2}\left [1-\tanh \left (\frac {y-\pi }{2 \sqrt {2} \epsilon }\right )\right ] - c_1(x,y,0); \\ \omega (x, y, 0) &=& 0. \nonumber \end{eqnarray}

Figure 16.

Plots showing the three coexisting phases 1 (blue), 2 (red) and 3 (green) in regions $\Omega _1$ , $\Omega _2$ and $\Omega _3$ , respectively. (a) The initial profile with a circular droplet of radius $R_0/L = 0.15$ of fluid 1 at the interface between fluids 2 and 3. ( $L = 2\pi$ is the side length of the simulation domain.) (b) The equilibrium profile of the droplet is a symmetrical lens, because we choose $(\sigma _{12}, \sigma _{13}, \sigma _{23}) = (1, 1, 1)$ ; the three contact angles are $\theta _1$ , $\theta _2$ and $\theta _3$ . (c) The equilibrium profile of the droplet is an asymmetrical lens if we choose $(\sigma _{12}, \sigma _{13}, \sigma _{23}) = (1.4, 0.8, 1)$ .

The initial and equilibrium configurations in three dimensions are a sphere and a lenticular biconvex lens, respectively. As time evolves in our (2-D) DNS, the initial circular droplet relaxes to its equilibrium-lens shape as shown in figures 16(b) and 16(c) for $(\sigma _{12},\sigma _{13},\sigma _{23}) = (1, 1, 1)$ and $(\sigma _{12},\sigma _{13},\sigma _{23}) = (1.4, 0.8, 1)$ , respectively. Now we verify the Young relations for liquid lenses at equilibrium (see, e.g., Boyer & Lapuerta Reference Boyer and Lapuerta2006; McHale et al. Reference McHale, Afify, Armstrong, Wells and Ledesma-Aguilar2022). The theoretical distance $d^{th}$ between the two triple-phase junctions in figures 16(b) and 16(c) are given by the following Young relations:

(7.13) \begin{eqnarray} d^{th} &=& (l_1 + l_2)^{-{\frac {1}{2}}};\nonumber \\ l_1 &=& \frac {2(\pi -\theta _3) - \sin (2(\pi -\theta _3))}{8A\sin (\pi -\theta _3)};\\ l_2 &=& \frac {2(\pi -\theta _1) - \sin (2(\pi -\theta _1))}{8A\sin (\pi -\theta _1)}.\nonumber \end{eqnarray}

The contact angles are related as follows:

(7.14) \begin{eqnarray} \frac {\sin \theta _1}{\sigma _{23}} &=& \frac {\sin \theta _2}{\sigma _{13}} = \frac {\sin \theta _3}{\sigma _{12}}. \end{eqnarray}

We calculate the distance between the triple-phase junctions $d^{sim}$ from our simulations (see figures 16 b and 16 c) and compare them with $d^{th}$ in table 3. The agreement between these values is good. Furthermore, we can compare the results of our DNS with the prediction of the Laplace law for pressure jumps, at equilibrium; these jumps are defined as follows:

Table 3.

The distance between two triple-phase junctions calculated from theory $d^{th}$ (see (7.14)] and numerical simulations $d^{sim}$ for the symmetrical lens (figure 16 b) and the asymmentrical lens (figure 16 c).

Figure 17.

Illustrations of the CHT (see the text) that we use to fit circles for the lens interfaces for (a) $(\sigma _{12}, \sigma _{13}, \sigma _{23}) = (1, 1, 1)$ (run $\mathcal {R}1$ ) and (b) $(\sigma _{12}, \sigma _{13}, \sigma _{23}) = (1.4, 0.8, 1)$ (run $\mathcal {R}2$ ).

Figure 18.

Pseudocolour plot of illustrative equilibrium-pressure profiles for (a) $(\sigma _{12}, \sigma _{23}, \sigma _{13}) = (0.6, 1, 0.8)$ (run RN3) and (b) $(\sigma _{12}, \sigma _{23}, \sigma _{13}) = (1.4, 1, 0.6)$ (run RN4).

(7.15) \begin{eqnarray} \frac {\sigma _{13}}{R_{13}} &=& P_1 - P_3 = P_1 - P_2 = \frac {\sigma _{12}}{R_{12}};\nonumber \\ P_2 - P_3 &=& 0; \end{eqnarray}

where $R_{12}$ and $R_{13}$ are the radii of curvatures of the interfaces between fluids 1 and 2 and between fluids 1 and 3, respectively. To calculate the radii of curvatures, we use the circle Hough transform (CHT) in MATLAB (see Atherton & Kerbyson Reference Atherton and Kerbyson1999). In figure 17, we illustrate the CHT of the images in figures 16(b) and 16(c); red circles are the best-fit circles to the curves. The CHT gives the coordinates of the centres and the radii of the red circles. We calculate the theoretical values of the pressure jumps from these radii of curvatures. To compare our DNS results with these theoretical values, we evaluate the pressure jumps from the following pressure-Poisson equation:

(7.16) \begin{eqnarray} \nabla ^2 P = \boldsymbol {\nabla } \cdot \left (-{\sum _{i=1}^{3}c_{i} \boldsymbol {\nabla } \mu _{i}}\right ), \end{eqnarray}

which we portray via pseudocolour plots of illustrative equilibrium-pressure profiles for runs RN1 and RN2 in figure 18. We present our DNS results in table 4 for various runs. These results are in good agreement with their theoretical counterparts. We obtain similar results from 3-D DNSs of the CHNS3 model: in figure 19(a) we give an isosurface plot, with $c_1 = 0.5$ , for a 3-D lenticular biconvex lens. We then calculate its Gaussian curvature $\kappa$ by implementing, in MATLAB, the algorithm described in Meyer et al. (Reference Meyer, Desbrun, Schröder and Barr2003). The isosurface plot of $\kappa$ is shown in figure 19(b); clearly, $\kappa$ is constant throughout the lens surface, except at the edges. So, in figure 19(c), we present the PDF of $\kappa$ to find out the most probable value of $\kappa$ . We consider the values of $\kappa$ with the highest probability, namely, $0.59, -0.04, 1.23$ ; the average value is $\kappa \simeq 0.593$ . The Gaussian curvature $\kappa = 1/R_{G}^2$ for a sphere of radius $R_G$ (see, e.g., Nothard et al. Reference Nothard, McKenzie, Haines and Jackson1996). The symmetric 3-D lens in figure 19(a) is a combination of two surfaces that are parts of spheres of equal radii. Then we follow the Laplace law in three dimensions to evaluate the pressure jumps,

(7.17) \begin{eqnarray} \frac {2\sigma _{13}}{R_{13}} &=& P_1 - P_3 = P_1 - P_2 = \frac {2\sigma _{12}}{R_{12}};\nonumber \\ P_2 - P_3 &=& 0; \end{eqnarray}

where $R_{13} = R_{12} = R_G = 1/\sqrt {\kappa }$ . The theoretical values of the pressure jumps are given in table 4 (see run RN-3D); we solve the pressure-Poisson equation (7.16) numerically; we find good agreement between the theoretical and numerical values of these jumps.

Table 4.

Illustrative comparisons of theoretical and our DNS results for Laplace-pressure jumps for different lens shapes, ${\Delta P} \equiv P_1 - P_2 = P_1 - P_3$ . There is good agreement between these results. While evaluating $\Delta P$ , we calculate the values of $P_1$ , $P_2$ and $P_3$ at points that are far from the interface.

Figure 19.

(a) The isosurface plot of $c_1 = 0.5$ , at equilibrium, illustrating a lens in three dimensions. (b) The isosurface plot of the corresponding Gaussian curvature $\kappa$ . (c) The PDF of the Gaussian curvature.

Figure 20.

Illustrative results from our DNSs of liquid-lens mergers in the CHNS3 model in (a,b) two dimensions and (c,d) three dimensions: pseudocolour plots of $\omega$ , with overlaid velocity vectors in (a) the viscous regime and (b) the inertial regime; the $c_1 = 0.5$ contour (magenta line) indicates the lens interface; $\omega$ is normalised by its maximal absolute value for ease of visualisation. Isosurface plots of $c_1$ (green) and $|\omega |$ (brown) for (c) the viscous regime and (d) the inertial regime.

We turn now to an overview of our recent study (Padhan & Pandit Reference Padhan and Pandit2023b ) that has shown how to use DNSs of the CHNS3 system (7.7)–(7.9) and their 2-D counterparts to study the spatiotemporal evolution of the merger of liquid lenses in both two dimensions and three dimensions. In figure 20 we present illustrative results from these DNSs in two dimensions (figure 20 a,b) and three dimensions (figure 20 c,d). In two dimensions we give pseudocolour plots of $\omega$ , with overlaid velocity vectors in (figure 20 a) the viscous regime and (figure 20 b) the inertial regime; the $c_1 = 0.5$ contour (magenta line) indicates the lens interface; $\omega$ is normalised by its maximal absolute value for ease of visualisation. Isosurface plots of $c_1$ (green) and $|\omega |$ (brown) for (figure 20 c) the viscous regime and (figure 20 d) the inertial regime. In the viscous regime, i.e. at large values of the Ohnesorge number $Oh$ , a vortex quadrupole dominates the flow in the region of the neck in both 2-D and 3-D lens mergers, shown in figures 20(a) and 20(c), respectively. In the small- $Oh$ inertial regime (figures 20 b and 20 d) for two dimensions and three dimensions, respectively, this quadrupole moves away from the region of the neck with the passage of time. We quantify the growth of the height $h(t)$ of the neck, in the low- $Oh$ case, in figure 21 that contains a log–log plot of $h$ versus the time $t$ . This plot shows clearly the crossover from the viscous regime with $h(t) \sim t$ , at early times, to $h(t) \sim t^{2/3}$ , at late times. The exponent $2/3$ is typical of inertial-regime neck growth; the early-time growth is similar to that found in the viscous case (because the early-time quadrupolar configuration is similar to that in the viscous case). The inset shows the velocity vectors near the neck region at a representative time. In figure 21(b) we show the profile of the $y$ -component of the velocity field $u_y(L/2, y)$ , along the vertical direction; the arrows show how the profiles flatten as $t$ increases.

Figure 21.

(a) Log–log plot of the neck height $h$ versus time $t$ . The axes are scaled by their respective viscous length scales for $Oh = 0.08$ . The plot shows the crossover in the neck growth from the viscous regime with $h(t) \sim t$ to the inertial regime with $h(t) \sim t^{2/3}$ . The inset shows the velocity vectors near the neck region at a representative time. (b) The profile of the $y$ component of the velocity field $u_y(L/2, y)$ along the vertical direction; the arrows show the direction of the evolution of the profiles.

7.6. Active CHNS model

We consider the following incompressible CHNS equations (also called active model H) to study active turbulence in systems of contractile swimmers (see, e.g., Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015; Padhan & Pandit Reference Padhan and Pandit2023a ; Cates & Nardini Reference Cates and Nardini2024; Padhan et al. Reference Padhan, Vincenzi and Pandit2024 a; Padhan & Voigt Reference Padhan and Voigt2025) in two spatial dimensions:

(7.18) \begin{eqnarray} \partial _t \phi + (\boldsymbol {u} \cdot \nabla ) \phi &=& M \nabla ^2\left (\frac {\delta \mathcal F}{\delta \phi }\right )\, ; \end{eqnarray}

(7.19) \begin{eqnarray} \partial _t \omega + (\boldsymbol {u} \cdot \nabla ) \omega &=& \nu \nabla ^2\omega + \frac {3}{2}\epsilon \nabla \times (\nabla \cdot \boldsymbol {\Sigma }^A) - \alpha \omega \, ;\,\, \end{eqnarray}

(7.20) \begin{eqnarray} \nabla \cdot \boldsymbol {u} &=& 0; \end{eqnarray}

where $\omega$ is the vorticity field; $\nu$ , $\alpha$ and $M$ are the kinematic viscosity, bottom friction and mobility, respectively. Here $\mathcal F$ is the LG variational free-energy functional

(7.21) \begin{eqnarray} \mathcal F[\phi, \nabla \phi ] = \int _{\Omega } \left [\frac {3}{16} \frac {\sigma }{\epsilon }(\phi ^2-1)^2 + \frac {3}{4} \sigma \epsilon |\nabla \phi |^2\right ], \end{eqnarray}

in which the first term is a double-well potential with minima at $\phi = \pm 1$ . The scalar order parameter $\phi$ is positive (negative) in regions where the microswimmer density is high (low); in the interfaces between these regions, $\phi$ varies smoothly, over a width $\epsilon$ . The free-energy penalty for an interface is given by the bare surface tension $\sigma$ . In the inherently non-equilibrium active model H all terms in the stress tensor do not follow from $\mathcal F$ . In particular, we must include the stress tensor $\boldsymbol {\Sigma }^A$ , which has the form of a nonlinear Burnett term and has the components (see, e.g., Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015; Das, Bhattacharjee & Kirkpatrick Reference Das, Bhattacharjee and Kirkpatrick2020; Bhattacharjee & Kirkpatrick Reference Bhattacharjee and Kirkpatrick2022; Padhan & Pandit Reference Padhan and Pandit2023a ; Padhan et al. Reference Padhan, Vincenzi and Pandit2024a )

(7.22) \begin{eqnarray} \Sigma ^{A}_{ij} = -\zeta \left [\partial _i \phi \partial _j \phi - \frac {\delta _{ij}}{2} |\nabla \phi |^2\right ], \end{eqnarray}

where $\zeta$ , the activity coefficient, can take both positive and negative values: $\zeta \lt 0$ ( $\zeta \gt 0$ ) for contractile (extensile) swimmers. We emphasise that the free-energy functional used in this active model is a mathematical construct without a direct physical origin. The activity is introduced phenomenologically. For instance, an active-stress term with an effective negative surface tension coefficient has been incorporated into the model to capture the coarsening-arrest mechanism in contractile systems (see, e.g., Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015; Cates & Nardini Reference Cates and Nardini2024).

7.7. Generalised active CHNS model for an active self-propelling droplet

To study active, self-propelling droplets, we follow Padhan & Pandit (Reference Padhan and Pandit2023a ) and use two scalar fields $\phi$ and $\psi$ , with $\psi$ an active scalar, in the active-matter sense (see, e.g., Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013). (The terminology used in active matter and conventional fluid dynamics differs slightly. In fluid dynamics, both $\phi$ and $\psi$ are considered active scalars because they influence the velocity field $\boldsymbol {u}$ . However, in active matter, only $\psi$ is regarded as active, whereas $\phi$ is not.) We employ the following free-energy functional,

(7.23) \begin{align} \mathcal {F}[\phi, \nabla \phi, \psi, \nabla \psi ] &= \int _{\Omega } \frac {3}{16}\left (\frac {\sigma _1}{\epsilon _1}(\phi ^2-1)^2 + \frac {\sigma _2}{\epsilon _2} (\psi ^2-1)^2 \right ) - \beta \phi \psi \nonumber \\ &\quad + \frac {3}{4}\left ( \sigma _1 \epsilon _1 |\nabla \phi |^2 + \sigma _2 \epsilon _2 |\nabla \psi |^2\right ) \textrm{d}\Omega , \end{align}

where $\Omega$ is the region we consider. This model allows for interfaces of $\phi$ and $\psi$ , with (bare) positive interfacial tensions $\sigma _1$ and $\sigma _2$ and widths $\epsilon _1$ and $\epsilon _2$ , respectively; the coupling constant $\beta \gt 0$ , so there is an attractive coupling between $\phi$ and $\psi$ . Experiments on active droplets are carried out confined planar domains, so we use the following generalisation of the 2-D active incompressible CHNS equations given in § 7.6:

(7.24) \begin{eqnarray} \partial _t \phi + (\boldsymbol {u} \cdot \nabla ) \phi &=& M_1 \nabla ^2 \left ( \frac {\delta \mathcal F}{\delta \phi }\right ); \end{eqnarray}

(7.25) \begin{eqnarray} \partial _t \psi + (\boldsymbol {u} \cdot \nabla ) \psi &=& M_2 \nabla ^2 \left ( \frac {\delta \mathcal F}{\delta \psi }\right ); \end{eqnarray}

(7.26) \begin{eqnarray} \partial _t \omega + (\boldsymbol {u} \cdot \nabla ) \omega &=& \nu \nabla ^2 \omega -\alpha \omega +[\nabla \times (\mathfrak {S}^{\phi } + \mathfrak {S}^{\psi })]; \end{eqnarray}

(7.27) \begin{eqnarray} \nabla \cdot \boldsymbol {u} &=& 0; \quad \omega = (\nabla \times \boldsymbol {u}); \end{eqnarray}

(7.28) \begin{eqnarray} \mathfrak {S}^{\phi } &=& -(3/2)\sigma _1 \epsilon _1 \nabla ^2 \phi \nabla \phi ; \end{eqnarray}

(7.29) \begin{eqnarray} \mathfrak {S}^{\psi } &=& -(3/2) \tilde \sigma _2 \epsilon _2 \nabla ^2 \psi \nabla \psi ; \end{eqnarray}

where, the vorticity, kinematic viscosity and the bottom friction are, respectively, $\omega$ , $\nu$ and $\alpha$ and we set the constant fluid density $\rho =1$ . We use the constant mobilities $M_1$ and $M_2$ for $\phi$ and $\psi$ , respectively; the $\phi$ interfacial stress $\mathfrak {S}^{\phi }$ (7.28) follows from $\mathcal F$ ; in contrast, for the active stress $\mathfrak {S}^{\psi }$ (7.29) from $\psi$ , we use the active-model-H formulation (see, e.g., Wittkowski et al. Reference Wittkowski, Tiribocchi, Stenhammar, Allen, Marenduzzo and Cates2014; Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015; Shaebani et al. Reference Shaebani, Wysocki, Winkler, Gompper and Rieger2020; Padhan & Pandit Reference Padhan and Pandit2023a ). It is important to note that (i) both $\omega$ and $[\nabla \times (\mathfrak {S}^{\phi } + \mathfrak {S}^{\psi })]$ are orthogonal to the 2-D plane and (ii) the mechanical surface tension $\tilde \sigma _2 \neq \sigma _2$ , which can be either negative or positive values, unlike $\sigma _1$ and $\sigma _2$ that are positive. For contractile swimmers $\tilde \sigma _2 \lt 0$ and for extensile swimmers $\tilde \sigma _2 \gt 0$ ; the former show arrested phase separation, whereas the latter display complete phase separation (see, e.g., Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015; Padhan & Pandit Reference Padhan and Pandit2023a ). We will show in § 7.9 that the activity

(7.30) \begin{equation} A= |\tilde \sigma _2|/\sigma _2, \end{equation}

is the most important control parameter here.

Figure 22.

(a–c) Pseudocolour plots of the active-scalar field $\phi$ , at three representative times, which increase from left to right, for the activity parameter $\zeta = 0.1$ (see (7.18)–(7.22)); (d–f) pseudocolour plots of the vorticity $\omega$ corresponding, respectively, to the pseudocolour plots $\phi$ in (a–c).

7.8. Active CHNS turbulence

Figure 23.

(a) Log–log plots of the energy spectra $E(k)$ versus the wavenumber $k$ , for the activities $|\zeta | = 0.1$ and $|\zeta | = 1.5$ in (7.18)–(7.22); power-law regimes in these spectra are consistent with the dashed line $E(k) \sim k^{-5/3}$ ; (b) plot of the root-mean-squared velocity $u_{rms}$ versus $|\zeta |$ ; (c) log–linear plots versus $|\zeta |$ of (c) the mean coarsening length scale $L_c \equiv \langle \mathcal L(t)_t\rangle$ and (d) the integral length scale $L_I$ .

Turbulence in active fluids, which include dense bacterial suspensions, has garnered considerable attention over the past decade (see, e.g., Wensink et al. Reference Wensink, Dunkel, Heidenreich, Drescher, Goldstein, Löwen and Yeomans2012; Dunkel et al. Reference Dunkel, Heidenreich, Bär and Goldstein2013; Bratanov et al. Reference Bratanov, Jenko and Frey2015; Alert et al. 2021). Many models of active fluids consider systems of polar active swimmers (e.g., Jain et al. Reference Jain, Rana, Ramaswamy and Perlekar2024; Rana et al. Reference Rana, Chatterjee, Ro, Levine, Ramaswamy and Perlekar2024) or Toner–Tu type systems and their generalisations (see, e.g., Toner & Tu Reference Toner and Tu1998; Toner, Tu & Ramaswamy Reference Toner, Tu and Ramaswamy2005; Rana & Perlekar Reference Rana and Perlekar2020; Alert et al. 2021; Mukherjee et al. Reference Mukherjee, Singh, James and Ray2021; Gibbon et al. Reference Gibbon, Kiran, Padhan and Pandit2023; Kiran et al. Reference Kiran, Gupta, Verma and Pandit2023, 2024). In a recent paper, Padhan et al. (Reference Padhan, Vincenzi and Pandit2024a ) have demonstrated that a new type of active-scalar turbulence occurs in active-model H (see (7.18)–(7.22)), whose stochastic version has been studied in the context of MIPS that has been discussed at very low Reynolds numbers by Tiribocchi et al. (Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015) and Cates & Tailleur (Reference Cates and Tailleur2015). We give an overview of the work of Padhan et al. (Reference Padhan, Vincenzi and Pandit2024a ) in figures 22 and 23, which examines activity-induced turbulence in (7.18)–(7.22) by increasing $\zeta$ in this active model H; positive values of $\zeta$ are used for contractile swimmers, whereas negative values of $\zeta$ are appropriate for extensile swimmers. Padhan et al. (Reference Padhan, Vincenzi and Pandit2024a ) concentrate on $\zeta \lt 0$ , which yields activity-induced turbulence that suppresses phase separation. It has been suggested in Padhan et al. (Reference Padhan, Vincenzi and Pandit2024a ) that this model, with $\zeta \lt 0$ , might be applicable to a dense suspension of Chlamydomonas reinhardtii.

In figures 22(a) and 22(c) we show greyscale plots of the active-scalar field $\phi$ , at three representative times, which increase from left to right, for the activity parameter $|\zeta | = 0.1$ (see (7.18)–(7.22)); figures 22(d) and 22(f) contain pseudocolour plots of the vorticity $\omega$ corresponding, respectively, to the pseudocolour plots $\phi$ in figures 22(a) and 22(c). These plots indicate that, as time increases, the activity induces spatiotemporal chaos and turbulence; eventually the systems reaches a non-equilibrium statistically steady state in which coarsening is arrested by active turbulence (much as it is arrested by conventional fluid turbulence as we have discussed in § 6.3). We can characterise the statistical properties of this turbulence using the spectra and lengths that we have defined in (6.6) for phase separation in the binary-fluid case. In figure 23(a) we present log–log plots of the energy spectra $E(k)$ versus the wavenumber $k$ , for the activities $|\zeta | = 0.1$ and $|\zeta | = 1.5$ in (7.18)–(7.22). Clearly, the energy is spread out over a large range of $k$ as it is in fluid turbulence; furthermore, the power-law regimes in these spectra are consistent with $E(k) \sim k^{-5/3}$ (indicated by the dashed line). The root-mean-squared velocity $u_{rms}$ grows with $|\zeta |$ (figure 23 b). In figures 23(c) and 23(d) we present log–linear plots versus $|\zeta |$ of the mean coarsening length scale $L_c \equiv \langle \mathcal L(t)\rangle _t$ and the integral length scale $L_I$ , respectively. Figures 23(b) and 23(d) quantify the enhancement of turbulence in (7.18)–(7.22) with increasing $|\zeta |$ ; and figure 23(c) characterises coarsening arrest by this form of active turbulence.

7.9. Activity-induced droplet propulsion

Figure 24.

Activity-induced droplet propulsion: pseudocolour plots of $\psi$ (the magenta contour shows $\phi = 0$ ), at various times $t$ and activities $A$ ( $t$ increases from the top row to the bottom row): (a)–(d) $A = 0$ (complete phase separation in $\psi$ and no droplet propulsion); (e)–(h) $A = 0.15$ (rectilinear droplet propulsion); and (m)–(p) $A = 1$ (chaotic droplet propulsion). (i)–(l) For $A = 0.15$ : vector plots of $\boldsymbol {u}$ , with the overlaid $\phi = 0$ contour line (magenta) and the pseudocolour plot of $\omega$ , normalised by its maximal value (velocity vectors have lengths proportional to $|\boldsymbol {u}|$ ).

Figure 25.

(a) The multifractal time series for the scaled perimeter-deformation parameter $\Gamma (t)$ for various values of the activity. (b) The multifractal spectrum for representative activity $A = 1.5$ . We present the spectrum for a monofractal time series (given in blue) to show the robustness of the multifractal spectrum. The inset shows the plot of generalised exponent $\tau (q)$ as a function of the order $q$ for the representative value $A=1.5$ ; the deviation from the linearity suggests the multifractality of $\Gamma (t)$ .

We now use (7.23)–(7.29) (see § 7.7) to demonstrate activity-induced droplet propulsion in a model that has been studied in detail by Padhan & Pandit (Reference Padhan and Pandit2023a ). This model employs two scalar fields $\phi$ and $\psi$ ; both affect the velocity field $\boldsymbol {u}$ by which they are advected; but only $\psi$ is active in the parlance of active matter. Negative and positive values of $\phi$ and $\psi$ lead, respectively, to low and high densities of these scalars. We begin with the following initial data: a circular droplet, of radius $R_0$ and centre at $(x_{0,1}, x_{0,2}) = (\pi, \pi )$ :

(7.31) \begin{eqnarray} \boldsymbol {u}(\boldsymbol {x},t=0)&=&0;\nonumber \\ \phi (\boldsymbol {x}, t=0) &=& \tanh {\left (\frac {R_0 - \sqrt {(x_1-x_{0,1})^2 + (x_2-x_{0,2})^2}}{\epsilon _1}\right )};\nonumber \\ \psi (\boldsymbol {x}, t=0) &=& \begin{cases} \psi _0(\boldsymbol {x}) &\text {for}\ |\boldsymbol {x}|\leqslant R_0 ;\\ -1 &\text {for}\ |\boldsymbol {x}| \gt R_0 ; \end{cases} \end{eqnarray}

where $\psi _0 (\boldsymbol {x})$ is distributed uniformly and randomly on the interval $[-0.1,0.1]$ . We then monitor the spatiotemporal evolution of $\phi$ , $\psi$ and the normalised $\omega$ , which we depict in figure 24 via pseudocolour plots with an overlaid $\phi =0$ contour; the pseudocolour plots of $\omega$ also have superimposed vector plots of the velocity field $\boldsymbol {u}$ . The non-dimensional Weber numbers $\textit{We}_1 = R_0 U_0^2 /\sigma _1$ and $\textit{We}_2 = R_0 U_0^2/\sigma _2$ , Cahn numbers $\textit{Cn}_1 = \epsilon _1/R_0$ and $\textit{Cn}_2 = \epsilon _2 / R_0$ , Peclet numbers $\textit{Pe}_1 = R_0 U_0 \epsilon _1/(M_1 \sigma _1)$ and $\textit{Pe}_2 = R_0 U_0 \epsilon _2/(M_2 \sigma _2)$ , Reynolds number $\textit{Re} = R_0U_0/\nu$ , where $U_0 = {\left \lt U_{CM}(t) \right \gt _t}$ , with $U_{CM}$ , the speed of the droplet’s centre of mass (subscript $CM$ ), the order-parameter couplings $\beta _1^{\prime } = \beta \epsilon _1/\sigma _1$ and $\beta _2^{\prime } = \beta \epsilon _2 / \sigma _2$ , and the friction $\alpha ^{\prime } = \alpha R_0 / U_0$ , all affect the detailed dynamics of this initial droplet. However, most important of all these control parameters is the activity $A$ (7.30).

In figure 24 we exhibit activity-induced droplet propulsion in this model by illustrative pseudocolour plots of $\psi$ , with the $\phi = 0$ contour shown in magenta; we show such plots at different representative times, which increase from top to bottom, and three values of $A$ , which move from low to high values, as we move from left to right. In figure 24(a–d) we show the spatiotemporal of this droplet for the case of vanishing activity $A = 0$ ; there is no droplet propulsion and the system proceeds towards complete phase separation of $\psi$ , inside the $\phi = 0$ contour, by the formation of alternating annuli of regions with $\psi \gt 0$ and $\psi \lt 0$ , which is reminiscent of the phase separation of oil and water in a microfluidic droplet (see Moerman et al. Reference Moerman, Hohenberg, Vanden-Eijnden and Brujic2018). Figure 24(e–h) show that, when $A = 0.15$ , the system displays rectilinear propulsion of an active droplet; this is driven by the formation of an oscillating dipole that is visible clearly in figure 24(i–l), where we show, for $A = 0.15$ , vector plots of the velocity field $\boldsymbol {u}$ , with the $\phi = 0$ contour line (magenta), overlaid on a pseudocolour plot of the vorticity $\omega$ normalised by its maximal value. Finally, we show in figure 24(m–p), where $A = 1$ , chaotic droplet propulsion, which is characterised by significant fluctuations inside the droplet and on its boundary; the former suppress phase separation within the droplet and lead to diffusive or superdiffusive meandering of the centre of mass of the droplet (see Padhan & Pandit (Reference Padhan and Pandit2023a ) for details). The fluctuations of the boundary can be quantified by using the scaled droplet perimeter $\Gamma (t)$ (see (7.11) and figure 14 in § 7.3). In figure 25(a) we show that $\Gamma (t)$ has multifractal time series for various values of the activity $A$ ; we characterise this in figure 25(b) by showing the multifractal spectrum for the representative value $A = 1.5$ ; for comparison we present this spectrum for a monofractal time series in blue; the inset shows a plot of the generalised exponent $\tau (q)$ as a function of the order $q$ ; the deviation of the red curve from linearity quantifies the multifractality of $\Gamma (t)$ .