Ghost-patterning and non-patterning in a draining film model

Patterns can form when the uniform state of any system is unstable so that some non-uniform motif grows in amplitude. Here, we identify an alternative way to form non-trivial structures, which we call “ghost-patterns”. Ghost-patterns emerge from noisy initial conditions when all non-uniform modes decay in amplitude except for one non-trivial motif which fails to decay. Hence, in seeking structured states, it is not necessary to find positive growth rates. We demonstrate ghost-patterns in an idealized non-equilibrium model intended to emulate draining thin-film suspensions.


Introduction
The emergence of complexity from simple systems has implications for the origins of life and the structure of the universe. Examples range from the structures formed by protein molecules, through the intricate patterns formed during embryonic morphogenesis and geological structures up to the clustering of galaxies.
Even in a far simpler experimental system -irregular quartz particles in a draining film of liquid-selforganized patterns have been reported [1] with a variety of morphologies. This system does not form large-scale inhomogeneities in the absence of flow and hence its patterning is an intrinsically non-equilibrium phenomenon, driven by the flux. Just how simple the building-blocks of a system can be, whilst yielding complexity and forming large-scale structures, remains an open question.
Motivated by this question, we have constructed an idealised model of a particle-laden fluid flowing over a structureless surface, as observed experimentally [1]. This is similar to Landau-Levich flow of a suspension of particles, which can form large-scale patterns (spinodal textures) if the particles are attractive [2]. But the existence of pattern-formation in attractive systems is well understood [3][4][5][6][7][8][9][10][11]; we are interested here in the regime in which capillary attractions are negligible (due to mean particle sizes an order of magnitude smaller than film thickness [1] 1 ) and the system is effectively infinite (far from any reservoir). The patterns thus formed [1] have a distinctly different morphology (cf. [2]). Patterning in this type of system occurs in many everyday fluids a e-mail: r.m.l.evans@leeds.ac.uk 1 Thus relaxing the assumption in ref. [1] that rare large particles are required to nucleate band formation.
such as sugar particles forming patterns around the top of a honey jar, fat particles in yoghurt or the dirt left behind on a window by rain water. Despite the range of systems involved there appear to be only two main types of pattern formed by draining particulate films [1]: one with channels parallel to the flow direction and the other with ridges perpendicular to it.
A fundamental explanation for these phenomena requires a model that is not specific to any one of the above mixtures, but omits extraneous system-specific features. Only by finding the simplest idealized model capable of reproducing the patterns can we discover the features that are essential to the pattern formation process. Thus, our idealized two-fluid model, described in the next section, incorporates only the simplest features of all draining thin-film suspensions: excluded volume, driven flow, and differential friction between the fluid components and substrate.
In a still simpler system -a single non-volatile liquid draining down a substrate-surface waves (variations in film thickness) have been observed [12] to form patterns with similar morphology to the particle distributions observed in the aqueous quartz suspensions [1]. But the similarity is likely to be superficial, as the film thickness remains fairly uniform during the patterning reported in [1]. The morphology of the film thickness at the leading edge of a particle-laden fluid, flowing down a substrate, has been studied in refs. [13][14][15], using models of the fluid dynamics that treat the suspension as a continuum, which we shall also do. While those model accurately reproduce the experimentally observed front or shocks in well-mixed, ergodic colloidal suspensions, they do not generate the patterning of the particle distribution observed throughout the films observed in [1]. We shall instead model the non-ergodic case of granular suspensions. And, in order to investigate mechanisms that are distinct from the well studied surface wave phenomena [16,17], we impose a condition of constant film thickness on our model.
In terms of modelling draining thin films, the investigation returned a null result; the model's features are not sufficient to generate the experimentally observed patterns. However, the model reveals a profound fact -that structure can emerge from initial randomness without the growth [18] of ordered modes, arising instead from the decay of all other modes. The emergent "ghost-patterns" are revealed when everything that is not the pattern decays.
In the next section, we introduce the idealised model system. Its linear stability is investigated analytically in sect. 3. This analysis reveals that the model system possesses two regimes of behaviour: one in which all fluctuations decay, leaving a uniform state, and one in which all fluctuations are neutrally stable, to linear order. Then, in sect. 4 and appendix B, numerical simulations of the model are presented, in which the neutral regime is analysed beyond linear order. In this regime, the numerics reveal that, while all spatial modes decay, the decay halts for the mode describing a structure of vertical channels, which then persist indefinitely.

The model
Our highly idealized model of a particle-laden fluid draining over a solid substrate describes the system as two interacting fluids, one representing the solvent and the other a course-grained representation of the particles and their interactions. In common with refs. [2,[13][14][15]19,20], we treat the set of particles as a continuum with a locally defined concentration, in order to investigate the most relevant features of the macroscopic dynamics. Detail at the level of individually resolved particles (as reviewed in [21]) is unnecessary in the present context, to discover the features essential to a universal process. This approach is valid in the limit where the length scales under consideration are much larger than the typical sizes of the particles (as is the case in ref. [1], where the thickness and width of the film are, respectively, one and three orders of magnitude larger than mean particle sizes).
We shall concentrate, in particular, on the case where the particles' behaviour is granular, rather than colloidal, so that particulate diffusion is negligible. Granular media have been successfully modelled as a continuous fluid in ref. [20] in the absence of a solvent. In the presence of a solvent, the motion of the granular "fluid" is strongly damped. Our model has no thermodynamic instability and therefore remains uniform in the absence of an external driving force.
As we are interested in the formation of large-scale patterning in the two macroscopic dimensions of the thin film, we construct a simplified two-dimensional model by coarse-graining over the third (z) dimension and neglecting variations in film thickness. The remaining dynamical variables in the two-dimensional space are the z-averaged solvent velocity v s and the z-averaged particulate velocity v p (which are both two-component vector fields), the local particulate volume fraction φ (a scalar field) and an effective pressure field p. All of the information about the dissipative interactions is then assumed to be accounted for in the effective viscosities and friction coefficients and all of the effects of excluded volume by an osmotic pressure.
Since the draining process is quite slow and the film thickness and particle size small, the flow is assumed Stokesian for simplicity 2 . So all stresses are linear in the true, three-dimensional velocities, and therefore also in the coarse-grained velocities v s and v p . It would be difficult to solve the full fluid-dynamical equations in the z direction, bounded by a static substrate and a free surface, in order to find the thickness-averaged drag forces on the solvent and particles. But, knowing that the solution must be linear in v s and v p , we can simply parametrize our two-dimensional model by three effective drag coefficients ζ s , ζ p and ζ sp , for relative motion between the two fluids and the fixed substrate or each other respectively, and by effective viscosities η s and η p . We shall assume that the film thickness remains approximately constant over the length and time scales considered, so that the effective drag coefficients and viscosities are constant.
Thus the momentum equation for the force-density on the solvent field (in terms of the Lagrangian time deriva- where F s is the body-force on the solvent per unit volume of solvent (so that (1 − φ)F s is the body-force on the solvent per unit volume of space), and the final equality is the Stokesian force balance condition. The total stress in the solvent is given by the usual Newtonian expression and body forces on the solvent arise from gravity g and drag, where the constant parameter ρ s is the mass density of pure solvent. The factor φ in the last term arises because the drag force from the particulate fluid must be proportional to the local amount of that fluid. Similarly, balancing force-density on the particle field gives where F p is the net body-force on the particles per unit volume of the particles and the constant ρ p is the mass density of the particles' material (e.g., quartz). Hence particles at concentration φ(x, t) have net mass per unit volume φρ p . The extra term ∇Π is the gradient of the osmotic pressure Π, which represents the direct excludedvolume interactions between particles. The osmotic pressure is a function Π(φ) of concentration, defining the equation of state of the particulate fluid. The non-osmotic contribution to its stress is assumed to be Newtonian for simplicity, We do not attempt to justify eq. (5) except by virtue of its simplicity. Recall that this is an effective relation in the idealized 2D system, which is notionally an approximate coarse-grained and thickness-averaged simplification of the 3D system. In reality, the effective inter-particle interactions are mediated by the fine details of the 3D solvent flow field. Recall also that our intention is to find the simplest physics capable of producing patterns, not to faithfully model the system in detail. A similar simplifying assumption was made for a simple granular medium in ref. [20]. The body forces on the particulate fluid are given by We see that the drag forces between solvent and particles (terms involving ζ sp in eqs. (3) and (6)) explicitly respect Newton's third law of motion where they enter the force density equations (1) and (4). Note that, while osmotic pressure Π is internal to the particulate fluid only, the pressure p acts on both fluids to respect incompressibility of the overall system, determined by continuity of their joint flux, By conservation of particles, evolution of the particulate concentration respects the continuity equation, Combining eqs. (1), (2), (3), (4), (5) and (6) gives and Equations (7), (8), (9) and (10), together with a prescription of the equation of state Π(φ), below, constitute the equations of motion for our model. These two vector and two scalar equations are sufficient to determine the two vector fields v s and v p and two scalar fields φ and p. The meanings of the symbols used here are summarized in table 1.
In order to create a simple model of the steric interactions between particles (in a similar vein to models in refs. [20,22]) the osmotic pressure is modelled by the continuous, piecewise doubly differentiable function, chosen for simplicity and efficiency of numerical calculations as shown in fig. 1. Here φ u and φ l are upper and lower thresholds of a plateau and φ max is the volume fraction at which the osmotic pressure diverges. The aim of this function is to provide a crude equation of state for a suspension in the granular regime, where diffusion is negligible, so that osmotic pressure is purely steric, and lacks any ideal-gas pressure. Hence the osmotic pressure does not increase until the volume fraction approaches random close packing, where it diverges. Thus the suspension has no tendency to spread out at lower concentrations (which would make patterns dissipate). Another divergence at φ = 0 is included to prevent the occurrence of negative concentrations. As in ref. [20], this simplified treatment of a granular medium neglects static friction of contacts.
Note that the factors of φ multiplying stress divergence, in eqs. (1) and (4), do not appear inside the divergence (as in [19]), because our model describes a granular suspension that is non-diffusive, unlike the models in ref. [19]. To see that this is correct, consider a case in which gravity is switched off, and a zero-velocity initial condition with a static arrangement of granular particles that is inhomogeneous but of sufficiently low concentration to be non-interacting (with Π = 0 everywhere). In the absence of diffusion, this inhomogeneous state should remain static forever, with the pressure field p uniform. Moving the factor of φ inside the divergence would produce a finite value of the spatial derivative on the RHS of those equations, and thereby wrongly generate finite time derivatives on the LHS, effectively describing an entropic (ideal-gas-like) contribution to the particulate osmotic pressure.
To reduce the number of free parameters, the equations of motion are non-dimensionalized in terms of the reduced quantitiest = ηp ηs ,ρ = ρp ρs and g = ρ s g ηs ζsp , wherev pd is the particulate fluid's drift velocity associated with the average of the initial volume To reduce numerical artefacts, the equations are also transformed into the co-moving frame of the initial drift velocityv pd , yieldinĝ

Linear stability analysis
To establish the linear stability of the system with respect to pattern formation, the equations of motion, eqs. (13), (14), (15), (16), are linearized about a uniform reference state, and the growth rate R of the Fourier modes of the linear perturbations in concentration 3 found to be where q is the magnitude of wavenumber of the mode in question and the subscript zero indicates that quantities are to be evaluated in the uniform state. A full derivation is given in appendix A.
Inspection of eq. (17) reveals that the uniform state is stable as long as the first derivative of the osmotic pressure at the average volume fraction is positive, since the other quantities in eq. (17) are positive. It is clear from fig. 1 that this is true at both low and high, average, volume fractions, so in these limits all forms of pattern are expected to decay to the uniform state. This is indeed observed to be the case in sect. 4, where the full non-linear equations of motion are numerically time-stepped. At intermediate volume fractions, however, the gradient of the osmotic pressure vanishes, so the linear stability analysis predicts neural stability and the results will be strongly influenced by non-linearities.

Numerical evolution of the non-linear equations of motion
Equations (13)-(16) were discretized on a rectangular grid so that the time evolution of the system could be calculated numerically. Details of the methods are given in appendix B. For various parameter values, the system was initialized in a near-uniform state of chosen mean volume fraction, with local small-amplitude random fluctuations of φ. Other initializations were also tested; see sect. 4.2.

Decay toward the uniform state
To check the validity of the calculations, the numerical algorithm was first tested for mean volume fractions where the equation of state (eq. (11)) has positive first drivative, Π (φ) > 0, corresponding to finite positive bulk modulus of the particulate fluid. In that regime, eq. (17) predicts that all fluctuations decay, and those on the longest length scales decay most slowly. To make quantitative comparison, the standard deviation of the volume fraction (a measure of the magnitude of the volume fraction fluctuations) was measured, and its decay with time is shown in fig. 2 for a selection of different parameter values.
Following the early-time decay of higher-frequency modes, the data demonstrate exponential decay of the slowest fluctuations, evident as a straight line on the logarithmic scale of fig. 2. Straight-line best fits to the linear section are shown as dashed lines in the figure. The measured rate of the exponential decay of the slowest mode is given in table 2 for the various parameter values. For comparison, theoretical values are given, calculated from eq. (17) for a wavevector q corresponding to the longest wavelength that can fit into the simulation cell. It is clear that the values obtained from linearization and from numerical time-stepping are in good quantitative agreement, with discrepancies below 2%. The small but statistically significant differences are almost certainly artefacts of the numerical discretization.
Note that the simulations used a non-zero value of the gravitational field g, so that the data in fig. 2 were obtained from continuously flowing systems. Nevertheless, the decay rate, in this regime, is predicted in eq. (17), as confirmed by the numerics, to be independent of both g and the ratio of particle to solvent densityρ, i.e. it is independent of the imposed flow field 4 .

Channel patterns
When the osmotic pressure at the mean volume fraction has zero gradient Π (φ), the linear stability analysis in sect. 3 predicts neutral stability of initially uniform states. Thus the behaviour of the system in this region is entirely controlled by non-linearities, requiring numerical investigation. The zero-gradient condition models the absence of diffusion in a granular medium, and therefore the absence of a driving force to smooth out inhomogeneities in particulate concentration.
The numerics reveal that the system evolves away from its initial state towards a state where vertical channels span the system. Figure 3 shows how these channels form from an initially random state. During this process the variance of volume fractions within the system initially falls, then stabilizes and remains constant once the channels span the whole system. Meanwhile, we find that no components of the structure factor (i.e., the spatial Fourier transform of the concentration φ) grow in amplitude; some components decay to zero amplitude while others decay to a finite asymptote. This is a "ghostpatterning" process, in the sense that the structure's amplitude is no larger than that of the initial randomness. The ghost-patterning is visible in the inserts of fig. 3 due to re-normalization of the grey-scale at each time.
The snapshots in fig. 3 show that obstructions, formed by slow-moving regions with high particulate volume frac-  fig. 2. The inserts show snapshots of the spatial distribution of the particle concentration φ at the times indicated on the plot. In these snapshots the volume fraction is normalised such that black corresponds to the highest volume fraction in the system at that time and white the lowest. (Using a constant normalization would yield a continuous reduction in contrast, obscuring the features.) The parameters used for this example were:ζ p = 2,ζs = 5,η = 7,ρ = 1.5, g = −15ŷ, φ 0 = 0.45, φmax = 0.64, φ l = 0.1 and φu = 0.54. tion, develop randomly in the system. These obstructions gradually lengthen due to build-up of particulate fluid behind them. Meanwhile the solvent that is able to pass through the obstruction drags some of the particles with it, elongating the region of higher volume fraction in front of the obstruction. In the absence of diffusion, the vertical features do not dissipate in the lateral direction.
It is worth considering how the values of the solvent and particle velocities vary throughout the system while the channels are developing. These are shown in fig. 4, along with the pressure and volume fraction, at two different stages of the process.
It is clear that the solvent travels significantly faster than the particles on average, due to its lower friction coefficient with the substrate. Note that, in the early stages ( fig. 4(a)), when the "obstructions" (high-φ regions) have just begun to spread in the vertical direction, the solvent flows around obstructions while the particle fluid reduces its velocity below the drift velocity (indicated by upwardpointing red vectors) as it encounters the osmotic pressure barrier of the obstruction's excluded volume. Meanwhile the particulate fluid travels faster in the regions where the volume fraction is lower, advected by the fast-moving solvent. The difference in the particulate velocity in the two regions (dark and light in fig. 4) is very clear once the channel pattern is well established, since there are long columns of cells where the particle velocity is significantly lower than the initial drift velocity.
Notice also the long-wavelength pressure variations. In the early stages high-pressure regions develop just behind the growing obstacles due to the retardation of solvent flow. This high pressure then drives the solvent outward around the obstacles where it meets least resistance. It is then drawn back by the low-pressure region on the trailing end of the obstacle.
Simulations performed with different initial states reveal that, while all initial configurations lead to a vertically channeled state, meaning that it is truly an absorbing state, the width, and spacing, of the channels are found to depend on the size of the obstructions in the system's initial state. Note that structure on the scale of the spatial discretization is present in fig. 3 only because it was present in the initial state. To demonstrate that the channels' stability is not an artefact of the spatial discretization, and to further quantify its evolution, we measure the evolution of its (two-dimensional) Fourier transform, starting from an initial state that is smoothly varying and contains just two Fourier modes, as depicted in fig. 5(a). In fig. 5(b)-(f), the results are shown of numerically time-stepping the initial state, using parameter values different from those in fig. 3 to demonstrate the robustness of the phenomenon. The concentration φ is shown in fig. 5(b) at times t = 1, 10, 100, 1000 where, after an initial movement of material, both vertically and horizontally, removing vertical variations and rarefying regions with concentration exceeding φ u (see fig. 1), the smoothly varying concentration field remained invariant during subsequent dynamics. Its standard deviation is shown in fig. 5(c).
The such that the modes initially present in eq. (18) are (m, n) = (2, 0) and (1,3). We see that the modes present in the ghost-pattern are not decoupled from the dynamics a priori, but initially evolve, until the system self-organizes into a state in which the pattern becomes decoupled. Notice that the (2, 0) mode ( fig. 5(d)), initially present (which has no variation in the y direction), initially decays, but eventually de-couples from the dynamics, and remains at finite amplitude at late times. Most modes that are not initially present (e.g., mode (3, 0) shown in fig. 5(d)) remain at negligible amplitude, but the (4, 0) mode grows due to it coupling to (2,0), and survives to late times due to its channel-like structure (without y-dependence).
Meanwhile, all modes with y-dependence (i.e., non-zero index n) decay at late times, as shown in fig. 5(e), including the initially present (1, 3) mode, and the (3, 3) mode which appears transiently as the flow distorts mode (1, 3). We have similarly Fourier-analyzed the horizontal component of the solvent velocity, which is small for all Fourier modes not present in the concentration distribu- tion, and decays monotonically with time for all modes except the (3, 3) mode, shown together with the (2, 0) mode in fig. 5(f).
We find that the ghost-pattern formation process and final morphology are robust with respect to variation of the various parameters of the model. For example, the effect of varying gravity is demonstrated in fig. 6, for a state with an initial state with weak random concentration fluctuations. Only the relative viscosityη affected the intensity of the fluctuations 5 , with higher viscosities resulting in fainter patterns. This is because, unlike the other forces, the viscous force opposes gradients in the flow field. Hence a high viscosity means that fast flowing, solvent-rich regions, cannot exist close to slow flowing particle-rich regions.
However, the reduced friction coefficients,ζ p andζ s , and relative viscosityη did have a profound effect on rate of ghost-pattern formation. Because this process is driven by the flow, the density and gravitational field strength also played a roll in determining this rate, with the strongest driving forces resulting in the system finding its steady state most quickly, as one might expect.

Concluding remarks
Our simple two-fluid model was created, to better understand the universal processes of pattering observed in draining thin films of particulate suspensions [1]. In the absence of osmotic forces (the granular regime), when ∂Π ∂φ = 0, the numerical solutions revealed an initial decay 5 Naturally the absolute value in the steady state was found to depend on the intensity of the fluctuations in the initial state.
of fluctuations towards a finite asymptote, corresponding to a faint structure of vertical channels, somewhat similar in shape to the channels observed in reference [1] in the non-granular regime, but lacking their intensity. The appearance of vertical structures in the model can be understood quite simply. The flow field smooths out any inhomogeneities in the volume fraction that occur in the flow direction and in doing so reduces the overall variance of concentration in the system, but does not continue to smooth out inhomogeneities in the direction perpendicular to the flow, as there is no net force in this direction to move the particles, once other (perpendicular) structure has disappeared. This demonstrates a novel behaviour that may exist more widely in dynamical systems: selforganization into a state in which particular non-trivial spatial modes are decoupled from the dynamics at late times, despite initially being coupled and time-dependent.
While the emergent structures bear some similarity to the vertical channels observed experimentally in draining thin-film suspensions, their amplitude is much lower, and the idealized model does not reproduce the horizontal bands seen experimentally. We infer that they require more than the simple ingredients of excluded volume and unequal substrate friction included in this model. Adapting the model [23] to include shear-thinning, shearthickening and concentration-dependence of the effective viscosity (all of which are phenomena observed in colloidal suspensions [24,25]) does not qualitatively alter the results [23]. We therefore speculate that the experimentally observed patterning requires a constitutive relation in which the particulate fluid's stress depends explicitly on the splay of its velocity field, as indicated by experimental observations of convergence-dependent jamming [26]. Alternatively, it might be necessary to include static contact friction between the particles and substrate, described by a finite yield stress, that might allow fast flowing, solventrich regions to exist close to slow flowing, particle-rich regions.
More significantly, the emergence of structuring from the decay of the non-structured initial fluctuations (rather than from an instability [18]) is an entirely novel class of process. This formation of morphology by ordered decay can usefully be called "ghost-patterning".
The new phenomenon of ghost-patterning is distinct from the phenomenon of critical slowing down, with which it shares some features. Critical slowing down occurs near a critical point of condensed matter [5], where restoring forces on inhomogeneities become small, so that fluctuations are long-lived, and the system decorrelates from its initial state slowly (as a power-law of distance from the critical point). In simple dynamical systems, similar slowing down occurs near a pitchfork bifurcation or on a saddle-node ghost [27]. In our ghost-patterning system, on the other hand, the dynamics remain fast whilst vertical inhomogeneities exist. They decay quickly and can also drive rearrangement and decay of the horizontal structure, which is non-trivially coupled to the fast modes via the flow dynamics and the osmotic pressure. Only once the system self-organizes into a channel pattern does it cease to evolve. Any structure without vertical variation and not exceeding the overlap concentration (the threshold for excluded-volume interactions) is a true steady state of the dynamics (not just a slowly evolving one). This can be Ghost-patterns, if small in amplitude, could be easily erased by the addition of diffusion to the dynamics. However, diffusion is negligibly slow in granular media, so systems and parameter regimes exist for which the erasure is effectively absent. Note also that the initial randomness is not necessarily small in amplitude, as it may not be thermal in origin. Initial random inhomogeneities in a distribution of large, athermal, granular particles could lead to robust, large-amplitude ghost-patterns.
CAH was supported by a Doctoral Training Allocation from the EPSRC.

Conflicts of interest
There are no conflicts of interest to declare.

Author contribution statement
RMLE conceived the work; the model was jointly designed and the paper jointly written; CAH preformed the linear stability calculation and the numerics and programmed the algorithm.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A. Linear stability analysis
The equations of motion, (13)- (16), expressing the x and y components of the velocity fields separately, and dropping the hat notationρ etc., become Making the substitutions p = p 0 + α, v p = v pd + β, v s = v sd + γ and φ = φ 0 + δ, where p 0 , v pd , γ and φ 0 are constants, and expanding to first order in the perturbations α, β, γ and δ gives  and using orthogonality of the Fourier modes and, by conservation of matter, ∂Dq ∂t = 0 when q = 0 gives (after rearranging and dropping the subscript q), Let D = de iθ , remembering that d and θ will depend on q. Then, after separating real and imaginary parts, we have (A.18b) Hence the magnitude of the concentration fluctuations will not grow, given that ∂π ∂φ | φ0 is non-negative, which is always the case in this model. The fluctuations decay exponentially with a characteristic rate (A.19) which is independent of gravity.