Spatial Organization of Active Particles with Field Mediated Interactions

TL;DR

This paper investigates how active particles interacting with a fluctuating field and driven by nonequilibrium noise can self-organize into complex spatial patterns, revealing the influence of external driving on pattern formation.

Contribution

It introduces a model of active particles with internal states and analyzes how nonequilibrium noise influences their spatial organization and phase behavior.

Findings

Identification of a phase diagram with patterned and unpatterned states.

Demonstration of nonequilibrium drive's role in pattern formation.

Quantitative analysis of pattern boundary dependence on parameters.

Abstract

We consider a system of independent point-like particles performing a Brownian motion while interacting with a Gaussian fluctuating background. These particles are in addition endowed with a discrete two-state internal degree of freedom that is subjected to a nonequilibrium source of noise, which affects their coupling with the background field. We explore the phase diagram of the system and pinpoint the role of the nonequilibrium drive in producing a nontrivial patterned spatial organization. We are able, by means of a weakly nonlinear analysis, to account for the parameter-dependence of the boundaries of the phase and pattern diagram in the stationary state.

Equations

$$H_0=\int d^{2}x\left[\frac{r}{2}{\varphi }^2+\frac{c}{2}(\nabla \varphi {)}^2\right],$$

$$H_{\mathrm{int}}=\sum _{k=1}^N\frac{B}{2}{\left(\varphi ({\mathbf{r}}_k)-S_k{\varphi }_0\right)}^2,$$

$${\partial }_t\varphi (\mathbf{x},t)=-\Gamma \frac{\delta H}{\delta \varphi (\mathbf{x},t)}+\sqrt{2\Gamma T}\zeta (\mathbf{x},t),$$

$$⟨\zeta (\mathbf{x},t)\zeta ({\mathbf{x}}^{\prime },t^{\prime })⟩=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime }),$$

$$\frac{d{\mathbf{r}}_k}{dt}=-\mu \frac{\partial H}{\partial {\mathbf{r}}_k}+\sqrt{2\mu T}{\mathbf{\xi }}_k(t),$$

$$⟨{\xi }_k^{\alpha }(t){\xi }_{\ell }^{\beta }(t^{\prime })⟩={\delta }_{\alpha \beta }{\delta }_{k\ell }\delta (t-t^{\prime }),$$

$$S_k=-1\underset{\text{leavevmode}\gamma \text{leavevmode}}{\overset{\text{leavevmode}\alpha \text{leavevmode}}{\rightleftharpoons }}{S}_k=+1.$$

$$\displaystyle\begin{split}\phi_{ij}(t+\Delta t)=&\phi_{ij}(t)-\Delta t\big{[}r\phi_{ij}(t)-\nabla^{2}\phi_{ij}(t)\\ &+B\sum_{k=1}^{N}(\phi_{\bm{r}_{k}}(t)-S_{k}\phi_{0})\big{]}+G(0,2\Delta t),\end{split}$$

$${\nabla }^2{\varphi }_{ij}\equiv {\varphi }_{i+1,j+1}+{\varphi }_{i-1,j}+{\varphi }_{i,j+1}+{\varphi }_{i,j-1}-4{\varphi }_{ij},$$

$$\Delta H_{{\mathbf{r}}_k\to {\mathbf{r}}_k^{\prime }}=\frac{B}{2}({\varphi }_{i_k^{\prime },j_k^{\prime }}-{\varphi }_{i_k,j_k})({\varphi }_{i_k^{\prime },j_k^{\prime }}+{\varphi }_{i_k,j_k}-2{\varphi }_0{S}_k).$$

$$P_{(i_k,j_k)\to (i_k^{\prime },j_k^{\prime })}=\mu \Delta t\exp \left(-\frac{\Delta H_{{\mathbf{r}}_k\to {\mathbf{r}}_k^{\prime }}}{2}\right).$$

$$\sum _{(i_k^{\prime },j_k^{\prime })}{P}_{(i_k,j_k)\to (i_k^{\prime },j_k^{\prime })}+P_k^f<1,$$

$$P_k^f(\uparrow \to \downarrow )=\eta \Delta t\exp \left(-\frac{H_{k\downarrow }-H_{k\uparrow }}{2}\right),$$

$$P_k^f=\{\begin{array}{ll}\gamma \Delta t& \text{mbox}if{S}_k=+1\\ \alpha \Delta t& \text{mbox}if{S}_k=-1.\end{array}$$

$$\begin{array}{rl}{f}_{\mathrm{MF}}=\frac{r}{2}{\varphi }^2& +B\frac{\rho +\psi }{4}(\varphi -{\varphi }_0{)}^2+B\frac{\rho -\psi }{4}(\varphi +{\varphi }_0{)}^2\\ & +\frac{\rho +\psi }{2}\ln \frac{\rho +\psi }{2}+\frac{\rho -\psi }{2}\ln \frac{\rho -\psi }{2},\end{array}$$

$${\partial }_t{\rho }^{+}$$

$${\partial }_t{\rho }^{-}$$

$${\partial }_t\varphi$$

$${\partial }_t\rho$$

$$\begin{array}{rl}{\partial }_t\psi & =\mu {\nabla }^2\psi +\mu B\nabla \cdot \left[(\psi \varphi -\rho {\varphi }_0)\nabla \varphi \right]\\ & -\omega \psi +(2s-1)\omega \rho ,\end{array}$$

$${\partial }_t\varphi$$

$$\displaystyle\begin{split}\rho_{ij}^{+}(t+\Delta t)=\rho_{ij}^{+}(t)+&\Delta t\bigg{[}\mu\nabla^{2}\rho_{ij}^{+}\\ &+\mu B\,\nabla_{x}(\rho^{+}_{ij}(\phi_{ij}-\phi_{0})\nabla_{x}\phi_{ij})\\ &+\mu B\,\nabla_{y}(\rho^{+}_{ij}(\phi_{ij}-\phi_{0})\nabla_{y}\phi_{ij})\\ &-\gamma\rho^{+}_{ij}+\alpha\rho^{-}_{ij}\bigg{]}(t),\end{split}$$

$$\displaystyle\begin{split}\phi_{ij}(t+\Delta t)=\phi_{ij}(t)+&\Delta t\big{[}\nabla^{2}\phi_{ij}-r\phi_{ij}\\ &-B\rho^{+}_{ij}(\phi_{ij}-\phi_{0})\\ &-B\rho_{ij}^{-}(\phi_{ij}+\phi_{0})\big{]}(t),\end{split}$$

$${\nabla }_x{g}_{ij}\equiv \frac{1}{2}(g_{i+1,j}-g_{i-1,j}),$$

$${\nabla }_y{g}_{ij}\equiv \frac{1}{2}(g_{i,j+1}-g_{i,j-1}),$$

$${\rho }_h$$

$${\psi }_h$$

$${\varphi }_h$$

$$\tilde{r}=r+B{\rho }_0={\tilde{\xi }}^{-2}$$

$$M=\left(\begin{array}{ccc}-\mu k^2& 0& -\mu B{\rho }_0{\varphi }_0(2s-1)(B{\rho }_0/\tilde{r}-1)k^2\\ (2s-1)\omega & -\mu k^2-\omega & -\mu B{\rho }_0{\varphi }_0[(2s-1{)}^{2}B{\rho }_0/\tilde{r}-1]k^2\\ -B^2{\rho }_0{\varphi }_0(2s-1)/\tilde{r}& B{\varphi }_0& -k^2-\tilde{r}\end{array}\right),$$

$$q_1$$

$$q_2$$

$$(2s-1{)}^2>\frac{1}{B{\varphi }_0^2}{\left(1+\frac{B{\rho }_0}{r}\right)}^2\left(1+\frac{r}{B{\rho }_0}\right),$$

$${\sigma }_1$$

$${\sigma }_2$$

$${\sigma }_3$$

Full text

I Introduction

Interactions between objects mediated by an elastic medium are ubiquitous in the soft-matter realm. Examples include surfaces, colloids or proteins in soft-matter media such as critical binary mixtures Fisher and de Gennes (1978); Hertlein et al. (2008), liquid crystals Ajdari et al. (1991); Poulin et al. (1997), capillary interfaces Nicolson (1949); Hu and Bush (2005) and bio-membranes Goulian et al. (1993); Dan et al. (1993); Dommersnes and Fournier (1999); Bitbol et al. (2012); Van Der Wel et al. (2016). Effective interactions, already at the level of equilibrium systems, may lead to a rich phase behavior Poulin et al. (1997); Brown et al. (2000); Noguchi and Fournier (2017). Yet relaxing the constraint of detailed balance is expected to allow for an even richer phenomenology. Such nonequilibrium examples include biophysical agents such as bacteria or cells. These are good examples of active particles whose chemical action on the medium Miller and Bassler (2001), or whose hydrodynamic interactions with the medium Stenhammar et al. (2017), give birth to mediated (non local) interactions between particles that cannot be encapsulated in an effective energy picture. Down to even smaller scales, recent experiments on reactive proteins suggest interesting collective behavior made possible by the induced deformation of the lipid bilayer in which proteins are embedded in Fribourg et al. (2014); Sumino et al. (2014). This class of systems also encompasses synthetic objects that actively shape the medium in which they are confined. The latter include liquid crystal colloids endowed with a tunable degree of freedom Evans et al. (2013), whose interactions have been experimentally explored Criante et al. (2013) with a view to controling self-assembly at the micrometer scale.

The goal of the present work is to build a generic model belonging to the same family of systems as those just described, in particular by incorporating an externally switchable (active) degree of freedom. To do so, we shall retain the physical ingredients shared by these systems.

First, we resolve particles at the individual level and we assume their individual motion to be diffusive. We further endow our particles with a two-state, spin-like, internal degree of freedom that can be switched by an external drive (this is where activity will come into play). Second, we choose to describe the embedding medium by a coarse-grained field. Third, we consider a coupling between the medium and the degrees of freedom –both spin and position– of the particles, that will lead to mediated interactions. In order to focus on the latter, we omit from our model any direct interaction (hard-core, attractive or else) between particles. The dynamics of the medium itself is assumed to be local and purely relaxational (model A-like). The out-of-equilibrium nature of the system comes from the active conformation switch of the particle that breaks detailed balance. Beyond proteins in biomembranes, such externally driven conformational changes can also be found in synthetic soft-matter systems (see Fribourg et al. (2014); Grawitter and Stark (2018); Chen (2004) for recent references).

Even with the simplifying assumptions that have led to our model, our particles do evolve far from equilibrium, and no free-energy based method is available. Predicting collective phenomena thus requires to implement a variety of approaches, both numerical (Monte Carlo) and analytical (mean-field equations, noiseless reaction-diffusion equations). We present the details of the model and in particular its key parameters in Sec. II. Its stationary phase diagram is explored in Sec. III, by means of Monte Carlo dynamical simulations, both for our active system and its equilibrium counterpart (that we properly define). A variety of patterned phases emerge in some regions of our parameter space. The subsequent mean-field analysis of Sec. IV allows us to understand the phase boundaries of the phase diagram given by our simulations. This very good qualitative (and good quantitative) agreement between the solution of the mean-field partial differential equations (PDE) and dynamical simulations suggests that the mean-field approach might also prove powerful to describe the physical nature of our patterned phases. Thus, in Sec. V, we embark into a linear and a nonlinear analysis of these equations, which allows us to describe the pattern content of our physical problem. Particular emphasis is placed on the extra mathematical difficulty of dealing with a conserved mode in a pattern forming system, a question that was hitherto sidelined in the existing literature of active inclusions in membranes Reigada et al. (2005a, b). As it will turn out, this is the existence of the conserved mode (expressing that particles are conserved, regardless of their internal spin state) that gives birth to a rich phenomenology of patterns. Our analysis of such patterns will draw from recent theoretical work Riecke (1992a); Matthews and Cox (2000); Cox and Matthews (2003); Winterbottom (2006). In our final two sections, we discuss the role of the nonequilibrium drive in producing patterned phases. To this end we introduce, in Sec.VI, a model that interpolates between the active system and its equilibrium counterpart. This allows us to probe the robustness of patterns with respect to a partial restoration of reversibility (via tunable coupling to the same thermal bath the particles and the field are in contact with). To further pinpoint at the microscopic level the processes by which entropy is created in our active system, we establish a spatial map of entropy production that we superimpose to the patterns we obtain. The various regimes observed throughout our simulations can be reasonably rationalized using this versatile entropy production as a quantitative indicator. This last piece of our analysis can be found in Sec. VII.

II Model

The first ingredient of our model is a fluctuating field $\phi$ standing for the surrounding medium in which our particles are embedded. Our analysis is confined to a two-dimensional medium. We choose to use a free field with a Gaussian Hamiltonian endowed with the following features: the value of the field at rest and without particles is [math], and the field has a finite correlation length $\xi$. The Hamiltonian of the field then reads

[TABLE]

with $\xi=\sqrt{c/r}$. We assume the medium is in contact with a thermostat at temperature $T$. We further assume a separation of scales between the medium constituents and the particles, so that the medium can be described by a continuous field on a continuous space. However we retain the individual localized nature of the particle which we describe by their position $\bm{r}_{k}$. The value $\phi(\bm{r}_{k})$ of the field at the position $\bm{r}_{k}$ of particle $k$ is elastically constrained to the value $\pm\phi_{0}$ by the internal degree of freedom $S_{k}=\pm 1$ of the particle. This leads us to use the following interaction energy between $N$ particles and the field

[TABLE]

where $B$ is the strength of the particle-field coupling. Note that, as discussed in the introduction, particles experience no direct interactions (not even hard-core repulsion). The Ising spin variable $S_{k}$ refers to the two internal states the particle is assumed to be found in. More realistic models will of course be system-dependent. For instance, to describe conically shaped proteins that locally constrain lipid membrane curvature, a Helfrich Hamiltonian Helfrich (1973) should be used instead of Eq. (1). A description of the membrane thickness with the Landau-Ginzburg Hamiltonian (1) is perhaps better adapted to the description of protein-protein interactions experiencing hydrophobic mismatch interactions and coarse-grained packing interactions, already existing in pure one-component lipid bilayer West et al. (2009). Once energy functions are specified, we turn to the question of how to implement dynamical evolution. Regarding the background field itself, discarding possible conservation laws or hydrodynamic interactions (either or both could prove relevant in a variety of physical systems), we resort to a purely relaxational dynamics consistent with the contact to a thermal bath at temperature $T$:

[TABLE]

where $H=H_{0}+H_{\mathrm{int}}$, $\Gamma$ is a mobility coefficient, $T$ is the temperature in energy units and $\zeta(\bm{x},t)$ a Gaussian white noise with zero average. As far as particles are concerned, their (low Reynolds) motion is described by an overdamped Langevin equation:

[TABLE]

where $\mu$ is a mobility coefficient, and the $\xi_{k}^{\alpha}(t)$ are the components of independent Gaussian white noises with zero average. We use the simplifying assumption that $\bm{\xi}_{k}$ and $\zeta$ are independent. So far, at fixed spin variables, our dynamics is consistent with detailed balance. The nonequilibrium drive will arise from the dynamics the spins are endowed with. With an external source of energy (such as photons or ATP in biological systems) in mind, we introduce temperature and state independent flipping rates $\alpha$ and $\gamma$:

[TABLE]

For the purpose of benchmarking genuinely nonequilibrium effects, we shall later introduce a detailed-balance preserving spin-flip dynamics. The final simplifying step is to work in terms of dimensionless parameters. We introduce a characteristic size $a$ which will be used to spatially discretize the field $\phi$, we normalize energies by $T$, times by $a^{2}/(\Gamma c)$ and we absorb $c$ in a redefinition of the field $\phi$. We thus carry out the replacements $\bm{x}/a\to\bm{x}$, $\Gamma ct/a^{2}\to t$, $c\phi^{2}/T\to\phi^{2}$, $c\phi_{0}^{2}/T\to\phi_{0}^{2}$, $a^{2}r/c\to r$, $B/c\to B$, $T\mu/(\Gamma c)\to\mu$, $a^{2}\alpha/(\Gamma c)\to\alpha$ and $a^{2}\gamma/(\Gamma c)\to\gamma$. In a nutshell, rescaling time, space, fields and constants boils down to $a=c=T=\Gamma=1$. Our model being now defined, we present the results of our numerical Monte Carlo-based exploration of its properties.

III Numerical simulations

We perform Monte Carlo dynamical simulations on a two dimensional square lattice of size $L_{x}\times L_{y}$ with periodic boundary conditions. The Gaussian field is defined on each site $(i,j)$ and takes continuous real values $\phi_{ij}$. The field $\phi$ evolves according to the explicit stochastic Euler scheme corresponding to the dimensionless form of Eq. (3):

[TABLE]

where the discrete Laplacian is defined as

[TABLE]

$\bm{r}_{k}\equiv(i_{k},j_{k})$ being the position of the particle $k$ on the lattice, and $G(0,2\Delta t)$ being a Gaussian variable of mean [math] and variance $2\Delta t$.

At each time step, we update the field, and then the particles’ positions. The particles lie on the lattice sites. Since our model involves non-interacting particles, we a priori allow for several particles to occupy the same lattice site. We implement a tower sampling algorithm Krauth (2006) to choose what action a particle should do, namely, either jump on a neighboring site, or stay on the same site or flip its spin. The total energy variation when the particle $k$ moves from site $(i_{k},j_{k})$ to $(i^{\prime}_{k},j^{\prime}_{k})$ is given by

[TABLE]

The following jump probability $P_{(i_{k},j_{k})\to(i^{\prime}_{k},j^{\prime}_{k})}$ between times $t$ and $t+\Delta t$ implements a discrete version of the Langevin equation (5):

[TABLE]

According to our model, the spin flip probability of particle $k$, $P^{f}_{k}$, is fixed by the rates $\alpha$ and $\gamma$, except when we consider detailed-balance preserving flipping rates for the purpose of comparison to the out-of-equilibrium case. Each case will be specified below. We take $\Delta t$ small enough to ensure that the probabilities verify

[TABLE]

then the probability $P^{n}_{k}$ that particle $k$ neither jumps nor flips is given by $P^{n}_{k}=1-[\sum_{(i^{\prime}_{k},j^{\prime}_{k})}P_{(i_{k},j_{k})\to(i^{\prime}_{k},j^{\prime}_{k})}+P^{f}_{k}]$. We take $\Delta t=2\times 10^{-5}$.

III.1 Equilibrium benchmark

Though we are mostly interested in the active system, for the purpose of discussion and comparison we first study the system with spin flips that preserve the detailed balance condition. This is useful to sort out generic collective phenomena already present in our equilibrium Hamiltonian from those induced by activity. In this equilibrium benchmark system, the probability of a spin flip, say from spin up to spin down, between $t$ and $t+\Delta t$ is

[TABLE]

where $\eta$ is an inverse time-scale and where $H_{k\uparrow}$ (resp. $H_{k\downarrow}$) refers to the energy of the system when spin $k$ is up (resp. down).

When endowed with this equilibrium reversible dynamics the particles+field system already displays a nontrivial phase diagram explored with dynamical Monte Carlo simulations (see Fig. 1). When the coupling $B$ with the field is small, a paramagnetic fluid (the average value of the spins is zero) is observed. When the coupling is increased, the system displays a paramagnetic-ferromagnetic transition. When the density is small, further increasing the coupling yields a phase separation between a paramagnetic gas and a ferromagnetic liquid (Fig.1c). As a consistency check, we verified that the equilibrium phase diagrams do not depend on dynamical parameters. The specific simulations shown in Fig. 1 were performed with $\mu=5$, $\eta=1$.

III.2 Active system

In the active system, the flipping probabilities corresponding to Eq. (7) are given by

[TABLE]

If flipping rates are symmetric ($\alpha=\gamma$) the average magnetization is zero and the system cannot develop a homogeneous ferromagnetic state. In the asymmetric case, it is convenient to define the total flip rate $\omega=\alpha+\gamma$ and the mean fraction $s=\alpha/\omega$ of spin-up particles in steady state. The flipping rates are thus given by $\alpha=s\omega$ and $\gamma=(1-s)\omega$.

In the following, we explore the phase diagram of the system for $s=1/2$ (Fig. 2, top). When the coupling $B$ to the field is weak, the system remains homogeneous (and paramagnetic). At low densities, when increasing $B$, finite size clusters of both magnetization appear (see Fig. 3a). At higher densities, the phenomenology becomes richer. Increasing $B$ from the homogeneous phase, macroscopic stripes of both magnetization (Fig. 3b) are observed. As $B$ is further increased, the stripes harbor the continuous nucleation of small lumps of particles of opposite magnetization (Fig. 3c). These lumps grow, drift, then merge with adjacent bands of same magnetization. Increasing $B$ again, the proliferation of lumps leads to a system of micro-clusters (Fig. 3d). In the patterned phase (stripes or clusters), increasing the flipping frequency $\omega$ yields local mixing of the spins, which results in the homogeneization of the whole system (Fig. 4, top).

If now asymmetric flipping rates ($s\neq 1/2$) are considered, the phase diagram features similar transitions. The homogeneous phase is however ferromagnetic since the mean number of spins up and spins down is different. In addition, because of the breaking of the up-down symmetry, hexagonal patterns can be observed (see Fig. 5).

In the following section, we work out a mean-field analysis which predicts the transition between different regimes. The appearance of lumps lies at the transition between two different pattern forming regimes.

IV Mean-field behavior

IV.1 Equilibrium free energy

We briefly treat the equilibrium case where particles are allowed to flip while respecting detailed balance (equilibrium benchmark). Our goal is to determine the phase diagram of the system. We first identify a conserved field, $\rho=\rho^{+}+\rho^{-}$, and a non conserved field $\psi=\rho^{+}-\rho^{-}$. From the Hamiltonian $H$, we write down a mean-field free energy density:

[TABLE]

where the first three terms are directly inferred from the energy functional $H$, while the last two ones reflect the particles’ entropy. Since $\rho$ is the only conserved quantity, we minimize $f_{\mathrm{MF}}$ with respect to $\phi$ and $\psi$. We obtain $\phi=B\phi_{0}\psi/(r+B\rho)$ and $\psi=\rho\tanh(B\phi_{0}\phi)$. This imposes a self-consistent equation on $\phi=B\phi_{0}\rho\tanh(B\phi_{0}\phi)/(r+B\rho)$ similar to what is obtained for the magnetization in the Ising model. Searching for homogeneous phases, yields either $\phi=0$ (paramagnetic phase), or $\phi\neq 0$ (ferromagnetic phase). At low values of $B$, the system is uniform and there is a continuous paramagnetic–ferromagnetic transition at $B_{c}^{(\mathrm{MF})}=(1+\sqrt{1+4r\,\phi_{0}^{2}/\rho})/(2\phi_{0}^{2})$. At higher values of $B$, we numerically solve the double tangent construction on $f_{\mathrm{MF}}(\rho)$ (already minimized with respect to $\psi$ and $\phi$). The system undergoes a phase separation between a low density paramagnetic phase and a high density ferromagnetic phase. These mean-field predictions correspond to the continuous lines of Fig. 1 while the results of the Monte Carlo simulations are indicated by the dashed lines. We have checked that the agreement is all the better as we are working at large $\phi_{0}$.

IV.2 Mean field dynamics

We consider now the original model of interest where flips are fixed by an external and independent source of energy. Out of equilibrium, we can no longer rely on the free energy to construct the phase diagram. Since particles execute Brownian motions, we consider the noiseless limit of the Dean–Kawasaki equations Dean (1996) for the up and down particle densities. Taking spin exchange into account (and neglecting the corresponding Poisson noise as well), we arrive at the deterministic evolution equations for $\rho^{\pm}$:

[TABLE]

It will prove convenient to write these equations in terms of the conserved field $\rho$, and of the non-conserved field $\psi$. We also parametrize the rates $\alpha$ and $\gamma$ by means of $\omega=\alpha+\gamma$ and $s=\alpha/\omega$ (the latter being the steady-state fraction of spin up particles). The dynamical evolutions of the fields then read

[TABLE]

These three equations are the starting point of our analysis of the patterns that form in the steady-state of our system. It is important to note that $s=\frac{1}{2}$ will play a special role because then these equations are invariant upon the up-down symmetry $(\rho,\psi,\phi)\to(\rho,-\psi,-\phi)$. Before we embark in a detailed analytical study of their pattern content, we begin with a numerical solution of these nonlinear coupled PDEs.

IV.3 Numerical solution of the coupled partial differential equations

We shall show that the numerical solution of the nonlinear coupled PDEs (which are noiseless) is relevant to analyze the stochastic simulations to the extent that phases and phase boundaries are quite faithfully captured. The coupled PDEs are solved on a lattice of size $L_{x}\times L_{y}=150\times 150$. The three fields $\rho^{+}$, $\rho^{-}$ and $\phi$ are discretized in time and space; an explicit Euler scheme to update the three fields is implemented. The explicit Euler scheme is easy to implement and it converges in the domains of the phase diagram we are interested in. Our discretized equations take the following form:

[TABLE]

[TABLE]

and the discretized equation on $\rho^{-}$ is formally identical to the discretized equation on $\rho^{+}$ up to the exchange $\rho^{+}\leftrightarrow\rho^{-}$, $\phi_{0}\to-\phi_{0}$, $\alpha\leftrightarrow\gamma$. The discrete spatial derivatives of any field $g_{ij}$ are defined as

[TABLE]

and the Laplacian has already been defined in Eq. (9). We confirm that different initial conditions lead to same stationary density profiles. We check the conservation of total density, namely $(L_{x}L_{y})^{-1}\sum_{ij}\rho_{ij}=\rho_{0}$, along with the positivity of $\rho^{+}$ and $\rho^{-}$ on each site.

To ease comparison of the PDE solution with the Monte Carlo simulation, the PDE phase diagram is plotted Fig. 4 (bottom) for the same physical parameters as those of the Monte Carlo results of Fig. 4 (top). The results of the PDEs numerical solution match the results of the Monte Carlo simulations. The solution of the PDEs is also in good agreement with the predictions obtained from a weakly nonlinear analysis in Sec. V.2: we can observe either homogeneous patterns (Fig. 6, top), or spatially localized patterns (Fig. 6, bottom), depending on the parameters. The coming section is devoted to an analysis of these patterns.

V Pattern analysis

V.1 Linear stability analysis

By resorting to a linear stability analysis (LSA), the range of parameters for which a uniform stationary state is destabilized can be found. While LSA tells us about the first unstable mode, the question of which are the selected modes that eventually build up into patterns requires a full analysis of the nonlinear equations. The homogeneous and stationary solution to Eqs. (19), (20), (21) is characterized by the following values of the fields

[TABLE]

with

[TABLE]

where $\tilde{\xi}$ is the renormalized correlation length of the field $\phi$. We set $\rho_{1}=\rho-\rho_{h}$, $\psi_{1}=\psi-\psi_{h}$ and $\phi_{1}=\phi-\phi_{h}$ and we expand Eqs. (19), (20), (21) to linear order in the $\rho_{1}$, $\phi_{1}$, $\psi_{1}$ fields. We expand the fields in Fourier modes $\sim e^{i\bm{k}\cdot\bm{x}}$ and we arrive at a linear system for the Fourier components $\partial_{t}(\rho_{1},\psi_{1},\phi_{1})^{T}=M(\rho_{1},\psi_{1},\phi_{1})^{T}$

with

[TABLE]

with $k=\|\bm{k}\|$. The eigenvalues of $M$ can be shown to be always real which excludes oscillating patterns close to the threshold. We denote them by $\sigma_{i}$, with $\sigma_{1}<\sigma_{2}<\sigma_{3}$. Solving $\det M(k)=0$ yields the modes for which temporal growth is marginal. In practice, we have $\det M=-\mu^{2}k^{2}Q(k^{2})$, where $Q(X)=X^{2}+q_{1}X+q_{2}$ is degree 2 polynomial, with

[TABLE]

Three different physical cases must be distinguished, depending on the roots $X_{-}$, $X_{+}$ of $Q$.

• •

case (i): $q_{1}\geq 0$ and $q_{2}\geq 0$, or $q_{2}<q_{1}^{2}/4$, then $Q$ has no real positive roots. One can show that the three eigenvalues of $M$ are negative: the homogeneous state is stable.

• •

case (ii): $q_{2}<0$ then $Q$ has only one positive root $X_{+}$. We set $k_{+}\equiv X_{+}^{1/2}$. In this regime, unstable modes go from $k=0$ to $k=k_{+}$ and we numerically observe either coarsening, or pattern formation depending on the sign of $q_{1}$, as shown in Fig. 8. When $\omega$ is non zero, we sit in regime (ii) where $q_{2}<0$ is equivalent to

[TABLE]

which surprisingly does not depend upon the dynamical parameters. Physically, the instability comes from the frustrated field $\phi$ whose value at rest and without particles is [math], different from $\phi_{h}$ in the presence of particles with non-symmetric flipping rates. This regime is referred to as type $\mathrm{II}_{s}$ in Cross and Hohenberg (1993).

• •

case (iii): $q_{1}^{2}/4\geq q_{2}>0$; patterns appear at finite wavelength (referred to as type $\mathrm{I}_{s}$ in Cross and Hohenberg (1993)). We have $X_{-},X_{+}>0$ and we set $k_{\pm}\equiv X_{\pm}^{1/2}$. The eigenvalue $\sigma_{3}(k)$ is positive for $k\in[k_{-},k_{+}]$. At the onset of instability, the only growing mode is indexed by $k_{c}$ with, at the threshold, $k_{c}=k_{-}=k_{+}$.

Note also that when $\omega=0$ we observe an equilibrium coarsening of the two populations of particles, under the condition $q_{1}<0$ (which is equivalent to sitting in the equilibrium ordered phase).

In summary, (ii) and (iii) are the two regimes where the homogeneous state is destabilized. For nonzero flipping rates, the only way to transition from regime (i) to regime (ii), or from regime (iii) to (ii), is by changing the equilibrium parameters, namely $r$, $B$, $\rho_{0}$, $\phi_{0}$ and $s$. By contrast, at fixed equilibrium parameters, we transition from regime (i) to regime (iii) by changing $\omega$ or $\mu$. In the following, we will focus on the transition caused by a change in the dynamics, and consequently, on instabilities starting at finite wavelength.

We begin our analysis with the simpler $s=1/2$ symmetric case, where the number of particles of each spin is identical in the steady-state. This ensures, after Eq. (33), that we are always in the pattern forming regime (iii). The matrix $M$ is now block diagonal and eigenvalues can be cast in a compact form:

[TABLE]

with

[TABLE]

For the purpose of discussion we use $\omega=\alpha+\gamma$ as the control parameter. Physically, we recall that for high flipping rates, the system remains homogeneous since particles locally efficiently mix, whereas for $\omega=0$, the system undergoes an equilibrium coarsening (see phase diagram Fig. 4). Solving $q_{1}^{2}-4q_{2}=0$ yields a critical value of $\omega$:

[TABLE]

below which the homogeneous system is no longer stable. To study the system close to this transition, we write $\omega=\omega_{c}-\varepsilon^{2}$, where the distance to the threshold $\varepsilon^{2}>0$ becomes our control parameter. Since we sit in regime (iii), destabilization occurs at a mode $k_{c}=k_{\pm}>0$ when $\varepsilon=0$. Thus, when $\omega<\omega_{c}$, $\sigma_{3}(k)\geqslant 0$ for $k\in[k_{-},k_{+}]$, with

[TABLE]

and where $\Omega=\omega/\mu$. The condition of existence of the $k_{\pm}$ modes (namely that $X_{\pm}$ are real) is given by

[TABLE]

At $\omega=\omega_{c}$ equality is achieved in Eq. (40) and this allows us to infer the critical wavelength of patterns $\lambda_{c}\equiv 2\pi/k_{c}=2\pi(\tilde{r}\omega_{c}/\mu)^{-1/4}$. This suggests that close to the threshold the patterns spatial periodicity $\lambda_{p}$ could be the combination $\lambda_{p}=2\pi/k_{p}=2\pi(\tilde{r}\omega/\mu)^{-1/4}$.

This prediction has been checked in simulations of a quasi 1D system of size $1000\times 10$ to force pattern formation along one direction, hence allowing us to achieve a good precision on the wavelength. We note on Fig. 9 that the prediction on the pattern periodicity $\lambda_{p}$ applies beyond the pattern formation threshold. Interestingly enough, $\lambda_{p}$ can be expressed as the geometric mean of the renormalized correlation length $\tilde{\xi}$ of the Gaussian field $\phi$ in presence of inclusions and of the typical diffusion length $\ell_{d}\sim\sqrt{\mu/\omega}$ of a particle between two flips. The formula $\lambda_{p}\sim(\tilde{\xi}\ell_{d})^{1/2}$, expresses, at the level of a cluster, the balance between accretion via interactions vs. loss by diffusion. It would certainly be interesting to see $\lambda_{p}$ emerge from a handwaving argument. Finally, it is worth noticing that close to threshold the selected wavelength does not depend upon the field mobility: it is only the particles’ mobility with respect to the spins’ flipping rate that matters. An estimate of the diffusion time of a particle over a characteristic correlation length $\tilde{\xi}$ of the field is $t_{d}\approx\tilde{\xi}^{2}/(2\mu)$. On the other hand, the correlation time $t_{\phi}$ of the field over a scale $\tilde{\xi}$ is given by $t_{\phi}\approx 1/(2\tilde{r})$. Hence particles are fast with respect to the field when $t_{d}\ll t_{\phi}$, or $1\ll\mu$. In Fig. 9, one can indeed see that the selected wavelength does not change, whether particles are slow or fast with respect to the field. We now turn to an analysis of the patterns that form beyond threshold.

V.2 Weakly nonlinear analysis

In order to gain insight into nonlinear effects at $s=1/2$, we can derive an amplitude equation for the fields by extending the approach of Swift and Hohenberg Swift and Hohenberg (1977). A direct, though naive, way of proceeding would be to extract the equations for the relevant fields for which we can find modes with exponential growth. In particular, in the basis where $M$ is diagonal, there is only one direction (corresponding to eigenvalue $\sigma_{3}$) along which we observe the temporal growth of the Fourier modes (see Fig. 10).

The eigen-fields are given by the LSA. By denoting $D=P^{-1}MP$ where the transformation matrix $P$ writes

[TABLE]

with

[TABLE]

we define $(U,V,W)^{T}=P^{-1}(\rho_{1},\psi_{1},\phi_{1})^{T}$, with $\rho_{1}$, $\psi_{1}$ and $\phi_{1}$ which are no longer infinitesimally small perturbations. The new fields $(U,V,W)$ now verify $U\sim e^{\sigma_{1}t}$, $V\sim e^{\sigma_{2}t}$ and $W\sim e^{\sigma_{3}t}$. However, we have to be careful: the existence of the conserved quantity $\rho$ implies that the mode $k=0$ is marginal for the field $U$ and has to be taken into account Coullet and Fauve (1985); Riecke (1992b); Matthews and Cox (2000). Indeed nonlinear terms couple modes $k\to 0$ with modes $k\sim k_{c}$ in products $O(UW)$. In the absence of a marginal growth, it would be correct to focus only on modes with exponential growth around $k_{c}$ and we could expand $\sigma_{1}$, $\sigma_{2}$ and $\sigma_{3}$ around $\varepsilon=0$ and $k\sim k_{c}$. This would lead to neglect terms of the form $V^{q}W^{p}$ and $U^{q}W^{p}$ (for $q\geqslant 1$ and $p\geqslant 0$) since they exponentially go to [math] when $k\sim k_{c}$. Interestingly, one of the erroneous conclusions we would arrive at is that square patterns could never be stable, in conflict with observations of PDEs’ solution. Now keeping both relevant fields $U$ and $W$, the evolution equations read

[TABLE]

where $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ are nonlinear operators that couple $U$ and $W$. To lowest order in $\varepsilon$, we find that $\mathcal{N}_{1}$ contains terms $\sim O(W^{2})$ and that $\mathcal{N}_{2}$ contains terms $\sim O(UW)$. Thus, $U$ will saturate to $O(W^{2})$, which renormalizes $O(W^{3})$ terms in $\mathcal{N}_{2}$. Of course, this previous analysis holds for the case $s=1/2$, where equations are invariant upon the $(\rho,\psi,\phi)\to(\rho,-\psi,-\phi)$ symmetry. If $s\neq 1/2$, new terms appear in nonlinear equations (44) and (45). To lowest order, new terms in $\mathcal{N}_{2}$ will take the form $O((2s-1)W^{2})$, $O((2s-1)^{2}W)$ and $O((2s-1)^{2}W^{3})$ such that the resulting equations remain consistent with the symmetry $(\rho,\psi,\phi,s)\to(\rho,-\psi,-\phi,\frac{1}{2}-s)$. These terms are directly responsible for the stability of hexagonal patterns Cox and Matthews (2003) as confirmed in the numerical simulations (see Fig. 5). We are now going to derive, in a pragmatic fashion, the amplitude equations for the fields when $s=1/2$. Our derivation is inspired by the methods presented in Winterbottom (2006).

We sit in the regime where patterns appear and we ask what the selected patterns beyond threshold are? Weakly nonlinear analysis begins by noticing that $\sigma_{3}\sim\varepsilon^{2}-a(k_{c}^{2}-k^{2})^{2}$, above the pattern threshold. We work in units of the the slow time scale by defining $T=\varepsilon^{2}t$; similarly in units of the large wavelength scale, we set $X=\varepsilon x$ which governs the evolution of the envelope of the fast growing patterns that develop at wavenumber $k_{c}$. The stationary homogeneous solution is perturbed when $\omega<\omega_{c}$. We expand the fields in a power series of the parameter $\varepsilon$. In the symmetric case $s=1/2$, the stable patterns are usually rolls and squares Cox and Matthews (2003). To study their relative stability in two dimensions, using $\psi_{h}=\phi_{h}=0$, the expansion for the fields reads

[TABLE]

with $Y=\varepsilon y$, and where $R(X,Y)$ is the large scale envelope of the marginal mode $k=0$ that has to be added in the expansion of the conserved field with the appropriate scaling to obtain a closure relation (see Coullet and Fauve (1985); Matthews and Cox (2000)). The functions $\rho_{n}$, $\psi_{n}$ and $\phi_{n}$ are expected to be products of slow dynamics envelopes and fast growing patterns. These considerations allow us to write differential operators with the chain rule, namely, $\partial_{x}\to\partial_{x}+\varepsilon\partial_{X}$, $\partial_{y}\to\partial_{y}+\varepsilon\partial_{Y}$ and $\partial_{t}\to\varepsilon^{2}\partial_{T}$. Next, we expand Eq. (19), (20), (21) to successive orders to get a closed set of equations. In the canonical case of the Swift-Hohenberg equation Swift and Hohenberg (1977); Cross and Hohenberg (1993), the closed relation for the lowest order amplitude is obtained to order $O(\varepsilon^{3})$. In our case of existence of a conserved quantity, we have to extract field evolution up to order $O(\varepsilon^{4})$ to get a closed system of equations. We are going to proceed recursively to extract the evolution of the fields. To order $O(\varepsilon)$, we find the following system:

[TABLE]

with

[TABLE]

and where $\nabla^{2}=\partial_{x}^{2}+\partial_{y}^{2}$. The solution of this system reads

[TABLE]

with $\lambda_{1}=B\phi_{0}/(\tilde{r}+k_{c}^{2})$ a simple scalar coefficient, $P_{1}$ and $Q_{1}$ are scalar functions of $X$ and $Y$, and where $c.c.$ stands for complex conjugate. To order $O(\varepsilon^{2})$, the system we arrive at is

[TABLE]

with $\zeta^{(2)}=(\zeta^{(2)}_{1},\zeta^{(2)}_{2},\zeta^{(2)}_{3})^{T}$ a vector which only depends on first order fields. The components of $\zeta^{(2)}$ read

[TABLE]

The solution of Eq. (54) reads:

[TABLE]

with $\lambda_{2}\equiv-B\lambda_{1}(\rho_{0}\lambda_{1}-\phi_{0})/2=\lambda_{1}^{2}(k_{c}^{2}+r)/2$ and where $Q_{1}^{*}$ is the complex conjugate of $Q_{1}$. At $O(\varepsilon^{3})$, we find the equation on $P_{1}$ (resp. $Q_{1}$) by collecting the terms proportional to $e^{ik_{c}x}$ (resp. $e^{ik_{c}y}$) in the two equations involving $\psi_{3}$ and $\phi_{3}$. A linear combination of these equations allows us to eliminate the second order amplitudes $P_{2}$ and $Q_{2}$ and to extract the time evolution on $P_{1}$ and $Q_{1}$. To order $O(\varepsilon^{3})$ we arrive at the following equations:

[TABLE]

with

[TABLE]

To order $O(\varepsilon^{4})$, we close the system with the time evolution of $R(X,Y)$, which is obtained by extracting coefficients of the mode $k=0$ in the $\rho_{4}$ equation. We obtain

[TABLE]

with

[TABLE]

We then perform a change of scale to fall back onto the canonical system found in Cox and Matthews (2003); Matthews and Cox (2000). Setting

[TABLE]

we define

[TABLE]

and we finally obtain

[TABLE]

with $\tilde{\nabla}\equiv\partial_{X}^{2}+\partial_{Y}^{2}$.

V.3 Roll and square stability

As we now deal with amplitude equations (79), (80) and (81) in a canonical form, the results obtained by Cox and Matthews (2003) are now directly transposable to our analysis. In particular, we can extract the stability boundaries of roll and square patterns, and predict modulational instabilities. The outcome of this analysis is that

• –

when $b_{2}+b_{3}>b_{1}$, rolls are unstable to one-dimensional disturbances (phase or amplitude modulation along the wave vector of patterns);

• –

if $b_{2}-b_{3}>b_{1}$ rolls undergo a two-dimensional instability, which is expressed through a transverse modulation of the rolls; see the dashed line in the phase diagram of Fig 11;

• –

if $g>1$, squares are unstable to rolls. Squares also undergo a modulational instability when $b_{2}>(1+g)b_{1}/2$ (we do not observe such patterns in the PDE solution); it turns out that in our model the condition for rolls to be unstable to squares is the same as the two-dimensional instability for rolls, and is thus described by the same dashed line in Fig 11.

Since we have $b_{3}<0$, the condition $b_{2}-b_{3}>b_{1}$ preempts $b_{2}+b_{3}>b_{1}$; it is shown in Cox and Matthews (2003) that the former then controls pattern formation. This explains why we observe the two-dimensional instability for rolls in Fig 12d. In our model, squares and transverse modulated rolls may exist separately at the same point of parameter space but they are ultimately selected by the geometry, the size, and the aspect ratio of the system (see magenta cross, Fig 11).

VI From equilibrium system to active system

So far, we have focused on the active system where spin flips are driven by a noise independent of temperature. Our analysis of the corresponding reaction-diffusion equations has shown the existence of a wealth of stationary patterns controlled by the values of the parameters of our model. These patterns simply do not exist in equilibrium when flips are controlled by temperature. To what extent does restoring a fraction of equilibrium spin flips within active flips suppresses the patterns we have obtained? Conversely, is adding a bit of activity over otherwise equilibrium flips sufficient to drive the system to a patterned stationary state? This section is about exploring the model system obtained by interpolating between fully active spin flips and equilibrium ones.

To implement both active and temperature controlled flips, the flipping rate $w(S_{k},\phi)$ for spin $S_{k}$ is now the sum of the active rate and the equilibrium rate which depends on the field value at the particle’s location $\phi_{k}\equiv\phi(\bm{r}_{k})$:

[TABLE]

and where $\eta$ is the equilibrium flipping rate if spins did not interact with the Gaussian field.

VI.1 Mean-field analysis

The mean-field equations are the same as Eq. (16) – (18) with $\alpha$ (resp. $\gamma$) changed into $w^{-}$ (resp. $w^{+}$):

[TABLE]

First, we search for a homogeneous stationary solution of this system. The self-consistent equation now verified by $\phi_{h}$ is more involved. We find that a homogeneous solution has

[TABLE]

where the function $G(\phi_{h})$ is given by

[TABLE]

We show a graphical solution of $G(\phi_{h})=\phi_{h}$ in Fig. 13.

We remark on Fig. 13(a) that the active fraction $s$ does not play a significant role when equilibrium flips are of the same order of magnitude as the active flip $\omega$. By contrast, when equilibrium flips are negligible, the homogeneous state is completely controlled by the fraction $s$ (see Fig. 13b). The most interesting regime is for $2\eta/\omega\sim 10^{-1}$ where the self-consistent equation has up to five solutions for some parameters (see Fig. 13c), unlike what we can observe in equilibrium (one or three solutions) or for the full active regime (one solution only). Finding out about the relative stability of these solutions comes first. This is the purpose of the following subsection.

VI.2 Linear stability analysis

We now perform a linear stability analysis of equations (83), (84) and (85), and study the stability for the different solutions of the self-consistent equation $\phi_{h}=G(\phi_{h})$. We restrict ourselves to the case $s=1/2$ which contains already rich physics. With this choice of $s$, we still have an up down symmetry for the spins, thus $\phi_{h}=0$ is always a homogeneous and stationary solution. We perform LSA around this solution and the results are shown on Fig. 14(a). The results are interesting when compared to the existence of other solutions of the self-consistent equations (see Fig. 14(b)). In particular, the homogeneous state develops patterns before other solutions for $\phi_{h}$ appear (case $\omega=0.45$ for instance). In addition, we observe in th explicit PDEs’ solutions that the final state of the system depends on initial conditions: even if the homogeneous state $\phi_{h}=0$ is not stable, the system may prefer creating patterns instead of having a full ferromagnetic order, which is also stable. To predict, at low cost, as a way of rationalizing our results, the final state of the system in this bi-stability regime, we can exhibit a mean field “free energy” $F$ whose minima are the possible homogeneous solutions for $\phi$ ($F$ is not a free energy since we are far from equilibrium). For a homogeneous density $\rho_{0}$, the evolution equation of the homogeneous field $\phi$ simply becomes

[TABLE]

with $F(\phi)=(r+B\rho_{0})\phi^{2}/2-\rho_{0}\ln\left(\omega+2\eta\cosh(B\phi_{0}\phi)\right)$ an even function of $\phi$ displayed on Fig. 14(c) for different values of $\omega$. The global minimum of the function $F$ corresponds indeed to the final state of the magnetization, namely $\langle\phi\rangle=0$ or $\langle\phi\rangle>0$, as observed in the PDE solution (see Fig. 15).

To sum up this section, we have seen that introducing a small amount of equilibrium in the dynamics of particle flip does not destroy the patterns. Furthermore, depending on its initial state, the system is now able to display either patterns, or ferromagnetic order, in striking contrast with a full active or equilibrium case where only one option is accessible. We now turn to the analysis of energy dissipation in the active system, and more precisely, we address the question of origin and of the location of entropy production.

VII Entropy production

Entropy production is a quantity that provides a measure of the degree of irreversibility of the dynamics of the system. In some cases it can simply be connected to the rate of energy dissipated by the system into the environment. We view entropy production as an elegant way to pinpoint the physical ingredients that are responsible for driving the system out of equilibrium. When spatially resolved, in the spirit of Nardini et al. (2017), entropy production may be used to connect the emerging structures, at a local level, with dissipation. While the nonequilibrium drive has been identified as the key ingredient for the generation of patterns, whether the genuinely nonequilibrium processes operate at the pattern boundaries, or within the bulk of the system, is a question of interest to us.

To estimate the total entropy production along a trajectory (in the whole phase space), we have to evaluate the probability of a trajectory relative to the probability of the time reversed trajectory Lebowitz and Spohn (1999). This question can be asked for various collections of degrees of freedom; we choose to focus on a single particle of position $X(t)\equiv X_{t}$ interacting with the Gaussian field $\phi(x,t)\equiv\phi_{t}(x)$. We restrict our derivation to a one dimensional system to keep the notation simple. The system evolves according to Eqs. (3) and (5) with $T=\Gamma=1$, and the spin flips from $+1$ to $-1$ (resp. from $-1$ to $+1$) with finite rate $\gamma$ (resp. $\alpha$). The spin $S_{t}$ jumps a finite number of times over an interval $[0,t_{F}]$, and $S_{t}$ is right continuous. The Hamiltonian $H$ will also be right continuous as a function of time. On an interval $[0,t_{F}]$ we define $t_{j}=jt_{F}/N$ and $t_{j+1}-t_{j}=t_{F}/N=\Delta t$. The probability of a noise history for the particle is given by

[TABLE]

Using the Itō convention, the probability of a trajectory $\{X(t)\}_{0<t<t_{F}}$ is equal to the probability of observing the corresponding noise history. We thus have

[TABLE]

where the set of points at which $H$ is discontinuous is of measure zero. At initial time $t_{I}=0$, we start with $X(0)=X_{I}$, $\phi(x,0)=\phi_{I}(x)$ and $S(0)=S_{I}$ and the system evolves to a final time $t_{F}$. We can define the time reversed noise history through $\tilde{\xi}(\tau)=\xi(t_{F}-\tau)$ and the reversed trajectory is then $\tilde{X}(\tau)=X(t_{F}-\tau)$ such that $\tilde{X}(0)=X(t_{F})\equiv\tilde{X}_{I}$. Entropy creation along a path is given by

[TABLE]

Similarly, since we have a Langevin equation for the field $\phi$, entropy creation of a field trajectory reads

[TABLE]

We also have entropy creation related to the realization of the sequence of flips. If we start from a down configuration and if we slice time into intervals of duration $\Delta t$ then entropy creation for a flip history writes

[TABLE]

where $\delta_{t_{F}}=-1$, [math] or $1$, depending on the initial and final values of the spin. We want to relate these previous results to the energy difference between the final time and the initial time. Let us start from the energy difference to see what terms appear in the calculation. We have

[TABLE]

and where the spin flips at time $t_{\alpha}$, with $t_{\alpha}^{+}$ and $t_{\alpha}^{-}$ denoting the times right after and right before the flip, respectively. We can now compute the entropy produced along any path:

[TABLE]

Dividing this result by $t_{F}$ and taking the limit $t_{F}\to\infty$ yields the entropy production rate $\sigma$. We immediately see that the first two terms in Eq. (104) vanish when divided by the total duration $t_{F}$ as $t_{F}$ is taken asymptotically large, since they are bounded. In the stationary state, the entropy production thus simplifies into

[TABLE]

with $N_{t_{F}}$ the number of flips in $[0,t_{F}]$. For a variable that flips between two states at fixed rates $\alpha$ and $\gamma$, this number is given by

[TABLE]

when $t_{F}\to\infty$, and thus scales like $O(t_{F})$. We can further simplify the expression of entropy production:

[TABLE]

The average $\langle S_{t_{\alpha}^{-}}\,\phi_{t_{\alpha}}(X_{t_{\alpha}})\rangle$ is however more complicated to compute because it depends on the whole dynamics. We are now considering two important limiting cases: (i) the time between two flips is large with respect to the particle–field dynamics, (ii) the time between two flips is small in that respect. In (i), we typically witness pattern formation. In this situation, the field at the particle’s location has the same sign as the spin before the flip, and thus scales like $O(\phi_{0}S_{t^{-}})$. We can further say that the field can equilibrate between flips, thus the field $\phi$ is equal to $\phi_{\mathrm{s.c.}}$ which satisfies the self-consistent equation $\phi_{\mathrm{s.c.}}=B\phi_{0}\rho_{0}\tanh(B\phi_{0}\phi_{\mathrm{s.c.}})/(r+B\rho_{0})$ (see Sec. IV.1) in the bulk of each microphase. In this case the entropy production rate reads:

[TABLE]

In regime (ii) of fast flipping, patterns disappear, and $S(t)$ and $\phi_{t}(X_{t})$ are almost uncorrelated. We can actually predict that the entropy production rate saturates when $\omega\to\infty$. We notice that for $N$ particles we have, by definition, $\psi(x,t)=\sum_{k=1}^{N}S_{k}(t)\delta(x-X_{k}(t))$ and thus for one particle $\psi(x,t)\sim S_{t}$. We thus approximate $\langle S_{t_{\alpha}^{-}}\,\phi_{t_{\alpha}}(X_{t_{\alpha}})\rangle$ with its continuum description, namely $\langle\psi_{t}\phi_{t}(X_{t})\rangle$. From the evolution equations of the fields, we compute $\phi\times\partial_{t}\psi+\psi\times\partial_{t}\phi$ (where $\partial_{t}\psi$ is given in Eq. (20) and $\partial_{t}\phi$ is given in Eq. (21)) to reconstruct a time derivative of a correlation, which is [math] in steady state. In the particular case of symmetric flips $s=1/2$, we have:

[TABLE]

which yields $\langle\psi\phi\rangle\underset{\omega\to\infty}{\to}B\phi_{0}\langle\psi^{2}\rangle/\omega$, from which we infer the produced entropy for $s=1/2$:

[TABLE]

In our numerical experiments, we start to measure entropy production at time [math] and we define the entropy production rate for a particle $k$ at time $t$ as

[TABLE]

where the $t_{\alpha k}$ are the time of flips of particle $k$. In Fig. 16 (left), we display the convergence of the entropy production rate towards its stationary state value for $100$ particles (out of $N=6479$ in the simulation), starting from a homogeneous state $\langle\phi\rangle=0$ and all particles spin up. Measuring entropy production for different values of $\omega$, we recover that entropy production rate saturates when $\omega\to\infty$, and we find the value we predicted in Eq. (110) for $\omega\to 0$. These results are displayed in Fig. 16 (right), where the graph also shows that the critical flipping rate predicted by the LSA matches the transition observed in entropy production.

Finally, another way to extract interesting information from our calculation for entropy production is to define a local entropy production rate or density of entropy production such that $\sigma=\int d^{2}r\sigma(\bm{r})$. Returning to a two-dimensional system, from Eq. (105), we can identify such an entropy production density for a system with $N$ particles:

[TABLE]

where the $(t_{\alpha k})_{\alpha\in\mathbb{N}}$ are the instants of flip of particle $k$, and where $\bm{X}^{k}_{t}$ is the position of particle $k$ at time $t$. We are now able to establish a map of the entropy production rate within the stationary state. In our simulations, though we observe diffusion of the whole pattern, the entropy production rate converges over a much smaller time scale, and we thus reach a “stationary state” before pattern blurring. In Fig. 17, we see that entropy production is localized within the bulk of the stripes. In other words, dissipation occurs in the bulk and not specifically at the boundaries of patterns. While the existence of patterns is a genuine nonequilibrium effect, one cannot interpret the role of the nonequilibrium drive in terms of a stabilizing effective surface tension at the boundaries of the ordered domains.

VIII Conclusion

Our goal was to explore and predict the emergence of collective phenomena in assemblies of active particles whose interactions are mediated by a fluctuating medium. We have done so on the basis of a minimal model in which particles diffuse while locally constraining the medium deformation. Activity is introduced by means of an internal degree of freedom that controls the interaction with the background field. This internal degree of freedom fluctuates independently of the bath temperature, and thus breaks the equilibrium nature of the dynamics of the whole system.

By means of Monte Carlo simulations and of a mean-field analysis of the dynamical equations, we have shown that this system displays a wealth of pattern formation regimes. When patterns appear, their wavelength is given by the geometric mean of the characteristic correlation length of the underlying elastic field and of the diffusion length of particles between two active flips. This geometric mean property is reminiscent of the typical wavelength emerging in crystal growth and in the Mullins-Sekerka instability Langer and Müller-Krumbhaar (1977). This coincidence might a posteriori be perceived as little surprising since we have at stake, in both systems, interactions favoring phase separation (a surface tension ingredient) competing with a diffusive process. In addition, as the number of particles is conserved in the system, patterns can be localized on a small fraction of the system size Matthews and Cox (2000); Cox and Matthews (2003). We have also examined how to interpolate between equilibrium dynamics and active dynamics for the flips since we reasonably expect that the flips might also feature temperature induced fluctuations. This interpolation has shown that the patterns could survive a moderate amount of equilibrium. Finally, we addressed the question of energy dissipation and entropy production in the active system. We have seen that entropy production vanishes for low flipping rates, as expected, and that it saturates for large flipping rates. We have also seen that entropy is produced within the bulk of the patterns, as opposed to other active systems where it is localized at the phase boundaries Nardini et al. (2017).

We are now at a stage where our model should be made more realistic. This may be achieved in a variety of directions. The Hamiltonian for the field can be adapted to specific systems we want to describe. Typically, we could use a Helfrich Hamiltonian to work on biological membranes. The field dynamics may also be changed. If the field now stands for a molecular density, we expect it to evolve according to a conserved dynamics (Cahn-Hilliard, Allen-Cahn). To focus on active proteins in the biological membrane, we also believe that hydrodynamic effects should be taken into account. This would certainly imply dealing with non-local equations, with the drag in a two-dimensional liquid layer (the lipid leaflet), and with the three-dimensional bulk liquid, which drives the system to another level of complexity, along with a (probably) richer behavior.

Acknowledgements.

R. Z. would like to thank Fridtjof Brauns for useful discussions. The authors are indebted to S. Fauve and J. Tailleur for their very useful feedback.