3.1. Linear-stability analysis

To determine the parameter values for which a uniform film with $h=h_i$ and $\psi =\psi _i$ is stable, where $h_i$ and $\psi _i$ are constants, we perform a linear-stability analysis. To do this, we substitute the following ansatz:

(3.1a)$$\begin{gather} h=h_i+\kappa \exp({{\rm i}(k_xx+k_yy)+\omega t}), \end{gather}$$

(3.1b)$$\begin{gather}\psi=\psi_i+\chi \exp({{\rm i}(k_xx+k_yy)+\omega t}), \end{gather}$$

into the coupled pair of (2.20) and then linearise in the amplitudes $\kappa$ and $\chi$, which are assumed to be small. This corresponds to a sinusoidal perturbation with dimensionless wavevector ${\boldsymbol k}\equiv (k_x,k_y)$ and dimensionless dispersion relation (growth rate) $\omega$. The aim is to obtain the dependence of $\omega$ on the wavenumber $k=|{\boldsymbol k}|$ and in particular the sign of $\omega (k)$. If $\omega (k)>0$ for any $k$, then the uniform film is unstable, since this indicates the perturbations in (3.1) grow over time. Taking together (2.20) and (3.1), Taylor expanding and then linearising in $\kappa$ and $\chi$, we obtain the system

(3.2)\begin{equation} \omega \begin{bmatrix} \kappa\\ \chi \end{bmatrix} = \begin{bmatrix} C_1 & C_2\\ C_3 & C_4 \end{bmatrix} \begin{bmatrix} \kappa\\ \chi \end{bmatrix} { = \bar {\boldsymbol {C}}\begin{bmatrix} \kappa\\ \chi \end{bmatrix} }, \end{equation}

where the matrix $\bar {\boldsymbol {C}}$ can be written as the product of two matrices, $\bar {\boldsymbol {C}}=\bar {\boldsymbol {Q}}\bar {\boldsymbol {J}}$, where the first is the linearised mobility matrix (cf. (2.2))

(3.3)\begin{equation} \bar{\boldsymbol{Q}} =h_i^3 \begin{bmatrix} 1 & \phi_i\\ \phi_i & \phi_i^2 + \dfrac{a^2\phi_i}{2{\rm \pi} h_i^2} \end{bmatrix}, \end{equation}

and the second is the following (largely thermodynamic in origin) matrix:

(3.4)\begin{equation} \bar{\boldsymbol{J}} ={-}k^2 \begin{bmatrix} k^2 + g''(h_i) + \phi_i \bar{F} & -\bar{F}\\ -\bar{F} & \dfrac{1}{\phi_i}\bar{F} \end{bmatrix}, \end{equation}

where $\bar {F} = (\phi_i \,f''(\phi_i)+ \epsilon '\phi _i k^2)/h_i=(K' - \alpha '\phi _i + 3\beta '\phi _i^3 + \epsilon '\phi _i k^2)/h_i$, and where we use $\partial {h}/\partial t = \omega \kappa \exp ({{\rm i}(k_xx+k_yy)+\omega t})$ and $\boldsymbol {\nabla }{h} = {\rm i}{\boldsymbol k}\kappa \exp ({{\rm i}(k_xx+k_yy)+\omega t})$, together with corresponding derivatives of $\psi$. The Jacobian matrix coefficients for our specific $g$ and $f$ are

(3.5)\begin{equation} C_1 ={-}h_i^{3}k^{2}\left[k^{2} + A'\left(-\frac{6}{h_i^{4}}+\frac{12}{h_i^{5}}\right) \right], \end{equation}

(3.6)\begin{equation} C_2 = 0, \end{equation}

(3.7)\begin{align} C_3 ={-}h_i^{3}k^{2}\phi_i\left[k^{2} + A'\left(-\frac{6}{h_i^{4}}+\frac{12}{h_i^{5}}\right) \right] + \frac{a^{2} \phi_i}{2 {\rm \pi}} k^{2}\left[K'-\alpha'\phi_i + 3\beta'\phi_i^3 + \epsilon'\phi_i k^{2}\right], \end{align}

and

(3.8)\begin{equation} C_4 ={-}\frac{a^{2}}{2 {\rm \pi}} k^{2}\left[K'-\alpha'\phi_i + 3\beta'\phi_i^3 + \epsilon'\phi_i k^{2}\right]. \end{equation}

Note that $C_3 = \phi _i (C_1 - C_4)$. We also observe that $C_2 = 0$, which originates from the multiple cancellations of terms that occur when the functional derivates of the free energy are substituted into the dynamical equations (see Appendix A), and it is an expression of the fact that bulk concentration gradients do not drive osmotic flow. In order to find the relationship between $\omega$ and $k$, we simply find the eigenvalues of the matrix

(3.9)\begin{equation} \begin{bmatrix} C_1-\omega & C_2\\ C_3 & C_4-\omega \end{bmatrix} . \end{equation}

The first eigenvalue is

(3.10)\begin{equation} \omega_h ={-}h_i^{3}k^{2}\left[k^{2} + A'\left(-\frac{6}{h_i^{4}}+\frac{12}{h_i^{5}}\right) \right], \end{equation}

and the second is

(3.11)\begin{equation} \omega_\psi ={-}\frac{a^{2}}{2 {\rm \pi}} k^{2}\left[K'-\alpha'\phi_i + 3\beta'\phi_i^3 + \epsilon'\phi_i k^{2}\right]. \end{equation}

It is because $C_2=0$ and $C_3 = \phi _i(\omega _h - \omega _\psi )$ that we obtain the simple results $\omega _h=C_1$ and $\omega _\psi =C_4$.

By writing the Jacobian matrix as the product of two symmetric matrices (shown in (3.3)–(3.4)) and noting also that $\bar {\boldsymbol {Q}}$ is positive definite allows one to write the linear problem as a generalised eigenvalue problem and therefore to directly deduce that (i) all eigenvalues are real, and that (ii) the stability thresholds do not depend at all on the mobility coefficients in $\bar {\boldsymbol {Q}}$, which is what one must expect, given the overall gradient dynamics structure in (2.1).

Note also that $\omega _h$ is the dispersion relation that is obtained for a pure liquid film, in the limit $\psi \to 0$. Moreover, $\omega _\psi$ is the dispersion relation that is obtained when considering just the colloids evolving in a liquid film that remains constant in thickness for all time. In other words, there is no coupling between film-height fluctuations and colloidal concentration fluctuations at the linear level. It is only at the nonlinear level that they are coupled. If we were to replace the binding potential in (2.4) with one including a coupling between $\psi$ and $h$, then our linear decoupling would no longer occur, as in the models of Náraigh & Thiffeault (Reference Náraigh and Thiffeault2010) and Thiele et al. (Reference Thiele, Todorova and Lopez2013). Such a coupling may even push an otherwise stable systems over the instability threshold.

Because of the decoupling just mentioned, the standard pure-liquid-film condition for $\omega _h(k)>0$ also applies here, i.e. one must have

(3.12)\begin{equation} \left.\frac{\partial^2g}{\partial h^2}\right|_{h_i}=A'\left(-\frac{6}{h_i^{4}}+\frac{12}{h_i^{5}}\right)>0. \end{equation}

The limit of linear stability for the liquid film is thus the locus in the phase diagram of the equation

(3.13)\begin{equation} \left.\frac{\partial^2g}{\partial h^2}\right|_{h_i}=0. \end{equation}

From (3.12) we see that this linear-stability threshold is simply $h_i=2$, i.e. for $h_i>2$ the film is linearly unstable. Note, however, that, when $h_i$ becomes large, the growth rate becomes very small and the unstable wavenumbers also become very small, i.e. corresponding to large length-scale instabilities.

Similarly, for the colloids to be stable, with $\omega _\psi (k)\leq 0$ for all $k$, we must have

(3.14)\begin{equation} \left.\frac{\partial^2f}{\partial \phi^2}\right|_{h_i,\psi_i}=K'-\alpha'\phi_i + 3\beta'\phi_i^3>0. \end{equation}

The corresponding limit of linear stability is just the locus in the phase diagram of the equation

(3.15)\begin{equation} \left.\frac{\partial^2 f}{\partial \phi^2}\right|_{h_i,\psi_i}=0. \end{equation}

Owing to the fact that we derive our dynamics from the free energy functional, as we show in the following subsection, this condition for the spinodal (onset of linear instability) for colloidal phase separation is precisely the same as that we obtain from the standard thermodynamic definition of the spinodal.

3.2. Bulk phase diagram for the colloids

When the colloids agglomerate, the system separates into two phases, one with colloid concentration $\phi _a$ and the other with concentration $\phi _b$. These two phases are in thermodynamic coexistence. The values of $\phi _a$ and $\phi _b$ vary as a function of the temperature and the other parameters in the model (i.e. $K'$, $\alpha '$ and $\beta '$). A plot of $\phi _a$ and $\phi _b$ in the temperature vs concentration plane yields the binodal curve. In the same diagram, we can also plot the spinodal curve from (3.15). These two curves meet at a single point in the phase diagram, namely the critical point (Hansen & McDonald Reference Hansen and McDonald2013). The binodal curve separates the single-phase and two-phase-coexistence regions of the phase diagram, while the spinodal curve is the linear-stability boundary: i.e. the system is linearly unstable if we choose the parameters so that the system is initially uniform and inside the spinodal curve. The phase diagram for our model is shown in figure 2, and here we discuss how it is constructed.

Figure 2.

Bulk colloid phase diagram in the plane of dimensionless temperature $K'/\alpha '=k_BT/(\alpha a^3)$ vs colloid concentration $\phi$, for the three values of $\beta '/\alpha '=\beta /\alpha = 0.8$, 1 and 2. The solid lines are the binodals and the dashed lines are the corresponding spinodals. The circles identify coexisting colloid concentrations for the particular temperature $K'/\alpha '=0.15$ and $\beta '/\alpha '=1$, that are referred to in § 4.1.

The binodal curve (coexisting states) can be determined by considering the thermodynamic requirements for two different phases to coexist: namely, they must be in thermal, mechanical and chemical equilibrium. This means that for phases $a$ and $b$ the temperature, pressure and chemical potential must be simultaneously equal

(3.16a–c)\begin{equation} T_a = T_b, \quad p_a=p_b, \quad \mu_a=\mu_b. \end{equation}

The first condition is equivalent to the requirement that the value of $K'$ be the same in both phases – see (2.18b). Taking the second and third conditions, using the expressions in (2.8) and (2.9) for $p$ and $\mu$ and remembering that the density of the colloids $\rho \propto \phi$, we obtain the following pair of simultaneous equations to be solved for the two coexisting concentrations $\phi _a$ and $\phi _b$:

(3.17)\begin{equation} \left.\frac{\partial f}{\partial \phi}\right|_{\phi = \phi_a}=\left.\frac{\partial f}{\partial \phi}\right|_{\phi = \phi_b}, \end{equation}

and

(3.18)\begin{equation} \left.\left[\phi\left(\frac{\partial f}{\partial \phi}\right)-f\right]\right|_{\phi = \phi_a} = \left.\left[\phi\left(\frac{\partial f}{\partial \phi}\right)-f\right]\right|_{\phi = \phi_b}. \end{equation}

Solving these for fixed $\alpha '$ and $\beta '$ and for a range of values of $K'$ yields the binodal curve. The corresponding spinodal is obtained from (3.15). These are shown in figure 2 for three different values of $\beta '/\alpha '=\beta /\alpha$. We plot the binodal and spinodal curves in the dimensionless temperature $K'/\alpha '=k_BT/(\alpha a^3)$ vs concentration $\phi$ plane. One can see that $K'/\alpha '$ and $\beta '/\alpha '$ are the relevant ratios to consider, by simply taking the spinodal condition obtained from (3.15) (see also (3.14)) and dividing through by $\alpha '$ to obtain the following expression for the spinodal:

(3.19)\begin{equation} \frac{K'}{\alpha'}=\phi_i - 3\frac{\beta'}{\alpha'}\phi_i^3. \end{equation}

From our dynamical equations, we identify the spinodal as a linear-stability threshold. However, from the thermodynamic point of view, it is the line in the phase diagram at which the isothermal compressibility $\chi _T$ diverges. This compressibility can be evaluated as (Hansen & McDonald Reference Hansen and McDonald2013)

(3.20)\begin{equation} \chi_T=\frac{1}{\phi}\frac{\partial\phi}{\partial p}. \end{equation}

From this expression, together with (2.8) and (2.9), one can easily show that that the spinodal condition in (3.15) is precisely the line in the phase diagram where $\chi _T\to \infty$.

Inspecting figure 2, we see that increasing the ratio $\beta '/\alpha '$ moves the critical point to the left and downwards in the phase diagram. Recall that $\beta '$ is the parameter originating from the steric repulsions between the colloids and $\alpha '$ is related to the strength of the attraction between the colloids. Thus, adjusting the ratio $\beta '/\alpha '$ for fixed $K'/\alpha '$ varies the value of $\phi$ in the dense colloidal phase, while at the same time also shifting the critical point. As long as $\beta '/\alpha '\sim {{O}}(1)$, the precise value makes very little difference qualitatively to the behaviour of our model. Thus, henceforth, we fix $\beta '/\alpha '=1$ (so that our phase diagram roughly matches some of the typical colloid phase diagrams displayed in chapter 12 of Hansen & McDonald Reference Hansen and McDonald2013, and references therein) and vary just the ratio $K'/\alpha '=k_BT/(\alpha a^3)$, which is equivalent to either varying the temperature or the strength of the attraction between the colloids. On the other hand, if modelling a particular experimental system, then adjusting the values of $\beta '/\alpha '$ and $K'/\alpha '$ to values different to those used here may be more appropriate.

3.3. Dispersion relations

Having obtained the dispersion relations (growth rates) in (3.10) and (3.11), we now discuss the different possible stability regimes. As mentioned already, the liquid film is unstable for $h_i>2$ and for all $A'>0$, where we have $\omega _h(k)>0$ for a band of wavenumbers $0< k<\sqrt {2}k_{h}$, where

(3.21)\begin{equation} k_{h}=\frac{\sqrt{3A'h_i(h_i - 2)}}{h_i^3} \end{equation}

is the fastest-growing mode, i.e. $\omega _h(k_h)$ is the maximum growth rate. The largest possible value of $k_h$ is $k_h^{max}=2\sqrt {60A'}/125\approx 0.124\sqrt {A'}$, which occurs when $h_i=2.5$.

Similarly, the colloids are unstable when the concentration takes a value inside the spinodal, where $f''(\phi )<0$ (see (3.14) and figure 2) and also where $\omega _\psi (k)>0$ for a band of wavenumbers $0< k<\sqrt {2}k_\psi$, where

(3.22)\begin{equation} k_\psi=\sqrt{\frac{-(K'-\alpha'\phi_i+3\beta'\phi_i^3)}{2\epsilon'\phi_i}}, \end{equation}

which is the fastest-growing mode, i.e. $\omega _\psi (k_\psi )$ is the maximum growth rate. The largest possible value of $k_\psi =k_\psi ^{max}$ occurs when $\phi _i=(K'/6\beta ')^{{1}/{3}}$, which corresponds to a curve in the phase diagram $\propto \phi _i^3$ that passes through both the origin and the critical point and where

(3.23)\begin{equation} k_\psi^{max}=\sqrt{\frac{{\alpha'-\left(\left(\dfrac{9}{2}\right)^{2}K'^2\beta'\right)^{1/3}}}{2\epsilon'}}. \end{equation}

Thus, we see that in general there are four possible stability scenarios, with four corresponding arrangements for the dispersion relations $\omega _h$ and $\omega _\psi$: (a) both the liquid and the colloids are stable (see figure 3a); (b) the liquid film is unstable, but the colloids are stable (see figure 3b); (c) the liquid stable, but the colloids are unstable (see figure 3c); and finally, (d) both are unstable (see figure 3d). In all cases we set $\beta '/\alpha ' = 1$, corresponding to the middle of the three phase diagrams in figure 2.

Figure 3.

Dispersion relation for four cases with $A' = 1$, $K' = 0.15$, $\epsilon ' = 0.4$ and $a^{2} = 100$. In case (a), both the film height and colloids are stable ($h_i = 1.9$ and $\phi _i = 0.15$). In case (b), the film height is unstable and colloids are stable ($h_i = 2.2$ and $\phi _i = 0.15$). In (c) the film height is stable and colloids are unstable ($h_i = 1.9$ and $\phi _i = 0.17$). In (d), both the film height and colloids are unstable ($h_i = 2.5$ and $\phi _i = 0.17$).

When both the liquid film and the colloids are unstable, then there are four possible arrangements of the dispersion relations, related to the question of whether $k_\psi >k_h$ or not and which lead to the fastest growth rate, i.e. whether $\omega _\psi (k_\psi )>\omega _h(k_h)$ or not. Recall that the fastest-growing modes $k_h$ and $k_\psi$ each correspond to initial fastest-growing wavelengths $\lambda _h = 2{\rm \pi} /k_h$ and $\lambda _\psi = 2{\rm \pi} /k_\psi$, respectively. Thus, suppose $\omega _\psi (k_\psi )<\omega _h(k_h)$, so that the fastest-growing mode is that corresponding to perturbations of the film height. When $k_\psi >k_h$, then the wavelength of the fastest growing film-height mode is longer than the corresponding fastest mode in the colloidal demixing. In contrast, when $k_\psi < k_h$, then the wavelength of the fastest film-height mode is shorter than the fastest colloidal demixing mode. However, when $\omega _\psi (k_\psi )>\omega _h(k_h)$, then the fastest-growing mode is that due to colloidal demixing and this can have a wavelength either longer or shorter than the fastest-growing film-height mode, depending on whether $k_\psi >k_h$ or not. All four of these cases are possible, depending on the choice of parameters in (3.10) and (3.11). Note that it is the presence of the two interfacial tensions $\gamma \propto 1/A'$ and $\epsilon '$ in the denominators of (3.21) and (3.22), respectively, that are the main players in determining the characteristic length scales $\lambda _h$ and $\lambda _\psi$.

4.1. Case when the film height instability is dominant and $k_h< k_\psi$

As a first case, we consider a situation where $\omega _h(k_h) > \omega _\psi (k_\psi )$, i.e. the film height instability grows faster and dominates over any colloidal modes, and also where $k_h< k_\psi$, so it is at longer wavelength. The relevant dispersion relations are shown in figure 4(a). Using the results in (4.4) to numerically evaluate the dispersion relations we can verify our numerical scheme through comparison with the analytical results in (3.10) and (3.11), which show excellent agreement.

Figure 4.

Results for a case with $A' = 2$, $K' = 0.15$, $\alpha '=1$, $\epsilon ' = 0.5$, $a^{2} = 2$, $h_i = 2.2$ and $\phi _i = 0.4$. (a) Dispersion relation calculated numerically via (4.4) using $\varepsilon _h = 10^{-7}$, $\varepsilon _\psi = 10^{-5}$ (symbols) compared with the analytic results in (3.10) and (3.11) (lines). (b) Measure of the evolution (cf. (4.5)).

In figure 4, we set the system length $L_x = 200$, with model parameters selected as $A' = 2$, $K' = 0.15$, $\alpha '=1$, $\epsilon ' = 0.5$ and $a^{2} = 2$, with $h_i = 2.2$ and $\phi _i = 0.4$ and, as previously mentioned, we fix $\beta '/\alpha '=1$ throughout. In figure 4(a) we see that the computational results (magenta dots for $h$ and black circles for $\psi$) lie on the analytical dispersion relation curve, with little error. Errors increase with system length $L_x$, since more growing modes fit in the domain (recall that all possible Fourier modes must satisfy our periodic boundary conditions), which then couple/interfere at an earlier stage in the dynamics. Through numerical experimentation we regularly see excellent agreement between the numerical and analytical values for the film-height eigenvalue $\omega _h$, whereas for the colloids there can be deviations of the numerical results from the analytically obtained curve $\omega _\psi$; they do, however, always follow the correct trend. This comparison makes for a very sensitive test of the numerics.

Next, we consider the evolution over time of the fastest-growing mode through our result in (4.5). In figure 4(b) we plot the logarithm of the maximum difference (as a function of position) between the initial and evolving heights. For this case, the film-height instability dominates, which is demonstrated by the slope of $0.0055$, matching the value of $\omega _h(k=k_h)$ from figure 4(a) with negligible error (cf. the slower maximum growth rate of the colloids is $\omega _\psi (k=k_\psi )=0.0013$). Additionally, we can predict that the system exits the linear regime when $\varepsilon _h \textrm {e}^{\omega _h t} \approx 1$, i.e. around $t \approx 3000$, which corresponds well with the region where linearity starts to be lost in figure 4(b). Note that the first few points of the curves in figure 4(b) do not correspond to the fastest-growing mode. This is because for simplicity we have chosen to quantify the evolution by considering the maximum change from the initial condition. The randomness of the initial condition (4.1) dictates that this does not initially coincide with a location demonstrating the growth of the fastest-growing mode. Thus, different random seeds exhibit minor differences, until the fastest-growing mode dominates over the effect of the choice of randomness. We further note that our result for $\psi$ in figure 4(b) shows that, after $t\approx 1500$, it does not grow with rate $\omega _\psi (k_\psi )$, due to being coupled to $h$; i.e. after $t\approx 1500$ the fastest-growing modes for both $h$ and $\psi$ are those with wavenumber $k_h$.

From the wavenumbers $k_h$ and $k_\psi$ of the fastest-growing modes in the system, the corresponding wavelengths $\lambda = 2{\rm \pi} /k$ and the length of the system $L_x = 200$, we can predict the number of peaks that initially form in the system. Using the values of $k_h$ and $k_\psi$ in (3.21) and (3.22), i.e. at the peaks in $\omega _h$ and $\omega _\psi$, together with the figure 4 parameter values, we have $\lambda _h \approx 42$ and $\lambda _\psi \approx 15$, so we expect at early times approximately 4 peaks in the film height, and 13 for the colloids. This is confirmed in figure 5, which displays the full time evolution of both $h$ and $\phi$.

Figure 5.

Panel (a) shows a waterfall plot of the local film height over time and (b) shows the corresponding local concentration of colloids. Both use a logarithmic scale in $t$. Panel (c) shows the final equilibrium profiles. These are for $A' = 2$, $K' = 0.15$, $\alpha ' = 1$, $\epsilon ' = 0.5$, $a^{2} = 2$, $h_i = 2.2$ and $\phi _i = 0.4$. In (c), the dashed black horizontal lines denote the two coexisting $\phi$ values, indicated in figure 2. Panel (d) shows the free energy of the system against time.

Having observed the agreement between our analytical and numerical results in the linear regime, it is interesting to explore the longer-time dynamics of our system. Figure 5 shows an example evolution, where we display the film height $h(x,t)$ in figure 5(a), the colloidal concentration $\phi (x,t)$ in figure 5(b), the final profiles of the system in figure 5(c) and the free energy of the system against time in figure 5(d). We can see that the free energy of the system always decreases over time. At early times we see the exponential growth or decay of modes as predicted by the linear theory. This corresponds to a decay of the small-wavelength modes, so the profiles initially appear to become smoother, but longer-wavelength unstable modes grow in amplitude, so that after some time the system starts forming peaks. As mentioned, the number of peaks is determined by the fastest-growing modes and can be predicted by considering the dispersion relations. In figure 5(a) these peaks become visible at $t\approx 2.5\times 10^{3}$ for the film height $h$, with the subdominant instability in the colloid profile taking longer to become visible, at around $t\approx 10^4$, as can be seen in figure 5(b).

The dynamics after the initial linear-growth regime can be complex, with evolution due to capillarity of the film height, spinodal decomposition and coarsening of the colloids, as well as dynamics driven by the coupling between the two. In the case shown in figure 5(a), multiple small droplets initially form and then subsequently coalesce, with translation occurring as a consequence of the coalescence events. The evolution is somewhat similar to that typically observed for films of a pure liquid evolving from a uniform film to form droplets and then ultimately a single droplet, as the liquid dewets from a surface. However, the major difference in the case here is that there are two sudden translations of the droplet around $t=2.3\times 10^7$ and $t=3.2 \times 10^{10}$, both of which arise through the coupling to the colloid concentration profile, which is exhibiting coarsening via Ostwald ripening – i.e. diffusion of colloids from one dense region to another through the low density regions (Lifshitz & Slyozov Reference Lifshitz and Slyozov1961; Wagner Reference Wagner1961) – to reduce the number of colloid-rich regions from 3 to 2, and then from 2 to 1.

Consider now the colloid evolution displayed in figure 5(b): the phase diagram in figure 2 shows that, for a given temperature ($K'/\alpha '$), there are two equilibrium values of $\phi$ that occur in a bulk system. These are indicated by the circles in figure 2, with concentrations $\phi _1 = 0.08$ and $\phi _2 = 0.6$. Their respective values being the intersection points between the horizontal temperature line and the corresponding binodal curve. These are the values that the colloidal concentration wants to evolve towards, with the interfacial tension term in our model (the term with coefficient $\epsilon$ in (2.3)) penalising interfaces and thus driving the system to have as few regions of dense/sparse colloids as possible. This is indeed observed, with regions of high density having a clear plateau that corresponds to the dense state point, and these regions evolving in a similar way to droplets through translation and coalescence or joining through Ostwald ripening. Whether aggregation occurs through translation and joining or through the Ostwald mode depends on the state point, size of the droplets and distance between them. The situations in which each coarsening mode dominates can be understood and determined using the methods presented in Pismen & Pomeau (Reference Pismen and Pomeau2004), Glasner & Witelski (Reference Glasner and Witelski2005), Glasner et al. (Reference Glasner, Otto, Rump and Slepčev2009), Dai (Reference Dai2010), Pototsky, Thiele & Archer (Reference Pototsky, Thiele and Archer2014) and Henkel, Snoeijer & Thiele (Reference Henkel, Snoeijer and Thiele2021).

For visual ease, and to highlight the relation to the phase diagram, we show the final profile of the system in figure 5(c). In this final state we observe that all colloids gather into a single region, with a single droplet too. Reaching equilibrium is slow for the film height, and even more so for the colloidal coarsening. To verify this profile as being final, we have separately considered the time evolution of the centre of mass and the free energy, both of which reach a plateau. The approach to equilibrium is faster when considering a 3-D situation (two spatial variables $x$ and $y$) since aggregation and droplet motion can occur in different directions, accelerating the process.

4.2. Case when the colloidal instability is dominant and $k_h< k_\psi$

A different case that is interesting to compare with the one in the previous subsection is a situation where both the film height and colloid profiles remain linearly unstable, with the fastest-growing mode for the film height being still at a longer wavelength than that of the colloids, but we swap the relative growth rates of the instabilities, i.e. $\omega _\psi (k=k_\psi ) > \omega _h(k=k_h)$, so that the colloid concentration fluctuations grow the fastest. We use the following parameter values to achieve this situation: $A' = 1$, $K' = 0.13$, $\alpha '=1$, $\epsilon ' = 0.5$, $a^{2} = 10$, $h_i = 2.5$ and $\phi _i = 0.4$. This situation is clear from the dispersion relations displayed in figure 6(a), where once again the analytical and numerical results show excellent agreement. The evolution of the profiles (not shown) is similar to that displayed in figure 5, except the dominant early-time linear-stage dynamics, which instead corresponds to the growth of colloidal concentration fluctuations. An interesting feature of the dynamics in this case can be observed from figure 6(b), where $L_x=300$. From (4.5) we expect a linear-growth regime at early times, which we do indeed observe, with the dominant colloid instability growing with a steeper slope than in the case shown in figure 4. However, the line in figure 6(b) for the film height has a sharp turn, departing at an early stage from its corresponding slower linear-growth regime. It is clear that the coupling to the colloidal demixing drives this later growth, which occurs at a rate greater than the rate predicted by linear analysis and therefore must be due to the nonlinear coupling terms. This difference from the case displayed in figure 4 can be traced back to the top-right zero in the matrix in (4.2), meaning that nonlinear couplings are required for fluctuations in $\psi$ to influence $h$.

Figure 6.

Results for a case with $A' = 1$, $K' = 0.13$, $\alpha '=1$, $\epsilon ' = 0.5$, $a^{2} = 10$, $h_i = 2.5$ and $\phi _i = 0.4$. (a) Dispersion relation calculated numerically via (4.4) using $\varepsilon _h = 10^{-7}$, $\varepsilon _\psi = 10^{-5}$ (symbols) compared with the analytic results in (3.10) and (3.11) (lines). (b) Measure of the evolution (cf. (4.5)).

4.3. Two strongly coupled cases where the long-wavelength colloid instability dominates

It is now of interest to consider cases where the coupling between the film height and the local colloid concentration is stronger. Inspecting the free energy in (2.14), we see that terms involving both $h$ and $\psi$ have coefficients $K'$, $\alpha '$, $\beta '$ and $\epsilon '$. To uniformly increase the coupling between the two fields without changing the colloidal bulk phase behaviour (rather than independently changing the relative strength of colloidal interactions, interfacial tension, etc.), we increase all of these parameters whilst maintaining the values of the ratios $K'/\alpha '$, $\beta '/\alpha '$ and $\epsilon '/\alpha '$ the same. The parameter values we now use are: $A' = 1$, $K' = 11$, $\alpha ' = \beta ' = 100$, $\epsilon ' = 4000$, $a^{2} = 100$ and $h_i = 2.5$. First, we consider the case $\phi _i = 0.4$. The dispersion relations (not displayed) show that the colloidal instability dominates at longer wavelength (smaller $k$) than the film-height instability. From the dispersion relations in (3.10) and (3.11) we obtain the fastest-growing wavelengths to be $\lambda _h = 52$ and $\lambda _c = 115$, so in a system of length $L_x=300$ we should observe around three maxima to initially develop in the colloid concentration profile and six in the film height. This simulation is displayed in figure 7(a–c), where we see that the early-time dynamics agrees with this calculation.

Figure 7.

Waterfall plots (using linear $t$) of (a) the film height, and (b) the local colloid concentration, for $A' = 1$, $K' = 11$, $\alpha ' =\beta ' = 100$, $\epsilon ' = 4000$, $a^{2} = 100$, $h_i = 2.5$ and $\phi _i = 0.4$. Panel (c) displays the final equilibrium profiles. Panels (d–f) show similar results, but with $\phi _i = 0.3$, corresponding to a decrease in the total amount of colloids in the system.

The simulation reaches an equilibrium state, displayed in figure 7(c), where the film height and colloid concentration profiles coarsen into a single drop. The equilibration process is much faster than in the cases considered in the previous subsections as there are fewer initial wavelengths, and the approach to equilibrium of the colloids is quicker because the case we consider here in figure 7(a–c) has all terms in the free energy that control the colloids multiplied by a factor of $100$. The influence of the much stronger coupling can be observed in figure 7(c), where either side of the droplet (peak in the film height) there are deep depressions in the film height that correspond in location to the interfaces between low and high density regions of the colloids. The influence of the strong coupling can also be observed during the evolution in the waterfall plot in figure 7(a), where these depressions appear as deep trenches. It is instructive to compare the equilibria of figure 7(c) with figure 5(c). In the latter, there is weak coupling, and hence the system evolves to a state where there is a droplet sitting on a precursor film at $h\approx 1.5$, corresponding to the ideal height based on the binding potential of our model, and the droplet being of a shape dictated by conservation of mass from that given by the initial condition and surface tension of the liquid. Similarly, the colloid concentration profile evolves to plateaus at high and low density phases with values given in the phase diagram in figure 2, with a minimum number of interfaces. By contrast, the highly coupled situation displayed in figure 7(c) has distorted profiles. In particular, we see strong depletion regions where the liquid height is significantly below the ideal precursor-film value, which occurs so as to accommodate the mass of colloidal particles in a single region within the droplet, and also the required interfaces between high and low density states. However, we would caution against ascribing too much significance to details of the ordering observed in the precursor film, since the most natural interpretation of films of this thickness is as a monolayer adsorbed on the surface, where the thin-film equation is arguably not a good description (Yin et al. Reference Yin, Sibley, Thiele and Archer2017).

As mentioned, in figure 7(c) we observe that the colloids end up within the single final droplet on the surface, i.e. the film height and the colloidal concentration profiles are in phase in the final result. This is not always the case; the outcome largely depends on the total amount of colloids within the system. In order to change the final result from being in phase to anti-phase (an example of the latter is shown in figure 7(f)) we can control the initial colloid concentration to a different location on the phase diagram in figure 2. In the case shown in figure 7(d–f), we decrease the initial average concentration to $\phi _i=0.3$, and achieve an anti-phase final result. We observe that the system is generally anti-phase when the total concentration of colloids is low (i.e. lower $\phi _i$, corresponding to locations on the left-hand side of the phase diagram within the spinodal in figure 2), and in phase for higher concentrations (i.e. higher $\phi _i$, corresponding to locations on the right-hand side within the spinodal). We investigate this matter further in § 5, with the aid of bifurcation diagrams.

4.4. Cases including stability of either film height or colloid profiles

As previously discussed, there are also situations where either a uniform film height or colloid concentration can be linearly stable. The case where both are stable is obvious, because in this situation any small perturbations in the initial profiles decay over time and the profiles both become flat. We also analysed some cases with profiles where one of the film height or colloid concentration profiles is stable, while the other unstable. Typical dispersion relations for cases like these are shown in figure 3(b,c). Examples of equilibrium profiles from our dynamic simulations are shown in figure 8.

Figure 8.

Final equilibrium profiles corresponding to cases where (a) the colloids are stable (and the film height is unstable) and (b) the film height is stable (with unstable colloids). The parameter values are the same as those for figures 3(b) and 3(c), respectively.

In figure 8(a), we have stable colloids and unstable film height. The film height evolves to form a droplet, while the concentration of the colloids remains roughly the same as the initial concentration $\phi _i$. As a result, the colloids follow the overall film height, as can be seen in the field $\psi$, so that the colloid profile ends up in phase with the film height profile. When the film height is stable, the system does not form a droplet. However, if the colloids demix, the film height is still affected by the unstable colloids. For example, in figure 8(b), we see that the colloids have separated into high and low density regions, and the interfaces between these distort the film height such that, together, the profiles minimise the overall free energy. We also notice that the computation time for this case is much longer. This is because when the film is stable, the colloids lose the driving from the film-height evolution, so the system reaches equilibrium in a much longer time. To put this another way: when the film height remains roughly constant, the only process governing the time evolution of the colloids is diffusive aggregation and this is a much slower process than the liquid dewetting dynamics.

4.5. Dynamics leading to asymmetrical final profiles

In our explorations of the model, we have also observed cases where the final profiles are asymmetrical. A typical example is displayed in figure 9, which occurs when $A' = 1$, $K' = 0.13$, $\alpha ' =\beta ' = 3$, $\epsilon ' = 0.5$, $a^{2} = 2$, $h_i = 2.5$ and $\phi _i = 0.41$. This corresponds to the dimensionless temperature $K'/\alpha '=0.043$, which is lower than in the cases considered previously. Here, the coexisting colloid concentrations are $\phi _1 = 0.00055$ and $\phi _2 = 0.77$ (cf. figure 2), so there is a greater difference between these than in the cases considered previously. The time evolution for a system of size $L_x = 60$ is displayed in waterfall plots in figure 9(a,b), with the corresponding final profiles in panel (c) and dispersion relations in (d). From the dispersion relations we see that the colloid mode is by far the fastest growing in the early stages. Thus, in the dynamics we see this instability dominating the overall dynamics and making the system exit the linear regime relatively quickly. We also see the colloid concentration variations coupling strongly to the film height $h$. Compared with the examples discussed previously, where we observe symmetric final profiles with the colloids distributed evenly outside the droplet (see e.g. in figure 5c), in the present case the coupling between $h$ and $\phi$ is much stronger. Therefore, the colloidal concentration variations are mirrored in the film-height variations over the surface. In figure 9(a,b) the initial number of peaks that form in $h$ and $\phi$ are accurately predicted by the dispersion relation. Since the system evolves to minimise the free energy, plots of the free energy always decreases over time (not displayed), with each decrease in the total number of colloid agglomerates (bumps) visible in figure 9(b) corresponding to a step-like drop in the total free energy. For some asymmetric cases we observe that, due to small numerical round-off errors, the final droplet can very slowly slide at constant velocity over the smooth surface. Such translations do not change the free energy. The sliding speed depends on choice of grid spacing and the direction can change with a different set of initial randomness. Of course, this cannot be a genuine feature, as confirmed with the aid of bifurcation diagrams in § 5.4.

Figure 9.

Panels (a,b) show waterfall plots over time (using $\log t$) of the film height and the local colloid concentration, respectively, for $A' = 1$, $K' = 0.13$, $\alpha ' =\beta ' = 3$, $\epsilon ' = 0.5$, $a^{2} = 2$, $h_i = 2.5$ and $\phi _i = 0.41$. Panel (c) shows the final equilibrium profiles and (d) shows the dispersion relation for this system.

5.1. Bifurcation diagram for a simple case with small coupling

In the previous sections, we showed examples of the dynamics of our system for a range of different parameters, and we noticed that the final film height and colloid concentration profiles equilibrate to being either in phase (peaks, or high density regions occur at similar positions) or anti-phase (where a peak in the film height corresponds to a lower density of colloids). We now investigate systematically the dependence of the equilibrium profiles on the system length $L_x$. We do this by producing bifurcation diagrams for some of the cases previously considered. The bifurcation diagrams are generated using our in-house numerical codes developed in Matlab employing a spectral method, see, for example, Lin et al. (Reference Lin, Tseluiko, Blyth and Kalliadasis2018) and Blyth et al. (Reference Blyth, Tseluiko, Lin and Kalliadasis2018); Blyth, Lin & Tseluiko (Reference Blyth, Lin and Tseluiko2023). Equations (2.20) are rewritten as a dynamical system for the Fourier coefficients of $h$ and $\psi$, which is truncated at a sufficiently high wavenumber to ensure accuracy. Steady-state solution profiles correspond to the fixed points of this dynamical system, resulting in a system of nonlinear equations for the non-zero-wavenumber Fourier coefficients of $h$ and $\psi$. Two additional equations are obtained from requiring fixed average values of $h$ and $\psi$. Additionally, we pin the solutions such that $h$ has a local extremum in the middle of the domain, i.e. we impose the condition $h_x|_{x=L_x/2}=0$. This necessitates the introduction of an additional unknown parameter, to maintain consistency between the number of imposed equations and unknowns. To achieve this, we introduce a fictitious wave speed, $c$, resulting in the addition of terms $ch_x$ and $c\psi _x$ to the right-hand sides of (2.20a) and (2.20b), respectively. This speed is trivially found to be zero in all continuation computations. Various solution branches of the resulting system of nonlinear equations for different parameter values are obtained by initially starting from small-amplitude nearly sinusoidal solutions on a domain size close to a cutoff wavelength obtained from the linear stability analysis. Numerical pseudo-arclength continuation is then performed, initially with respect to the domain-size parameter, and subsequently with respect to other relevant parameters, see, for example, Tseluiko, Baxter & Thiele (Reference Tseluiko, Baxter and Thiele2013), Lin et al. (Reference Lin, Rogers, Tseluiko and Thiele2016), Engelnkemper et al. (Reference Engelnkemper, Gurevich, Uecker, Wetzel and Thiele2019) and Tseluiko et al. (Reference Tseluiko, Alesemi, Lin and Thiele2020) for more details on numerical continuation techniques for liquid-film problems.

Our bifurcation diagrams depict the equilibria of our system when varying $L_x$ and keeping all other parameters fixed. Each realisation of the equilibrium profiles is represented by a point in each diagram (one for $h$ and one for $\psi$) corresponding to the $L^2$-norm of the profile, subtracting the profile average value $h_i$ or $\psi _i$, i.e.

(5.1)\begin{equation} \| h - h_i \| = \sqrt{\int_{0}^{L_x}| h - h_i | ^2 \, \textrm{d}x}, \end{equation}

and

(5.2)\begin{equation} \| \psi - \psi_i \| = \sqrt{\int_{0}^{L_x}| \psi - \psi_i | ^2 \, \textrm{d}x}. \end{equation}

In all the bifurcation diagrams shown here, blue lines represent continuation results starting from the known initial linear instability in the film height and red lines are similarly from the colloid instability (see § 3), termed here the film-height mode and the colloid mode. However, due to the coupling between the film height and the local colloid concentration, an instability in either affects both profiles, and hence both red and blue lines appear on the bifurcation diagrams for the norm of the equilibrium profile for both $h$ and $\psi$. Sections of the branches corresponding to stable and unstable solutions are shown with solid and dashed lines, respectively. Bifurcation points to side branches are illustrated using squares.

The bifurcation diagrams displayed do not show all possible branches, and hence not all possible equilibria of our system, but those included are instructive for understanding the dynamics of our system. Similarly, not all bifurcation points where the stability of solution branches changes are displayed, but some are illustrated with squares. In later cases in this section, some significant points where branches collide are discussed.

We begin by presenting the bifurcation diagrams in figure 10(a,b), which correspond to the dynamical simulations discussed in § 4.1. The parameter values are $A' = 2$, $K' = 0.15$, $\alpha '=1$, $\epsilon ' = 0.5$, $a^{2} = 2$, $h_i = 2.2$ and $\phi _i = 0.4$. Solution branches corresponding to inhomogeneous profiles begin on the $x$-axis at the points corresponding to where the system is first big enough to fit one unstable wavelength. We can read from figure 4(a) the largest wavenumbers when the system is unstable (the roots of $\omega _h(k)=0$ and $\omega _\psi (k)=0$) as $k_{h0}=0.215$ and $k_{\psi 0}=0.535$. Note that $k_{h0}\equiv \sqrt {2}k_h$ and $k_{\psi 0}\equiv \sqrt {2}k_\psi$, which can be obtained from (3.21) and (3.22), respectively. The corresponding wavelengths are $\lambda _{h0}=2{\rm \pi} /k_{h0}=29.2$ and $\lambda _{\psi 0}=2{\rm \pi} /k_{\psi 0}=11.7$, and are accurately captured in the bifurcation diagrams, figure 10(a,b), corresponding to the circular blue and red points on the $x$-axis, where the norms equal zero.

Figure 10.

Bifurcation diagrams and final states for cases shown in § 4.1, where $A' = 2$, $K' = 0.15$, $\alpha '=1$, $\epsilon ' = 0.5$, $a^{2} = 2$, $h_i = 2.2$ and $\phi _i = 0.4$. Panel (a) shows the film-height $L^2$-norm corresponding to two main branches of solutions for varying system size $L_x$. These originate from instabilities in either the film height (blue) or in the colloid local concentration (red), shown with circles. Squares represent locations of bifurcation points. Stable and unstable solutions are shown with solid and dashed lines, respectively. Similarly, (b) shows the $L^2$-norm of the corresponding $\psi$ profiles. Panels (c,d) show equilibrium profiles from continuation at the final $L_x=200$ point in (a,b): (c) shows the film-height branch (blue lines in a,b), and (d) the colloid instability branch (red dashed lines in a,b).

As $L_x$ is increased, the equilibrium profiles change. For $h$ in figure 10(a), we see that profiles on the two different branches originating from either the film-height or colloid instability are very similar, since both norm lines lie very close to each other. However, for $\psi$ in figure 10(b), the norms are very different and correspond to very different profiles. To illustrate this, we plot representative profiles for the points at $L_x=200$ in figure 10(c,d). In figure 10(c) the $\psi$ profile is largely flat with a small central bump, whereas in figure 10(d) the $\psi$ profile has a larger central bump and shoulders extending into the flat-film region, being significantly different from the flat solution, and thus has a much larger $L^2$-norm. This substantial difference is seen in figure 10(b), where the $L^2$-norm at $L_x=200$ on one branch is roughly twice the value on the other. Additionally, we see in the representative profiles in figure 10(c,d) that the maxima in both $h$ and $\psi$ occur together, which is the case all along these solution branches. Hence, for this set of parameter values, we conclude that profiles of any system length $L_x$ have the colloids located in phase with the centre of the droplet. We note that the solutions of the colloid-mode branch are linearly stable, whereas the solutions of the thin-film mode branch are all linearly unstable. In addition, there are a number of bifurcation points on the thin-film branch. Although we have not checked the nature of all these bifurcation points, our preliminary analysis shows that most likely they correspond to transcritical bifurcations. The side branches passing through these points are not shown here, as the solutions of these branches are linearly unstable. We note also that the equilibrium profiles in figure 10(d) agree well (up to horizontal translation) with the final profiles after the dynamic evolution of our hydrodynamic model in figure 5, giving confidence to both results.

5.2. Bifurcation analysis for a case with strong coupling

In § 4 we noted that the system can evolve to situations where the colloid local concentration profile can end up being either in- or anti-phase with the film height. We now investigate how this is manifested in the bifurcation diagrams showing the various possible equilibrium solutions, and how the transition occurs. First, we display in figure 11(a,b) bifurcation diagrams for the (stronger coupling) set of parameters used in figure 7(d–f), where the final profiles correspond to the film height and the colloid concentration profiles being anti-phase.

Figure 11.

(a,b) Bifurcation diagram for parameters as in figure 7(d–f), where now a second blue line branch corresponding to a double-droplet profile is depicted. The solid lines correspond to stable solutions, whereas the dashed lines correspond to unstable solutions. (c–e) Equilibrium profiles from the film-height mode (c) first film-height branch and (d) second film height branch; (e) equilibrium profiles from the colloid film mode.

In figure 11(a,b) the two blue lines are branches of solutions corresponding to instabilities in the film height. The first, as in figure 10, corresponds to one wavelength (single droplet) in the system and the second, starting at double the wavelength of the first branch, corresponds to profiles of two wavelengths (two droplets). As a result, the starting point of the second film-height branch is at twice the value (system length $L_x$) of the starting position of the leftmost branch. An example of a single wavelength solution on the first branch is shown in figure 11(c), and an example of a two wavelength solution is shown in figure 11(d) (recall that the domain is periodic in $x$). In the bifurcation diagram, we also display (red line) the continuation along the colloid-mode branch, that crosses the first film-height branch at the first bifurcation point, which turns out to correspond to a transcritical bifurcation (as can be more clearly seen in the inset, discussed in more detail later). The starting point of this branch is located to the left, at smaller $L_x$, of the second film mode branch and they do not intersect. This corresponds to an anti-phase profile, with a representative example shown in figure 11(e). This profile is the same as the one shown in our dynamical simulation in figure 7(f) with the droplet at a different location since in dynamical simulation we used a randomly perturbed initial condition.

The eigenvalues (not displayed) corresponding to solution branches displayed in figure 11(a,b) show that the part of the first blue film-height branch before the first bifurcation point is stable, and so is indicated as a solid line, while the rest of the branch is unstable, indicated as a dashed (blue) line. The second blue film-height branch is entirely unstable, so is all dashed.

The red colloid branch starts out unstable until reaching the leftmost turning point, except for the small part near where it crosses the film-height branch at the transcritical bifurcation point. This part is magnified in figure 11(a) in the top inset and the stable part is highlighted by the thicker red solid line. The branch becomes stable for a very small section between the turning point and the first bifurcation point on the film-height branch (corresponding to a transcritical bifurcation). The colloid branch is stable from the leftmost turning point until $L_x = 219.5$, where it becomes unstable again for a very short section (see the bottom inset).

When we increase the total average colloid concentration, the zero (neutral wavenumber) in the dispersion relation for the colloids shifts to smaller $k$, making the starting point for the branch of solutions corresponding to the colloid instability shift to larger $L_x$ (smaller $k$ corresponds to larger wavelengths). A sufficient increase can make the colloid branch start to the right of the second film-height branch. For such model parameters (here changing $\phi _i$ from 0.3 to 0.4 and keeping all other values the same), as we increase the system size $L_x$, we find that the colloid-mode branch now crosses the second film-height branch, as shown in the corresponding bifurcation diagrams displayed figure 12(a,b). In such a case, the colloid-mode branch crosses the second film-height branch at a pitchfork bifurcation point. We show the profiles at the intersection point in figure 12(f). Note that, at this point, the solution has the period equal to half of the domain size, but this symmetry is broken on the colloid branch away from this pitchfork bifurcation point.

Figure 12.

Bifurcation diagrams and profiles for an in-phase case, parameters as in figure 11 but with $\phi _i=0.4$. (a,b) Show bifurcation diagrams for the $L^2$-norm of the film height and colloidal profiles, respectively. Note that there are many other branches in the bifurcation diagram; those displayed correspond to the one- and two-wavelength solutions as predicted by our linear-stability analysis. The solid lines correspond to stable solutions, whereas the dashed lines correspond to unstable solutions. Panels (c–e) show profiles on the three displayed branches at $L_x=300$: (c) shows the profile on the first film-height mode branch, (d) on the second film-height mode branch and (e) on the colloids branch. (f) Shows the profile when the colloid-mode branch terminates on the second film-height branch.

To obtain the equivalent colloid-mode branch for lower $\phi _i$, we can take an already computed equilibrium solution for $\phi _i=0.3$ (for instance at $L_x=200$), continue the solution in $\phi _i$ to the desired value $\phi _i=0.4$ at the same $L_x$ and use this to then continue forward and backward in $L_x$, to generate the full colloid-mode solution branch displayed in figure 12(a,b) (red line). Doing this, on continuing backward in $L_x$ we find that this branch meets the second branch of the film-height bifurcation curve at a pitchfork bifurcation point, rather than the flat state (horizontal axis of our figures). We also use continuation to take the system forwards to $L_x = 300$. Again, the eigenvalue plot (not displayed) suggests that the first branch of the film-height mode is stable until the first bifurcation point in figure 12(a,b), and the entire second branch for the thin film is again unstable.

From figure 12(a,b) we can see that the main branch of the colloid mode (thicker red line) intersects the first branch of the film-height curve through the first bifurcation point (which turns out to correspond to a transcritical bifurcation), and lands on the first bifurcation point of the second film-height branch at a pitchfork bifurcation. We show the profiles from each branch at $L_x = 300$ in figure 12(c–e). The profiles from the film-height branches (shown in figure 12c,d) are similar to the anti-phase case because the parameters related to the film height have not changed. However, from the colloid branch (figure 12e) we see from the profiles that the film height and colloid concentration profiles are now in phase. We now explore this transition between anti-phase to in phase further.

5.3. Transition from anti-phase to in-phase solution profiles

With the two sets of bifurcation diagrams in figures 11 and 12, we understand that there is a fundamentally different structure between the anti-phase and in-phase cases. Critically, the transition occurs depending on whether the colloid mode begins at higher or lower $L_x$ than the second film-height branch, or, equivalently, the largest value of $k$ such that the colloids are linearly unstable compared with half the value for the film (or twice the wavelength). This change modifies the structure of the bifurcation diagram to enable an intersection with the second film-height branch or not. The transition thus occurs at this point. Starting from a high value and decreasing $\phi _i$, eventually the start point of the colloid branch collides with the start point of the second film-height branch and the connecting branch is annihilated below this critical value.

Another way to visualise the transition is to plot the maximum difference in film height against the initial concentration $\phi _i$; see figure 13(a) for $L_x=200$. We observe a clear turning point, located at $\phi _i = 0.367$ (the middle red point) and this is the critical concentration for the system to go from anti-phase to in phase. In the $L_x=200$ profiles displayed from the colloids branch, we can see a clear height difference between drops. In figure 13(b) (corresponding to the left red point in figure 13a), the droplet in the middle is taller, with the colloid concentration higher in the second smaller droplet (wrapping around the two sides due to periodicity), making the system anti-phase. On the other hand, in figure 13(d) (corresponding to the right red point in figure 13a), the drop with the most colloids in is taller, and hence in phase. In figure 13(c) we plot the profile of the system when the concentration is at the turning point of figure 13(a) (where $\phi _i = 0.367$), where the heights of both droplets are the same.

Figure 13.

Transition between anti-phase and in phase at $L_x = 200$. (a) Turning point in the maximum difference in film height. (b–d) Profiles for $L_x=200$ at $\phi _i=0.3$, $\phi _i=0.4$ and the critical concentration $\phi _i=0.367$, with other parameters as in figure 11.

5.4. Bifurcation diagrams for a case with asymmetrical solutions

In § 4.5 we described an example situation with asymmetrical final profiles. Here, we map out the bifurcation diagram for this case, displayed in figure 14. From noting that the neutral wavenumber in the dispersion relation in figure 9(d) for the colloids is much bigger than that for the film height, we can predict that the starting point for the colloid branch is at a much smaller wavelength (i.e. much smaller system size) than the one for the film-height branch. Arguably, the most interesting feature that we observe in figure 14 is the pitchfork bifurcation point on the first colloid branch, indicated by a red square. The side branch that bifurcates from this point is displayed as a red solid line. Interestingly, the solutions on this side branch are stable and asymmetrical. The profiles for $L_x = 60$ for both the symmetrical and asymmetrical colloid branches are shown in figures 15(a) and 15(b), respectively.

Figure 14.

Bifurcation diagrams for the same parameters as in figure 9. Panel (a) shows the $L^2$-norm for $h$, while (b) shows this for $\psi$. The red dashed line is the solution branch corresponding to unstable equilibria with symmetrical profiles. The red solid line is for the stable equilibria, which to the right of the red square are asymetrical. The blue dashed line corresponds to the unstable thin-film branch of solutions.

Figure 15.

Profiles for a system of length $L_x=60$ and corresponding to the bifurcation diagrams in figure 14. Panel (a) shows the unstable symmetrical solution on the first colloid branch and (b) shows the stable asymmetrical solution on the side branch bifurcating from the first colloid branch.

In figure 15(b), we display an asymmetrical profile, corresponding to the colloids being slightly gathered together at one end of the droplet, with not enough of them in the system to span equally to the other end of the droplet. Note that this equilibrium profile is almost the exact replicate of the final equilibrium obtained from our dynamical simulation, shown in figure 9(c). Note too, however, that, with different random seeds, we may observe equivalent profiles that are either translations or a mirror image of the equilibrium profiles shown in figure 9(c), i.e. in our dynamical simulations, where the centre of the droplet ends up and on which side of the droplet the colloids gather depends on the random initial conditions. To confirm the stability of the asymmetrical state calculated in the bifurcation diagram for $L_x=60$, we set it as an initial condition for our dynamical code, finding that over time it does not change, i.e. it is indeed stable.