The effect of self-induced Marangoni flow on polar-nematic waves in active-matter systems

TL;DR

This paper explores how Marangoni flow affects density waves in active matter systems with polar-nematic symmetry.

Contribution

The study reveals how Marangoni flow influences the propagation and stability of polar-nematic density waves in active matter.

Findings

Density waves broaden and their speed changes with increasing Marangoni parameter.

Density waves can disappear via saddle-node or Hopf bifurcations.

Wave behavior depends on wavelength and mean density.

Abstract

We study the formation of propagating large-scale density waves of mixed polar-nematic symmetry in a colony of self-propelled agents that are bound to move along the planar surface of a thin viscous film. The agents act as an insoluble surfactant, i.e. the surface tension of the liquid depends on their density. Therefore, density gradients generate a Marangoni flow. We demonstrate that for active matter in the form of self-propelled surfactants with local (nematic) aligning interactions such a Marangoni flow nontrivially influences the propagation of the density waves. Upon gradually increasing the Marangoni parameter, which characterises the relative strength of the Marangoni flow as compared to the self-propulsion speed, the density waves broaden while their speed may either increase or decrease depending on wavelength and overall mean density. A further increase in the Marangoni parameter eventually results in the disappearance of the density waves. This may occur either discontinuously at finite wave amplitude via a saddle-node bifurcation or continuously with vanishing wave amplitude at a wave bifurcation, i.e. a finite-wavelength Hopf bifurcation.

Figures

Schematic representation of the system in a side view and b top view. Swimming agents represented by filled circles move along a flat liquid–air interface of a thin liquid film of uniform thickness *h*. On a coarse-grained level, swimmers of concentration *c*(*x*, *y*, *t*) act as a surfactant, inducing Marangoni flow with velocity $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}(x,y,z,t)$$\end{document}

Snapshots of fully developed polar-nematic travelling waves in 2D for different values of the Marangoni parameter illustrate the effect of Marangoni flow. Shown are a $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document}$ , b $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackag

Bifurcation diagram of the uniform steady states $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\rho _{1}^\textrm{s}, \rho _{2}^\textrm{s})$$\end{document}$ as a function of the mean density $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{docu

Two branches of travelling waves in a dry system ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document}$ ) at fixed $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document}$ and $\documentclass[

Three disconnected branches of finite-amplitude travelling waves for a dry system (Ma $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=0$$\end{document}$ ) at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document}$ a

Speed of travelling waves of different wavelength *L* as a function of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document}$ for $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document}$ . Star sy

Full text

Introduction

Meso-scale propagating density fronts and wavy patterns are known to spontaneously arise in systems of active matter, i.e. of self-propelled interacting agents (for comprehensive reviews, see [1–3]). These patterns emerge across a broad range of spatial scales, e.g. flocks of birds, fish, and insects on the metre-scale [4–9], auto-chemotactic suspension of microorganisms [10] and density waves in actin motility assays [11] on the micrometre scale. Driven permanently out of equilibrium by the agents’ motility, resulting self-organised complex spatio-temporal dynamic regimes may be indefinitely sustained by chemical fuel.

To model the characteristics of “dry” systems, where the motion of the fluid matrix supporting the self-propelled agents is not considered, the density distribution of the agents is described by the Boltzmann equation. It accounts for the microscopic rules governing binary collisions and random reorientations of the agents [12–14]. In the small-velocity regime, truncated forms of the Boltzmann equation can be derived by assuming dominant polar [12] or polar and nematic order parameter fields [13, 14]. Analysing the resulting set of closed hydrodynamic equations for these fields reveals that the observed spatio-temporal patterns in active matter originate in spontaneous symmetry breaking that occurs when the agent density surpasses a critical threshold [14].

The situation and resulting models are more complex for “wet” systems, defined by a two-way coupling of the motion of the self-propelled agents and the surrounding fluid medium. In this scenario, the Boltzmann equation for the agents’ density distribution is coupled to the Navier–Stokes equations that govern fluid motion. Here, we concentrate on a two-dimensional (2D) colony of self-propelled agents bound to move along the planar surface of a thin viscous film of liquid on a solid substrate. Then, the Boltzmann equation is extended by an advection term representing liquid flow at the film surface. Previously, we investigated the dynamics of a colony of non-interacting self-propelled agents on the surface of a liquid film with a deformable interface, particularly focusing on self-induced Marangoni flow [15, 16]. Such flows result from surface tension gradients that are due to concentration gradients of the agents.

When the average orientation of the self-propelled agents has a component orthogonal to the free surface, their motile force exerts an excess pressure similar to the hydrostatic pressure for a pending liquid film [15]. When the excess pressure dominates the stabilising surface tension, a Rayleigh–Taylor instability occurs [17, 18]. It is a long-wave stationary instability, i.e. its onset occurs at zero wavenumber and perturbations grow monotonically. In the classification of Ref. [19], this corresponds to a Cahn–Hilliard type. In a passive system (i.e. one that approaches thermodynamic equilibrium), such an instability results in the formation of drops that slowly coarsen. However, the situation changes if also a component of motile force exists that is parallel to the free surface. This results in the generation of intricate self-propelled dynamic density patterns and is studied in Ref. [16] by means of coupled thin-film and Smoluchowski equations. Note that in the context of the dynamics of an active Brownian particle, the Smoluchowski equation is essentially identical to the Boltzmann equation with an additional condition for the translational diffusivity constant which is linked to the temperature via the fluctuation–dissipation theorem [20].

In many practical scenarios, however, the orientation of the agents is always parallel to the liquid surface, which remains undeformed or weakly deformed at all times. An example of an active-matter system that meets these conditions is a two-dimensional bacterial bath, formed by stretching a dilute bacterial suspension on a wire mesh to create a free liquid film with a thickness of several micrometres [21–23]. In such cases, an isotropic disordered state may be unstable due to local interactions among the agents, driving the system towards spatio-temporal coherent dynamics and chaos.

Here, we consider a model system of swimmers moving parallel to the non-deformable free surface of a liquid film on a solid substrate, and study the dynamics of the relevant order parameter fields. In this way, we explore how the advective flow of the liquid matrix, induced by surface tension gradients affects the propagation of density waves. We expand a macroscopic continuum approach developed for dry active systems with local nematic alignment [13, 14]. In particular, we consider the wet active system that results from the incorporation of self-induced Marangoni flow into the hydrodynamic mean-field equations for the dominant polar and nematic order parameter fields.

Specifically, we demonstrate that upon gradually increasing the Marangoni parameter, which characterises the relative strength of the Marangoni flow as compared to the self-propulsion speed of the agents, density waves broaden, while their speed may either increase or decrease depending on wavelength and on overall mean density. This result echoes the influence of self-induced Marangoni flow on travelling reaction–diffusion waves confined in a liquid drop or film [24, 25]. There, the reaction is fuelled by an exothermic autocatalytic chemical reaction, such as the Belousov–Zhabotinsky or the iodate-arsenous acid reaction, and all reactants are non-motile. Similar to our observations, the inclusion of the Marangoni flow in autocatalytic chemical reaction systems leads to a change in propagation speed and front deformation. We also demonstrate two different bifurcation scenarios that result in the destruction of the waves at a relatively strong Marangoni flow. In the first scenario, a gradual increase in the Marangoni parameter results in the approach and ensuing collision of branches of stable und unstable finite-amplitude waves. The corresponding saddle-node bifurcation results in a discontinuous disappearance of the density waves. In the second scenario, a gradual increase in the Marangoni parameter results in a gradual decrease in the wave amplitude until it eventually vanishes in a continuous transition via a wave bifurcation, i.e. a finite-wavelength Hopf bifurcation, of a uniform state.

Governing equations

Consider a colony of self-propelled agents that are bound to move along a planar surface of a thin liquid carrier film of uniform and constant thickness h and dynamic viscosity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} that is supported by a solid bottom plate, as schematically shown in Fig. 1. The orientation of all agents remains parallel to the film surface at all times. Any two agents interact with each other via local binary collision rules. Thereby, any two colliding agents will align along their average pre-collision orientation when colliding at an acute angle and anti-align when colliding at an obtuse angle. Here, we adopt the kinetic Boltzmann approach developed for dry active matter in the absence of a fluid matrix [13, 14].Fig. 1. Schematic representation of the system in a side view and b top view. Swimming agents represented by filled circles move along a flat liquid–air interface of a thin liquid film of uniform thickness h. On a coarse-grained level, swimmers of concentration c(x, y, t) act as a surfactant, inducing Marangoni flow with velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}(x,y,z,t)$$\end{document} directed away from the region with higher density. The agents are advected by the surface flow field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}(x,y,t)=\textbf{u}(x,y,z=h,t)$$\end{document} , shown by blue arrows in (b). The orientations of individual agents are indicated by black arrows in b, while the coarse-grained particle density is colour-coded

To describe the interaction with the liquid, we neglect the size of the agents and assume that the main contribution to their motion is due to the advection by the liquid along the surface of the carrier film with the local velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}=(U_x,U_y)=\textbf{u}(x,y,z=h,t)$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}(x,y,z,t)$$\end{document} is the horizontal flow field in the bulk. We emphasise that the motion of solid particles of finite size trapped at an interface between two fluids is a complex hydrodynamic problem that (in the limit of small Reynolds number) involves the solution of the Stokes equation supplemented with no-slip boundary conditions at the particle surface and the continuity conditions for the flow field involving stresses at the interfaces as well as wetting conditions. For a comprehensive review of the motion of passive colloidal particles at liquid interfaces, see, for example, [26]. The complete hydrodynamic description becomes even more involved for self-propelled particles moving near penetrable liquid–liquid interfaces and can only be resolved numerically [27]. Here, we treat the particles as point-like implying that they are advected with the velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}$$\end{document} .

On the continuum level, the ensemble of self-propelled agents at the film surface is described by the one-particle probability density distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\textbf{r},\theta ,t)$$\end{document} that depends on the position in the film surface plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{r}=(x,y)$$\end{document} (that is parallel to the solid substrate) and the orientation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} of the agents. Its time evolution is described by the kinetic Boltzmann equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \partial _t \rho +\mathbf {\nabla }\cdot [(\textbf{U}+v_0 \textbf{e}(\theta ))\rho ]=d\mathbf {\nabla }^2 \rho + I_{\text {rot}}+ I_{\text {coll}}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0$$\end{document} is the self-propulsion speed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{e}(\theta )=(\cos (\theta ),\sin (\theta ))$$\end{document} is a unit orientation vector, d is the translational diffusivity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {rot}}$$\end{document} describes rotational diffusion, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {coll}}$$\end{document} is the collision term that describes the alignment of the agents. If the one-particle density distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\textbf{r},\theta ,t)$$\end{document} is normalised in such a way that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _A d\textbf{r}\int _0^{2\pi }d\theta \rho (\textbf{r},\theta ,t)=N$$\end{document} , where N is the total number of agents for a liquid film of surface area A, then their local number density per surface area is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c(\textbf{r},t)=\int _0^{2\pi }d\theta \,\rho (\textbf{r},\theta ,t)$$\end{document} . Then, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=A^{-1}\int _A c(\textbf{r},t)\,d\textbf{r}=NA^{-1}$$\end{document} is the mean number density.

The rotational diffusion term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {rot}}$$\end{document} in Eq. (1) depends on the one-particle density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\textbf{r},\theta ,t)$$\end{document} , while for binary collisions, the collision term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {coll}}$$\end{document} depends on the two-particle density function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2(\textbf{r}_1,\theta _1,\textbf{r}_2,\theta _2,t)$$\end{document} . However, in the mean field approximation one truncates the resulting hierarchy of equations assuming the closure relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2(\textbf{r}_1,\theta _1,\textbf{r}_2,\theta _2,t)=\rho (\textbf{r}_1,\theta _1,t)\rho (\textbf{r}_2,\theta _2,t)$$\end{document} [28]. This allows one to represent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {coll}}$$\end{document} as a quadratic functional of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} , as described in Ref. [13, 14].

Equation (1) must be supplemented with an equation for the advection flow field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}(\textbf{r},t)$$\end{document} . Similar to Ref. [15, 16] we assume that in the presence of the carrier liquid the self-propelled agents act as a surfactant, i.e. due to the solutal Marangoni effect, gradients in their density induce a Marangoni flow at the film surface. Note that translational diffusion and linear Marangoni effect both result from entropic contributions of an underlying energy functional as shown for insoluble surfactants in Ref. [29].

The flow field generated by a non-uniform concentration field c can be found by solving the Stokes equation in the fluid supplemented with the Marangoni stress boundary condition at the film surface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \partial _z \textbf{u}(x,y,z,t)\big \vert _{z=h}=- \frac{d\sigma }{dc} \mathbf { \nabla }c$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma (c)$$\end{document} is the concentration-dependent surface tension.

This problem has been solved in Ref. [30] by assuming that the deformation of the film surface can be neglected, which is valid at small capillary numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu U/\sigma _0\ll 1$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _0$$\end{document} is the surface tension of the bare liquid in the absence of surfactant. However, at large capillary number a complete description of the dynamics of a wet system would necessarily include the dynamics of film surface deformations. This will be addressed in a future study.

In particular, it was shown in Ref. [30] that for thick liquid layers of uniform thickness the surface velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}$$\end{document} induced by the Marangoni effect depends non-locally on gradients of the surfactant density c and can be represented as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{U}=-\frac{E h }{4\pi ^2 c_0 \mu }\int d^2k \frac{i\textbf{k} \hat{c}(\sinh ^2(|\textbf{k}|h)-(|\textbf{k}|h)^2)}{|\textbf{k}|h\sinh (2|\textbf{k}|h)-2(|\textbf{k}|h)^2}e^{i\textbf{k}\cdot \textbf{r}}. \end{aligned}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E=-c_0(d\sigma (c)/dc)\big \vert _{c=c_0}>0$$\end{document} is the Marangoni elasticity of the film surface that encodes the strength of change of the surface tension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma (c)$$\end{document} with c. Further, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{k}=(k_x,k_y)$$\end{document} is a 2D wave number and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{c}(\textbf{k},t)=4\pi ^2 \int d^2r\, c(\textbf{r},t)e^{-i\textbf{k}\cdot \textbf{r}}$$\end{document} is the Fourier transform of c.

In long-wave approximation [17, 31, 32], valid when the thickness of the liquid film is small as compared to the characteristic horizontal length scale of the variation in surfactant density, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h\ll 2\pi /|\textbf{k}|$$\end{document} , Eq. (2) reduces to the local relationship [17]

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{U}(\textbf{r},t)=-\frac{Eh}{4 c_0 \mu } \mathbf {\nabla } c(\textbf{r},t). \end{aligned}$$\end{document}

It is interesting to observe that the functional dependence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}$$\end{document} on c in Eq. (3) is identical to the one found for the velocity of a spherical solid colloidal particle transported at an interface by diffusiophoresis [26].

The uniform disordered state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\textbf{r},\theta )=c_0/(2\pi )$$\end{document} may become unstable to spatial and orientational perturbations, leading to a transition towards dynamic inhomogeneous ordered states including steady clusters or bands, and propagating waves as it was demonstrated previously for dry systems [12–14]. One of the key observations derived from the numerical analysis of the Boltzmann equation (1) is that despite the short-range interactions due to collisions, which typically occur when the distance between the agents is comparable with their size, the global dynamics of the system exhibits a long-range order with a much larger characteristic scale. This allows us to replace the Boltzmann equations for the one-particle density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (x,y,\theta ,t)$$\end{document} with a set of coupled hydrodynamic equations for certain coarse-grained order parameter fields that do not depend on the orientation angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} . As such, the Fourier transform of the density in the angular space is used.

The orientation order on a large scale is well described by the first few order parameter fields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho _q(\textbf{r},t)=\int _0^{2\pi } e^{iq\theta }\rho (\textbf{r},\theta ,t)\,d\theta . \end{aligned}$$\end{document}

In fact, for polar or polar-nematic alignment between the agents during collisions, the orientation order is well approximated by the first four dominant order parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2,3,4}$$\end{document} with all other order fields assumed to be negligibly small, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{q>4}=0$$\end{document} .

Our goal is to follow the closure scheme outlined in [13, 14] to derive a closed set of coarse-grained hydrodynamic equations for the density c and the dominant order parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}$$\end{document} for a wet system by including the advection velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}$$\end{document} . Thus, we also assume the following scaling and truncation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & |\textbf{U}(c)|\sim v_0,\quad \partial _{t,x,y}\sim \epsilon ,\quad c-c_0\sim \epsilon ,\nonumber \\ & \rho _{1,2}\sim \epsilon ,\quad \rho _{3,4}\sim \epsilon ^2,\quad \rho _{q>4}=0 \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \ll 1$$\end{document} is the dimensionless super-criticality, i.e. the distance from the onset of the order–disorder transition. From the physical point of view, the above scaling implies a slow dynamics with a long-wave spatial orientation order dominated by polar \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1$$\end{document} and nematic \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2$$\end{document} order parameter fields. The scaling (5) automatically restricts the applicability of the hydrodynamic approximation to relatively small magnitudes and slow variations of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{2}$$\end{document} . In consequence, similar to other weakly nonlinear approaches, the approximation is valid only in the vicinity of the order–disorder transition where the disordered state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c=c_0$$\end{document} becomes unstable.

The closed form of the dynamic equations for the density c and the complex-valued polar and nematic order parameter fields \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2$$\end{document} is obtained in several steps as detailed in Appendix A. First, the Boltzmann equation is Fourier-transformed in the angular space and truncated at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=4$$\end{document} , i.e. all order parameter fields with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q>4$$\end{document} are identically zero. Next, using scaling arguments (5), the fields \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{3,4}$$\end{document} are expressed in terms of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}$$\end{document} . Finally, after substituting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{3,4}$$\end{document} into the equations for c and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}$$\end{document} , the closed set of hydrodynamic equations is obtained as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \partial _t c + \mathbf {\nabla }\cdot (\textbf{U}c)=-\frac{v_0}{2}\left( \partial \rho _1^* + \partial ^* \rho _1 \right) +d\mathbf {\nabla }^2 c\nonumber \\ & \partial _t \rho _1+\mathbf {\nabla }\cdot (\textbf{U}\rho _1)=-(\alpha _0+c\alpha _1)\rho _1+\alpha _2\rho _1^*\rho _2 -\alpha _3 |\rho _2|^2\rho _1\nonumber \\ & -\frac{v_0}{2}(\partial c + \partial ^* \rho _2)+\gamma _1 \rho _2^*\partial \rho _2+d\mathbf {\nabla }^2 \rho _1 \nonumber \\ & \partial _t \rho _2+\mathbf {\nabla }\cdot (\textbf{U}\rho _2)=(-\beta _0+c\beta _1)\rho _2+\beta _2\rho _1^2-\beta _3 |\rho _2|^2\rho _2\nonumber \\ & -\beta _4|\rho _1|^2\rho _2-\frac{v_0}{2}\partial \rho _1 +\gamma _2 |\partial |^2 \rho _2 \nonumber \\ & -\gamma _3 \rho _1^* \partial \rho _2 -\gamma _4 \partial ^*(\rho _1\rho _2)+d\mathbf {\nabla }^2 \rho _2. \end{aligned}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial =\partial _x + i\partial _y$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf {\nabla }=(\partial _x,\partial _y)^T$$\end{document} , while the star indicates the complex conjugate. The system (6) is closed by Eq. (2) [or Eq. (3) in the long-wave limit] that relates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}$$\end{document} to the number density c. Note that we have retained translational diffusion with diffusivity d to ensure numerical stability and consistency with the Marangoni term. For ease of comparison, system (6) is non-dimensionalised using the scaling introduced in Ref. [14] for a dry system, i.e. when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}=0$$\end{document} . Namely, time is scaled with the characteristic duration of the ballistic flight path of the agents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau =r^{-1}$$\end{document} , where r is the rate of random orientation changes of the agents. The spatial coordinates are scaled with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0\tau $$\end{document} , and the one-particle density distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} is measured in units of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r/Dv_0$$\end{document} , where D is the size of the agent particles. In dimensionless units, the dynamics of a wet system in the long-wave limit is described by Eqs. (6) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0=1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}=-\text {Ma}\mathbf {\nabla }c$$\end{document} where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=E h/(4 v_0^2\mu \tau )$$\end{document} represents the Marangoni parameter that controls the relative strength of the Marangoni flow as compared to the self-propulsion speed.

The coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _i$$\end{document} in Eqs. (6) determine the properties of the spatially uniform states, while the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _i$$\end{document} are responsible for the coupling between the gradients of polar and nematic order parameter fields.

As outlined in Appendix A, the closure relations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{3,4}$$\end{document} and, therefore, the coefficients in Eqs. (6) depend on the average density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} , as well as on the microscopic parameters such as the statistics of the rotational diffusion, the effective interaction distance and the alignment rules associated with binary collisions. The behaviour of the system will therefore depend on changes in these microscopic parameters even if the average coarse-grained density is kept fixed. However, since our aim is to understand the effect of the self-induced Marangoni flow on the dynamics of the system on the coarse-grained level, we apply a semi-phenomenological approach as proposed in Ref. [14], i.e. we derive the governing equations by coarse graining but then pursue a parametric study on the continuum level. Accordingly, the obtained hydrodynamic Eqs. (6) are studied for selected physically meaningful control parameters, while all other parameters remain fixed. In particular, our main focus is the effect of the Marangoni number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} and of the mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} . Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} can be varied independently at fixed density for a given set of the microscopic parameters.

Following Ref. [14], we choose the coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _i$$\end{document} to correspond to a system that is known to feature several transitions in absence of Marangoni flow. When the mean number density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} is varied at fixed values of the coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _i$$\end{document} , transitions occur between a uniform disordered state and a uniform ordered nematic state and as well between the latter and a uniform ordered polar-nematic state of mixed symmetry. The parameter values employed in the following are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha _0= & 0.011,\quad \alpha _1=0.44,\quad \alpha _2=1.5,\quad \alpha _3=4.4,\nonumber \\ \beta _0= & 0.044,\quad \beta _1=0.59,\quad \beta _2=-0.16,\quad \nonumber \\ \beta _3= & 6.6,\quad \beta _4=-1.0,\nonumber \\ \gamma _1= & 1,\quad \gamma _2=0.85,\quad \gamma _3=0.23,\quad \gamma _4=3.7. \end{aligned}$$\end{document}

Coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _i$$\end{document} can be traced back to the microscopic parameters and the average density via the relations given in Appendix A. In fact, the values used here correspond to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.08$$\end{document} for the choice of the collision rules used in Ref. [14]. The translational diffusivity is varied between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=10^{-2}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=10^{-1}$$\end{document} .

Marangoni flow-induced change in polar-nematic waves

As shown in Ref. [14], the dry system, i.e. without the presence of Marangoni flow, develops polar-nematic travelling waves for mean densities above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=\beta _0/\beta _1=0.074$$\end{document} . Before embarking on a detailed analysis, we first demonstrate that self-induced Marangoni flow influences the propagation of these density waves. We use the long-wave variant [Eq. (3)] to calculate the Marangoni flow1 and numerically solve Eqs. (6) in 2D using central differences for spatial derivatives and a semi-implicit second-order Heun method for time integration. As initial condition, we use a uniform disordered state with small-amplitude random initial perturbations with a slight polar and nematic bias in the x-direction. Equations (6) are numerically integrated over a sufficiently long time interval of 5000 units at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} to allow the system to reach a steady state. Then, the average wave number associated with the Fourier-transformed density field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int |\textbf{k}| \hat{c}(\textbf{k},t)\,d\textbf{k}$$\end{document} fluctuates about a constant value and does not show any trend as a function of time. The density field c recorded over the subsequent time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t=1000$$\end{document} is used to extract the statistics of the speed of the waves (not shown). The dynamics of the system at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}>0$$\end{document} is obtained using naive parameter continuation. Namely, after the system is integrated for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t=6000$$\end{document} , the Marangoni parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} is abruptly incremented by 10 and Eqs. (6) are further integrated without setting a new initial condition.Fig. 2. Snapshots of fully developed polar-nematic travelling waves in 2D for different values of the Marangoni parameter illustrate the effect of Marangoni flow. Shown are a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} , b \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=10$$\end{document} , c \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=30$$\end{document} , and d \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=150$$\end{document} . Accompanying movies in the Supplementary Material further illustrate the wave dynamics. The mean density is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} , the diffusivity is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=0.1$$\end{document} , and a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$800\times 800$$\end{document} domain with periodic boundaries is used. The background colouring indicates the number density c as indicated in the colour bar. Black arrows indicate the vectorial polar order parameter field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\text {Re}(\rho _1),\text {Im}(\rho _1))$$\end{document} . No stable waves are found for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}>160$$\end{document} . All remaining parameters are given in the main text

Snapshots of fully developed travelling wave patterns are shown in Fig. 2 for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=0.1$$\end{document} in a square domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$800\times 800$$\end{document} with periodic boundaries. These waves appear to be long-time stable and propagate without persistent coarsening or splitting. Inspecting Fig. 2, one discerns that the inclusion of Marangoni flow results in a broadening of the polar-nematic wave in the direction of propagation as well as a moderate increase in its average propagation speed that can be estimated from the movies in the Supplementary Material (SM). The dependence of the propagation speed of the waves on the Marangoni parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} is addressed in detail in the next section.

An interesting feature of the travelling waves is their apparent stability with respect to coarsening, i.e. for any fixed value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} , the waves seemingly propagate without merging to form a broader wave of larger wavelength. Coarsening occurs only during the transient phase directly after an abrupt increase of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} when the system evolves into a new stationary state with a larger average wavelength. Thus, upon incrementally increasing the Marangoni parameter, the number of wave crests reduces from five in the dry system at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} to two at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=30$$\end{document} and finally only one at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=150$$\end{document} . Further increasing the Marangoni parameter beyond a certain critical value ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}\gtrapprox 160$$\end{document} ), the wave dissipates into a uniform nematic state.

Suppressed or completely arrested coarsening is a well-known phenomenon in active systems like matter driven by external forces, active matter (i.e. motile matter) and reaction–diffusion waves. Generally, there exist two physical mechanisms that underpin coarsening, namely (i) the volume mode (also known as Ostwald ripening or mass transfer) and (ii) the translation mode (also known as coalescence where two or more clusters merge into one) [31, 33, 34]. In active systems formed by non-motile agents, suppression and arrest of coarsening is observed in the presence of an external driving force studied, e.g. with convective Cahn–Hilliard models [35, 36], or is caused by non-equilibrium chemical potentials studied, e.g. with nonreciprocal Cahn–Hilliard models [37].

In systems of active matter, i.e. composed of motile agents, self-propulsion leads to the formation of clusters of agents that do not merge to form a single dense phase [38–40]. As discussed in [41], the exact physical mechanisms behind arrested coarsening in active matter systems are still unclear. However, on the phenomenological level, suppressed coarsening could be due to the non-variational nature of the evolution equations (6) which is inherent to all active matter systems. The question of the exact physical origin of the arrested coarsening requires a detailed stability analysis of various travelling wave states and will be addressed in future studies.

One-dimensional polar-nematic waves

To gain a deeper understanding of the influence of Marangoni flow on the characteristics of the density waves, we next analyse a one-dimensional (1D) approximation. Namely, we assume that waves like the ones shown in Fig. 2 are nearly translationally invariant in the direction orthogonal to the direction of propagation. In particular, we consider waves that travel in x-direction and are invariant in y-direction. This implies that all fields depend only on the x-coordinate, that the imaginary parts of the order parameter fields are zero \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Im}(\rho _1)=\text {Im}(\rho _2)=0$$\end{document} , i.e. Eqs. (6) with (3) become

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \partial _t c -\text {Ma} \,\partial _x ( c\partial _x c)=- \partial _x \rho _1+d\partial _x^2 c\nonumber \\ & \partial _t \rho _1-\text {Ma}\,\partial _x ( \rho _1 \partial _x c)=-(\alpha _0+c\alpha _1)\rho _1+\alpha _2\rho _1\rho _2\nonumber \\ & -\alpha _3 \rho _2^2\rho _1-\frac{1}{2}(\partial _x c + \partial _x \rho _2)+\gamma _1 \rho _2\partial _x \rho _2+d\partial _x^2 \rho _1\nonumber \\ & \partial _t \rho _2-\text {Ma}\,\partial _x (\rho _2\partial _x c)=(-\beta _0+c\beta _1)\rho _2+\beta _2\rho _1^2-\beta _3 \rho _2^3\nonumber \\ & -\beta _4\rho _1^2\rho _2-\frac{1}{2}\partial _x \rho _1 +(\gamma _2+d) \partial _x^2 \rho _2 -\gamma _3 \rho _1 \partial _x \rho _2 \nonumber \\ & -\gamma _4 \partial _x(\rho _1\rho _2). \end{aligned}$$\end{document}

Note that the system (8) is invariant under the transformation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(t,x,c,\rho _1,\rho _2)\rightarrow (t,-x,c,-\rho _1,\rho _2)$$\end{document} , i.e. it is of mixed parity. In the following, we analyse the 1D model by determining uniform steady states and their linear stability as well as travelling wave states.

Uniform states and their linear stability

First, we discus the uniform steady states of Eqs. (8) that exist at the present parameters [Eq. (7)] using the mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} as control parameter. The trivial state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1^\textrm{t}=\rho _2^\textrm{t}=0$$\end{document} exists for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} . Two other branches of uniform steady states exist in certain ranges of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} : The first branch bifurcates in a pitchfork bifurcation from the trivial state and exists for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0>\beta _0/\beta _1=0.074$$\end{document} . It corresponds to purely nematic uniform states

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho _1^\textrm{n}=0,\quad \rho _2^\textrm{n}=\pm \sqrt{(-\beta _0+c_0\beta _1)/\beta _3}. \end{aligned}$$\end{document}

The second branch starts and ends in respective pitchfork bifurcations on the branch of nematic states and exists for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} in the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.091\le c_0\le 0.2425$$\end{document} . It represents a mixed polar-nematic uniform state with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho _2^\textrm{pn}= & (2\alpha _3)^{-1}[\alpha _2\pm \sqrt{\alpha _2^2-4\alpha _3(\alpha _0+c_0\alpha _1)}].\nonumber \\ \rho _1^\textrm{pn}= & \pm \sqrt{(\beta _2-\beta _4 \rho _2^\textrm{pn})^{-1}[\beta _3(\rho _2^\textrm{pn})^3+(\beta _0-c_0\beta _1)\rho _2^\textrm{pn}]}, \end{aligned}$$\end{document}

Note that Eqs. (9), (10) are valid for dry systems as well as for wet systems, since a uniform density distribution does not induce a Marangoni flow. Figure 3 presents a bifurcation diagram featuring all uniform steady states. Both, polar order parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1$$\end{document} and nematic order parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2$$\end{document} , are given as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} . The three-dimensional representation together with the projections of the bifurcation diagram onto the ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0, \rho _1$$\end{document} )- and ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0, \rho _2$$\end{document} )-planes well represent the connections between the different states and the symmetries of the system. In the case of a dry system, the bifurcation diagram in Fig. 3 can be directly compared with a slice of the stability diagram obtained in Ref. [14] for varying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _2$$\end{document} . Indeed, at the here used value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _2=1.5$$\end{document} the critical density for the nematic-polar transition is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} as obtained in Ref. [14].

To determine the linear stability of the uniform states that is already indicated in Fig. 3, we use the harmonic perturbation ansatz

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} c= & c_0+\varepsilon \tilde{c}e^{\lambda t+ikx},\nonumber \\ \rho _1= & \rho _1^\textrm{s}+ \varepsilon \tilde{\rho }_1e^{\lambda t+ikx},\nonumber \\ \rho _2= & \rho _2^\textrm{s}+ \varepsilon \tilde{\rho }_2e^{\lambda t+ikx}, \end{aligned}$$\end{document}

where the superscript s stands for any of the above introduced uniform steady states t, n, pn, and the quantities with tilde are the scaled perturbation amplitudes. Linearising Eq. (8) in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \ll 1$$\end{document} about any of the uniform states one obtains an eigenvalue problem for the eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} and corresponding eigenvectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\rho }=(\tilde{c},\tilde{\rho }_1,\tilde{\rho }_2)^T$$\end{document} , namely

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda \varvec{\rho }= (\mathbf {\underline{J}}_0 + ik \mathbf {\underline{ J}}_1 + k^2 \mathbf {\underline{J}}_2) \varvec{\rho } = \mathbf {\underline{J}} \varvec{\rho }. \end{aligned}$$\end{document}

Note that the real and imaginary part of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} correspond to the growth rate and frequency of harmonic modes with wave number k, respectively. The propagation speed of travelling modes is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=\text {Im}(\lambda )/k$$\end{document} . The matrices within the Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf {\underline{J}}$$\end{document} are given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathbf {\underline{J}}_0= & \left( \begin{array}{ccc} 0 & 0 & 0\\ -\alpha _1\rho _1^\textrm{s} & -(\alpha _0+c_0\alpha _1)+\alpha _2\rho _2^\textrm{s}-\alpha _3(\rho _2^\textrm{s})^2 & \alpha _2\rho _1^\textrm{s}-2\alpha _3\rho _1^\textrm{s}\rho _2^\textrm{s} \\ \beta _1\rho _2^\textrm{s} & 2\beta _2\rho _1^\textrm{s}-2\beta _4\rho _2^\textrm{s}\rho _1^\textrm{s} & (-\beta _0+c_0\beta _1)-3\beta _3(\rho _2^\textrm{s})^2-\beta _4(\rho _1^\textrm{s})^2 \end{array} \right) ,\nonumber \\ \mathbf {\underline{J}}_1= & \left( \begin{array}{ccc} 0 & -1 & 0\\ -\frac{1}{2} & 0& -\frac{1}{2}+\gamma _1\rho _2^\textrm{s} \\ 0& -\frac{1}{2}-\gamma _4\rho _2^\textrm{s} & -(\gamma _4+\gamma _3)\rho _1^\textrm{s} \end{array} \right) ,\mathbf {\underline{J}}_2= \left( \begin{array}{ccc} -(d+\text {Ma}\,c_0) & 0 & 0\\ -\text {Ma}\,\rho _1^\textrm{s} & -d & 0 \\ -\text {Ma}\,\rho _2^\textrm{s} & 0 & -(d+\gamma _2) \end{array} \right) . \end{aligned}$$\end{document}

Fig. 3. Bifurcation diagram of the uniform steady states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\rho _{1}^\textrm{s}, \rho _{2}^\textrm{s})$$\end{document} as a function of the mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} in a 3D representation. Shown are the trivial state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^{\textrm{t}}=0$$\end{document} (red line) from which the purely nematic state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} [blue lines, Eq. (9)] emerge in a pitchfork bifurcation at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.074$$\end{document} . From the latter branch the mixed polar-nematic state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} [green lines, Eq. (10)] emerges in pitchfork bifurcations at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.2425$$\end{document} . Solid (dashed) lines indicate stability (instability) w.r.t. uniform perturbations. Considering spatially extended systems, the state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} first destabilises at a finite wave number in a wave instability (at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.086$$\end{document} , marked by a square symbol and a blue “W”). Note that for clarity the mirroring \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} -branch at negative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2$$\end{document} is not included. Parameters are given in Eq. (7)

Note that the Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf {\underline{J}}$$\end{document} is complex due to the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ik \mathbf {\underline{J}}_1$$\end{document} that is related to self-propulsion. In consequence, for any wave number k, up to three different eigenvalues may exist, each either complex or real. However, in the case of the trivial uniform state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^{\textrm{s}}=\rho _{1,2}^{\textrm{t}}=0$$\end{document} and in the case of the purely nematic state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^{\textrm{s}}=\rho _{1,2}^\textrm{n}$$\end{document} [Eq. (9)] the resulting characteristic polynomial \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det (\mathbf {\underline{J}}-\lambda \mathbb {I})$$\end{document} has only real coefficients. Therefore, all three eigenvalues are either real, or one eigenvalue is real and the remaining two form a complex conjugate pair. This feature is similar to isotropic systems where the Jacobian is always real. In contrast, for the mixed polar-nematic state with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^{\textrm{s}}=\rho _{1,2}^\textrm{pn}$$\end{document} the characteristic polynomial has complex coefficients resulting in three different, in general, complex eigenvalues.

The eigenvalue problem Eq. (12) is solved numerically for various mean densities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} and Marangoni parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} , to determine the regions of linear stability (where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Re}(\lambda )<0$$\end{document} ) for the three uniform steady states. Resulting stability diagrams for the nematic and the nematic-polar states in the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(c_0,\text {Ma})$$\end{document} are given in Fig. 4a, d, respectively, as explained in the following. As expected, the trivial state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1^\textrm{t}=\rho _2^\textrm{t}=0$$\end{document} is (independently of Ma) linearly stable at small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} and becomes unstable at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=\beta _0/\beta _1=0.074$$\end{document} . The instability is stationary and large-scale, i.e. the onset is at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0$$\end{document} . At the instability threshold, the above-mentioned pitchfork bifurcation occurs where the uniform purely nematic state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} [Eq. (9)] emerges. For slightly supercritical \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} , the leading eigenvalue of the trivial state is real and is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda (k)=-\beta _0+c_0\beta _1-(d+\gamma _2)k^2-k^2/[4(\alpha _0+c_0(\alpha _1+\beta _1)-\beta _0)]+O(k^4)$$\end{document} .Fig. 4a Linear stability of the purely nematic uniform state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} [Eq. (9)] in the plane spanned by mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} and Marangoni parameter Ma. Labels “s”, “w”, and “m” mark the linearly stable region, i.e. the region where all wave numbers are stable, the region of the wave instability (oscillatory instability of a finite band of wave numbers bound away from zero) and the region of mixed mode instability (unstable band of wave numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\le k<k_c$$\end{document} with real or complex eigenvalues depending on k), respectively, see main text. b, c give typical dispersion relations, i.e. growth rates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Re}(\lambda )$$\end{document} for the leading eigenvalues as a function of the wavenumber k, in region w at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} and in region m at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.095$$\end{document} , respectively, at two different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} as indicated near each curve. Dashed (solid) lines correspond to complex (real) eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} . d Linear stability of the mixed polar-nematic uniform state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} [Eq. (10)]. In the regions marked by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {I}$$\end{document} [ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {II}$$\end{document} ] at most one [two] eigenvalues have a positive real part for any given k. e and f give dispersion relations in region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {I}$$\end{document} at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document} and region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {II}$$\end{document} at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=30$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.18$$\end{document} , respectively

Following the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} -branch for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} when gradually increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} , we find a wave instability at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.086$$\end{document} , i.e. a small-scale oscillatory instability, sometimes also called a finite-wavelength Hopf or oscillatory Turing instability. It is marked by a blue square symbol and the letter “w” in Fig. 3. Beyond this point the uniform nematic state is always unstable for sufficiently large domain sizes. At the wave bifurcation, travelling nematic bands emerge as discussed in Ref. [14]. Figure 4b shows the dispersion relation, represented by the real parts of the two leading eigenvalues in the parameter range between the onset of the wave instability and the pitchfork bifurcation where the uniform \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} -state emerges from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} . Dashed (solid) lines correspond to complex (real) eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} . The dispersion relations are qualitatively similar in the entire region marked “w” in Fig. 4a. In this region, there exists a band of unstable wave numbers k; namely, the eigenvalue \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} is complex with a positive real part. The Marangoni flow suppresses the instability as shown in Fig. 4b.

At the onset of the wave instability (W-point in Fig. 3), the leading eigenvalue of the Jacobian is purely imaginary at some critical wave number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_c$$\end{document} , i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda (k_c)=i\Omega $$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=\Omega /k_c$$\end{document} is the speed of the propagating non-decaying wave that is at onset of infinitesimal amplitude. Further increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} beyond the W-point, the already unstable uniform nematic state at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} additionally undergoes the above described spatially uniform stationary instability where the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} -branch emerges. Beyond this point, the dispersion relation for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} -state features a real positive eigenvalue at small k that connects to a pair of complex conjugate eigenvalues at larger k, see Fig. 4c. Together they form a continuous band of unstable wave numbers. Note that in the vicinity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} the wavenumber band of the stationary instability becomes very narrow. Exactly at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} , it only contains the point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0$$\end{document} implying that a this single parameter value the dispersion relation corresponds to a conserved-Hopf instability [19, 42].

We call the combination of real and complex eigenvalues above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} “mixed mode” and mark the corresponding region in Fig. 4a by “m”. Here, large-scale perturbations are monotonically unstable, while perturbations of smaller scale are oscillatory unstable. However, both mode types form a common band of unstable wave numbers. Sometimes the transition point between real and complex modes within the band is called an “exceptional point”. Note, however, that such points only directly affect the primary bifurcations if they occur directly at zero crossing—then they correspond to a Takens–Bogdanov bifurcation. As they are detected in a linear analysis, their potential influence on the nonlinear behaviour needs additional analysis. Inspecting Fig. 4a, we see that an inclusion of Marangoni flow suppresses the wave instability of the nematic state. It completely disappears at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}\approx 100$$\end{document} . Interestingly, this further implies that in contrast to the case studied in [42], here the conserved-Hopf instability appears as a codimension-2 instability exactly where the line separating the s- and w-region in Fig. 4a ends.

It was shown already in Ref. [14] for the dry case that the time evolution of the system in the mixed mode region is dominated by nematic-polar wave patterns characterised by spatially distributed polar and nematic order parameter fields. Analysing the linear stability of the uniform polar-nematic states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} [Eq. (10)] for a spatially extended system, we find as expected that the state inherits the wave instability from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{n}$$\end{document} -branch it emerges from at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.091$$\end{document} . A typical dispersion relation at moderate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} is given in Fig. 4e and shows that the band of unstable wave numbers is bound away from zero. However, at lower and higher densities the unstable band can also start at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0$$\end{document} (not shown).

The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{1,2}^\textrm{pn}$$\end{document} -branch retains this instability till it ends at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.2425$$\end{document} . It even acquires a second instability mode in an intermediate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} -range, see the dispersion relation in Fig. 4f shown for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.18$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=30$$\end{document} close to the codimension-2 point in Fig. 4d, where the lines separating the regions “I” and “II” intersect. In these regions I and II, there exist one and two bands of unstable wave numbers, respectively.

From the stability diagram Fig. 4d, we see that an increase in the Marangoni parameter tends to suppress the wave instabilities: at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}\approx 100$$\end{document} , the II-region has been replaced by I-regions, i.e. the second wave mode has been suppressed, and a stable range of concentrations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} marked by “s” in Fig. 4d has started to develop around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.2$$\end{document} .

Note that at the codimension-2 point in Fig. 4d, where the two lines meet that separate regions I and II, the two wave instabilities have a simultaneous onset. Adjusting parameters that control the corresponding critical wave numbers, one might in the future be able to find resonances with interesting nonlinear behaviour.

Travelling waves in the dry system

The dispersion relations in Fig. 4e, f indicate that in regions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {I}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {II}$$\end{document} of Fig. 4d, respectively, one and two wave modes are unstable. As both modes correspond to bands of unstable wave numbers bound away from zero, the one-mode [two-mode] case corresponds to two [four] bifurcations where branches of travelling waves emerge when using the domain size as control parameter. In situations as in the first row of Fig. 4, there is one wave mode with two (Fig. 4b) or one (Fig. 4c) bifurcation(s). Depending on the specific parameter values, the bifurcations may be super- or subcritical; the branches may be connected or independent—this is determined by the fully nonlinear behaviour. Directly at each bifurcation, there exist harmonic waves of infinitesimal amplitude with propagation speed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_c=\text {Im}(\lambda )/k_c$$\end{document} and wavelength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_c=2\pi /k_c$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_c$$\end{document} corresponds to the neutral mode, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Re}(\lambda )(k_c)=0$$\end{document} . To gain knowledge about the branches of fully nonlinear travelling waves, one has to resort to numerical methods.

Finite-amplitude travelling waves correspond to steady states in the co-moving frame that moves with the wave velocity V. Setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c=c(x-Vt)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1=\rho _1(x-Vt)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2=\rho _2(x-Vt)$$\end{document} , we rewrite Eq. (8) as a boundary value problem with periodic boundary conditions and an integral condition that controls the mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^{-1}\int _0^{L} c(\chi )\,d\chi =c_0$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi =x-Vt$$\end{document} is the phase of the wave in the co-moving frame. We then employ numerical path-continuation techniques [43, 44] bundled in the software package auto07p [45] to determine families of travelling waves. As starting states, we use the analytically known small-amplitude harmonic waves, i.e. Eq. (11) with eigenvalue, eigenvector, and wavenumber of a neutral mode.

First, we consider the dry system ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} ) at fixed mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} (where the polar-nematic state does not yet exist) and examine the wave bifurcation where a branch of travelling waves emerges from the uniform nematic state. As determined in Sect. 4.1, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0>0.086$$\end{document} this uniform state is unstable in a certain range of wave numbers that depends on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} . When we fix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} and use the domain size L as control parameter, i.e. as principal continuation parameter, we encounter two disconnected branches of travelling waves that emerge where the dispersion relation in Fig. 4b has its two zero crossings. Figure 5 shows in panel (a) the propagation speed V, in panel (b) the amplitude of the density field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=\text {max}(c)-\text {min}(c)$$\end{document} , and in panel (c) the mean polar and nematic order parameter fields \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \rho _{1,2}\rangle = L^{-1}\int _0^L \rho _{1,2}(\chi ),\,d\chi $$\end{document} as a function of the domain size L. The wave with the larger [smaller] amplitude corresponds to a travelling density elevation [depression] both propagating in the uniform background, as shown in the insets of Fig. 5. When numerically integrating Eqs. (8) in time using a naive continuation scheme as described in Sect. 3, at all L the system converges to the state with the elevation (density hump, black line), as indicated by small star symbols in Fig. 5a, b. The travelling holes (red line) are never found in time simulations, which indicates that they are either unstable or linearly stable but with a very small basin of attraction.

One might expect that the two branches in Fig. 5 will be connected in a parameter region close to onset of the underlying wave instability when the unstable wave number band is very small. This is, however, for the present system not the case. Instead, a pre-existing branch of unstable nonlinear waves “touches down” to zero amplitude at the critical wavenumber at onset and roughly speaking splits into the two branches found above (not shown).Fig. 5. Two branches of travelling waves in a dry system ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=0$$\end{document} ) at fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=0.01$$\end{document} . Shown are a the wave speed V, b the amplitude of the density field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=\text {max}(c)-\text {min}(c)$$\end{document} as a function of the domain size L. Insets show the density profile at indicated points on the branches. Star symbols indicate stable wave states obtained by direct numerical time integration of Eq. (8). c gives the average polar and nematic fields \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \rho _{1,2}\rangle $$\end{document} for the solution branch highlighted by the star symbols

Second, we look at the dry system at the larger mean density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document} where also the uniform polar-nematic state exists. The linear stability results then indicate that overall three types of small-amplitude travelling wave states may emerge from the uniform states (one from purely nematic and two from mixed polar-nematic). The corresponding zero crossings of the dispersion relations (Fig. 4) correspond to three wave bifurcations and may give rise to up to three disconnected branches of finite-amplitude waves (the exact number depends on the fully nonlinear behaviour). At \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document} , we find three disconnected branches, see in Fig. 6a–c. The first one bifurcates at the critical k from the purely nematic state (red dashed line) and the other two at the critical k’s for the mixed polar-nematic state (solid black and dot-dashed blue line). The dot-dashed branch undergoes two saddle-node bifurcations; however, at large domain sizes three clearly distinguished travelling waves remain. They correspond to a single-hump wave (solid black line) as shown in Fig. 6(1,3), a double-hump wave (blue dot-dashed line) as shown in Fig. 6(2), and a travelling-hole state (red dashed line) as shown in Fig. 6(4). The double-hump state corresponds to a large travelling density cluster followed at a small distance by a small satellite peak. Direct numerical time integration of Eq. (8) using the naive continuation algorithm as described above only converges to single-hump waves (as indicated by star symbols in Fig. 6a–c) indicating their linear stability. The other types of travelling waves were not found in the time simulations.Fig. 6. Three disconnected branches of finite-amplitude travelling waves for a dry system (Ma \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=0$$\end{document} ) at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=0.01$$\end{document} . a Wave speed V, b amplitude of the density field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=\text {max}(c)-\text {min}(c)$$\end{document} and c average nematic order parameter field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \rho _2\rangle $$\end{document} as a function of domain size L. The four red filled circles marked (1) to (4) in panel (b) indicate loci of density profiles shown in the bottom panels. Star symbols indicate stable wave states obtained by direct numerical time integration of Eqs. (8)

For the purpose of simplicity, we do not discuss the more complicated bifurcation diagram for travelling waves in region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {II}$$\end{document} of Fig. 4d, where five wave bifurcations occur when using L as control parameter: four branches emerge from the mixed polar-nematic state and one from the purely nematic state.

Effect of the Marangoni flow on the travelling waves

Finally, we transition in a controlled manner from the dry to the wet system and analyse the effect of an increasing Marangoni flow on the propagation of polar-nematic waves. To do so, we follow the linearly stable nonlinear travelling waves found in Sect. 4.2 when increasing the Marangoni parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} from zero at various fixed domain sizes L. Figure 7a shows resulting branches of travelling wave states for five values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=50,60,70,100,200,300,400$$\end{document} at fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} (as in Fig. 5).Fig. 7a Speed of nonlinear travelling waves as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.09$$\end{document} at various fixed domain sizes L as indicated at the curves. Star symbols indicate linearly stable solutions found using direct numerical integration of Eqs. (8) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=400$$\end{document} . Square symbols with labels (1), (2) and (3) on the same branch indicate representative states that are shown in the panels on the right where black solid, red dashed and blue dotted lines show c, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _2$$\end{document} profiles, respectively. In (a), saddle-node bifurcation points are marked by green triangles, and the dashed red line corresponds to the loci of wave bifurcations in the (Ma,V)-plane where branches of travelling waves emerge from the purely nematic state [Eq. (9)]

Interestingly, the speed of the waves of larger periods, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L>100$$\end{document} , depends non-monotonically on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} . In complete agreement with the 2D case discussed in Fig. 2, the speed of the waves on the stable part of the branches increases with increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} , reaches a maximal value at a certain value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}_\textrm{max}(L)$$\end{document} , and then decreases again until the branch reaches the saddle-node bifurcation, marked by a green triangle. There, the stable part of the branch annihilates with the unstable part.

A typical scenario is for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=400$$\end{document} further illustrated by the density, polarisation, and nematic profiles on the right of Fig. 7: a stable wave in a dry system [state (1)] propagates at a speed of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=0.99$$\end{document} , and the speed slightly increases to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=1.03$$\end{document} at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=42$$\end{document} [state (2)] and then decreases to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=0.95$$\end{document} at the saddle-node bifurcation at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}=370$$\end{document} [state (3)]. The solitary waves (1) and (2) show pronounced peaks of density as well as of polar and nematic order and move in a nearly isotropic background, i.e. without polar or nematic order.

Note that the 2D travelling waves of maximal speed shown in Fig. 2c have an effective wavelength of about 400, which corresponds to the branch with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=400$$\end{document} in Fig. 7a. Stable waves are not found beyond the saddle-node bifurcation, as confirmed by direct numerical simulations of Eqs. (8), indicated by star symbols in Fig. 7a. At the saddle-node bifurcation, the branch of travelling wave states becomes unstable and turns towards smaller \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} . Finally, it terminates in a subcritical wave bifurcation on the dashed red line that indicates the loci of wave bifurcations in the (Ma,V)-plane. There, a continuous spectrum of branches of travelling waves of different L emerges from the purely nematic state [Eq. (9)]. It is obtained by continuing the zero crossing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Re}(\lambda )=0$$\end{document} of Eq. (12) in parameter space when gradually changing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=2\pi /k$$\end{document} .Fig. 8. Speed of travelling waves of different wavelength L as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document} . Star symbols indicate stable wave states found using direct numerical integration of Eq. (8) in a periodic domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=300$$\end{document} . The dashed red line corresponds to the loci of wave bifurcations in the (Ma,V)-plane where branches of travelling waves emerge from the mixed polar-nematic uniform state [Eq. (10)]. For the shown L, these points are indicated by empty square symbols

The speed V of the stable waves with the relatively small wavelength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=50$$\end{document} monotonically decreases with increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} until the saddle-node bifurcation is reached. If the Marangoni parameter is increased beyond the saddle-node point, in a time simulation the finite-amplitude wave relaxes to a uniform solution.

Interestingly, at larger mean densities, e.g. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0=0.12$$\end{document} and moderately small wavelength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\lesssim 300$$\end{document} , we find that the destruction of the wave by Marangoni flow may follow a different scenario. As shown in Fig. 8, e.g. for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=300$$\end{document} , the speed of stable waves can also monotonically increase with increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} and by a much larger percentage as at lower \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} ), until the point where the branch of the wave states ends in a wave bifurcation on the mixed polar-nematic uniform state. In this scenario, the wave amplitude gradually decreases until it eventually reaches zero at the supercritical wave bifurcation on the dashed red line. This is in contrast to the scenario of Fig. 7, where the finite-amplitude waves bifurcate from the trivial purely nematic state (given by the dashed line) via a sub-critical wave bifurcation.

Conclusion

We have studied the effect of self-induced Marangoni flow on the propagation of periodic and solitary density waves in a system of self-propelled agents with aligning interactions that move on the planar surface of a viscous liquid film. Thereby, the agents act as an insoluble surfactant in the purely entropic regime, i.e. with increasing density the surface tension of the liquid decreases linearly. This implies that emerging concentration gradients result in Marangoni forces that cause self-induced Marangoni flows. These influence the dynamics of the agents in such a way that advective flows in the liquid film and the dynamics of the active agents are intricately coupled. In this way, we have been able to quantitatively compare the dry case (no self-induced liquid flow) and the wet case (self-induced liquid flow).

The dynamics of a dry system with local nematic alignment is described by a truncated Boltzmann equation derived in Ref. [13, 14]. Here, we have adapted the model for wet systems by incorporating an advection term that results from the Marangoni flow. After establishing a simplified dynamical model that couples density, polar order parameter, and nematic order parameter fields, we have focused on the case of a long-wave approximation where the local flow field is proportional to the local density gradient. To characterise the degree of departure from the dry regime, we have introduced a Marangoni parameter that measures the relative strength of the Marangoni flow as compared to the self-propulsion speed of the agents.

We have then chosen agent interaction parameters that correspond to a regime where stable travelling density waves of mixed polar-nematic order are known to exist in the dry limiting case. As the Marangoni parameter is gradually increased, these waves broaden in the direction of propagation, while their speed either decreases or increases monotonically or changes non-monotonically depending on the mean concentration and wavelength. In addition, in the 2D case the density distribution homogenises in the direction transverse to the wave propagation, leading to the formation of almost perfectly parallel travelling density bands. We emphasise that such bands are well known to form in the absence of the Marangoni flow in dry as well as wet active-matter systems [4–8, 10]. As the strength of the Marangoni flow is further increased, the waves are destroyed and the system relaxes to a uniform polar or mixed polar-nematic state depending on the mean density. Depending on the wavelength of the nonlinear waves, this occurs either discontinuously via a saddle-node bifurcation of finite-amplitude travelling waves or continuously via a gradual decrease in the wave amplitude till it reaches zero at a wave bifurcation of a uniform state. The two scenarios correspond to cases of subcritical and supercritical wave bifurcations, respectively, that result in different non-equilibrium phase behaviour.

Our results have shown strong parallels between wet active-matter systems with local alignment and reaction–diffusion systems with advection [24, 25]. In both types of systems, the propagation of fronts and waves is similarly influenced by the addition of hydrodynamic flow, i.e. by “making the system wet”. Specifically, in autocatalytic chemical reaction systems, the Marangoni flow causes wave fronts to widen and increases their propagation speed. Similarly, in active-matter systems composed of self-propelled particles, we have observed that the propagation speed of polar-nematic waves increases with the Marangoni parameter, provided the wavelength is sufficiently large. The widening of the waves we have observed can also be attributed to the Marangoni advection flow, which moves away from regions of higher concentration. This analogy between the two systems can be understood from a phenomenological perspective: in the limit of small density variations, the dynamics of active-matter systems can be described by a closed set of reaction–advection–diffusion equations governing the dominant order parameter fields.