Spatial Organization of Active Particles with Field Mediated Interactions
Ruben Zakine Jean-Baptiste Fournier Frédéric van Wijland
cond-mat.soft cond-mat.stat-mech




 




Université de Paris, Laboratoire Matière et Systèmes Complexes (MSC), UMR 7057 CNRS, F-75205 Paris, France



...


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.




Spatial Organization of Active Particles with Field Mediated Interactions
    Ruben Zakine, Jean-Baptiste Fournier and Frédéric van Wijland
    December 30, 2023
=========================================================================





§ 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 <cit.>, liquid crystals <cit.>, capillary interfaces <cit.> and bio-membranes <cit.>.
Effective interactions, already at the level of equilibrium systems, may lead to a  rich phase behavior <cit.>. 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 <cit.>, or whose hydrodynamic interactions with the medium <cit.>, 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 <cit.>. 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 <cit.>, whose interactions have been experimentally explored <cit.> 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 <cit.> 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. <ref>. 
Its stationary phase diagram is explored in Sec. <ref>, 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. <ref> 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. <ref>, 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 <cit.>.
 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 <cit.>.
In our final two sections, we discuss the role of the nonequilibrium drive in producing patterned phases. To this end we introduce, in Sec.<ref>, 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. <ref>.





§ MODEL


The first ingredient of our model is a fluctuating field ϕ 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 0, and the field has a finite correlation length ξ. The Hamiltonian of the field then reads

    H_0=∫ d^2 x [r/2ϕ^2+ c/2 (∇ϕ)^2],

with ξ=√(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 r_k. The value ϕ( r_k) of the field at the position r_k of particle k is elastically constrained to the value ±ϕ_0 by the internal degree of freedom S_k=±1 of the particle. This leads us to use the following interaction energy between N particles and the field 

    H_int=∑_k=1^N B/2(ϕ( r_k)- S_kϕ_0)^2,

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 <cit.> should be used instead of Eq. (<ref>). A description of the membrane thickness with the Landau-Ginzburg Hamiltonian (<ref>) 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 <cit.>. 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:

    ∂_t ϕ( x,t) =-Γδ H/δϕ( x,t)+√(2Γ T) ζ( x,t),
       ⟨ζ( x,t)ζ( x',t')⟩ =δ( x- x')δ(t-t'),

where H=H_0+H_int, Γ is a mobility coefficient, T is the temperature in energy units and ζ( 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:

    d r_k/d t=-μ∂ H/∂ r_k+√(2μ T)ξ_k(t),
       ⟨ξ_k^α(t) ξ_ℓ^β(t')⟩=δ_αβδ_kℓδ(t-t'),

where μ is a mobility coefficient, and the ξ_k^α(t) are the  components of independent Gaussian white noises with zero average. We use the simplifying assumption that ξ_k and ζ 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 α and γ:

    S_k=-1 [ γ ] α  S_k=+1.

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 ϕ, we normalize energies by T, times by a^2/(Γ c) and we absorb c in a redefinition of the field ϕ. We thus carry out the replacements x/a→ x, Γ ct/a^2→ t, cϕ^2/T→ϕ^2, cϕ_0^2/T→ϕ_0^2, a^2r/c→ r, B/c→ B, Tμ/(Γ c)→μ, a^2α/(Γ c)→α and a^2γ/(Γ c)→γ. In a nutshell, rescaling time, space, fields and constants boils down to a=c=T=Γ=1. Our model being now defined, we present the results of our numerical Monte Carlo-based exploration of its properties.




§ NUMERICAL SIMULATIONS



We perform Monte Carlo dynamical simulations on a two dimensional square lattice of size L_x× L_y with periodic boundary conditions. The Gaussian field is defined on each site (i,j)  and takes continuous real values ϕ_ij. The field ϕ evolves according to the explicit stochastic Euler scheme corresponding to the dimensionless form of Eq. (<ref>):

    ϕ_ij(t+Δ t)=   ϕ_ij(t)-Δ t[rϕ_ij(t)-∇^2ϕ_ij(t)
       + B∑_k=1^N (ϕ_ r_k(t)-S_kϕ_0) ]+ G(0,2Δ t),

where the discrete Laplacian is defined as

    ∇^2ϕ_ij≡ϕ_i+1,j+1+ϕ_i-1,j+ϕ_i,j+1+ϕ_i,j-1-4ϕ_ij,

 r_k≡(i_k,j_k) being the position of the particle k on the lattice, and G(0,2Δ t) being a Gaussian variable of mean 0 and variance 2Δ 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 <cit.> 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'_k,j'_k) is given by 

    Δ H_ r_k→ r'_k=B/2(ϕ_i'_k,j'_k-ϕ_i_k,j_k)(ϕ_i'_k,j'_k+ϕ_i_k,j_k-2ϕ_0 S_k).

The following jump probability P_(i_k,j_k)→(i'_k,j'_k) between times t and t+Δ t implements a discrete version of the Langevin equation (<ref>):

    P_(i_k,j_k)→(i'_k,j'_k) =μΔ t exp(-Δ H_ r_k→ r'_k/2).

According to our model, the spin flip probability of particle k, P^f_k, is fixed by the rates α and γ, 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 Δ t small enough to ensure that the probabilities verify

    ∑_(i'_k,j'_k) P_(i_k,j_k)→(i'_k,j'_k)+P^f_k<1,

then  the probability P^n_k that particle k neither jumps nor flips is given by P^n_k=1-[∑_(i'_k,j'_k) P_(i_k,j_k)→(i'_k,j'_k)+P^f_k].
We take Δ t=2× 10^-5.





 §.§ 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+Δ t is

    P^f_k(↑→↓)=ηΔ t exp(-H_k↓-H_k↑/2),

where η is an inverse time-scale and where H_k↑ (resp. H_k↓) 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. <ref>).
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.<ref>c).
As a consistency check, we verified  that the equilibrium phase diagrams do not depend on dynamical parameters.
The specific simulations shown in Fig. <ref> were performed with μ=5, η=1.







 §.§ Active system





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

    P_k^f=γΔ t     S_k=+1
    αΔ t     S_k=-1.

If flipping rates are symmetric (α=γ) 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 ω=α+γ and the mean fraction s=α/ω of spin-up particles in steady state. The flipping rates are thus given by α=sω and γ=(1-s)ω.

In the following, we explore the phase diagram of the system for s=1/2 (Fig. <ref>, 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. <ref>a). At higher densities, the phenomenology becomes richer. Increasing B from the homogeneous phase, macroscopic stripes of both magnetization (Fig. <ref>b) are observed. 
As B is further increased, the stripes harbor the continuous nucleation of small lumps of particles of opposite magnetization (Fig. <ref>c). 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. <ref>d).
In the patterned phase (stripes or clusters), increasing the flipping frequency ω yields local mixing of the spins, which results in the homogeneization of the whole system (Fig. <ref>, top).









If now asymmetric flipping rates (s 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. <ref>).





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.






















§ MEAN-FIELD BEHAVIOR





 §.§ 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, ρ=ρ^++ ρ ^-, and a non conserved field ψ=ρ^+-ρ^-. From the Hamiltonian H, we write down a mean-field free energy density:

    f_MF=r/2ϕ^2   +Bρ+ψ/4(ϕ-ϕ_0)^2+Bρ-ψ/4(ϕ+ϕ_0)^2
         +ρ+ψ/2lnρ+ψ/2 +ρ-ψ/2lnρ-ψ/2,

where the first three terms are directly inferred from the energy functional H, while the last two ones reflect the particles' entropy.
Since ρ is the only conserved quantity, we minimize f_MF with respect to ϕ and ψ. We obtain ϕ=Bϕ_0ψ/(r+Bρ) and ψ=ρtanh(Bϕ_0ϕ). This imposes a self-consistent equation on ϕ=Bϕ_ 0ρtanh(Bϕ_0ϕ)/(r+Bρ) similar to what is obtained for the magnetization in the Ising model. Searching for homogeneous phases, yields either ϕ=0 (paramagnetic phase), or ϕ 0 (ferromagnetic phase). At low values of B, the system is uniform and there is a continuous paramagnetic–ferromagnetic transition at B_c^(MF)=(1+√(1+4r ϕ_0^2/ρ))/(2ϕ_0^2).
At higher values of B, we numerically solve the double tangent construction on f_MF(ρ) (already minimized with respect to ψ and ϕ). 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. <ref> 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 ϕ_0.




















 §.§ 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 <cit.> 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 ρ^±:

    ∂_tρ^+   =
    μ∇·[ρ^+∇∂ f_MF/∂ρ^+]+αρ^–γρ^+,
    
    ∂_tρ^-   =μ∇·[ρ^-∇∂ f_MF/∂ρ^-]
    -αρ^-+γρ^+
    
    ∂_tϕ   = ∇^2 ϕ - rϕ- B ρ^+(ϕ-ϕ_0)- B ρ^-(ϕ+ϕ_0).

It will prove convenient to write these equations in terms of the conserved field ρ, and of the non-conserved field ψ. We also parametrize the rates α and γ by means of ω=α+γ and s=α/ω (the latter being the steady-state fraction of spin up particles). The dynamical evolutions of the fields then read

    ∂_tρ   =
    μ∇^2ρ+μ B ∇·[(ρ ϕ-ψϕ_0)∇ϕ],
    ∂_tψ   =
    μ∇^2ψ+μ B ∇·[(ψϕ-ρϕ_0)∇ϕ]
               -ω ψ + (2s-1)ωρ,
    
    ∂_tϕ   = ∇^2 ϕ - rϕ- B ρϕ+ Bϕ_0ψ.

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=1/2 will play a special role because then these equations are invariant upon the up-down symmetry (ρ,ψ,ϕ)→(ρ,-ψ,-ϕ). Before we embark in a detailed analytical study of their pattern content, we begin with a numerical solution of these nonlinear coupled PDEs.






 §.§ 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× L_y=150× 150. The three fields ρ^+, ρ^- and ϕ 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:

    ρ_ij^+(t+Δ t)=ρ_ij^+(t)+   Δ t[ 
    μ∇^2ρ_ij^+
       +μ B ∇_x( ρ^+_ij(ϕ_ij-ϕ_0)∇_x ϕ_ij)
       +μ B ∇_y( ρ^+_ij(ϕ_ij-ϕ_0)∇_y ϕ_ij)
       -γρ^+_ij+αρ^-_ij](t),


    ϕ_ij(t+Δ t)=ϕ_ij(t)+   Δ t[∇^2ϕ_ij-rϕ_ij
       -Bρ^+_ij(ϕ_ij-ϕ_0)
       -Bρ_ij^-(ϕ_ij+ϕ_0)](t),

and the discretized equation on ρ^- is formally identical to the discretized equation on ρ^+ up to the exchange ρ^+↔ρ^-, ϕ_0→-ϕ_0, α↔γ.
The discrete spatial derivatives of any field g_ij are defined as

    ∇_x g_ij≡1/2(g_i+1,j-g_i-1,j),
       ∇_yg_ij≡1/2(g_i,j+1-g_i,j-1),

and the Laplacian has already been defined in Eq. (<ref>).
 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∑_ijρ_ij=ρ_0, along with the positivity of ρ^+ and ρ^- on each site.

 
To ease comparison of the PDE solution with the Monte Carlo simulation, the PDE phase diagram is plotted Fig. <ref> (bottom) for the same physical parameters as those of the Monte Carlo results of Fig. <ref> (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. <ref>: we can observe either homogeneous patterns (Fig. <ref>, top), or spatially localized patterns (Fig. <ref>, bottom), depending on the parameters. The coming section is devoted to an analysis of these patterns.










§ PATTERN ANALYSIS






 §.§ 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. (<ref>), (<ref>), (<ref>) is characterized by the following values of the fields

    ρ_h   =ρ_0,
    ψ_h   =(2s-1)ρ_0,
    ϕ_h   =(2s-1)Bρ_0/r̃ϕ_0,

with 

    r̃=r+Bρ_0=ξ̃^-2

where ξ̃ is the renormalized correlation length of the field ϕ. We set ρ_1=ρ-ρ_h, ψ_1=ψ-ψ_h and  ϕ_1=ϕ-ϕ_h and we expand Eqs. (<ref>), (<ref>), (<ref>) to linear order in the  ρ_1, ϕ_1, ψ_1 fields. We expand the fields in Fourier modes ∼ e^i k· x and we arrive at a linear system for the Fourier components ∂_t(ρ_1,ψ_1,ϕ_1)^T=M(ρ_1,ψ_1,ϕ_1)^T

with

    M=[                                                -μ k^2                                                     0                -μ Bρ_0ϕ_0(2s-1)(Bρ_0/r̃-1)k^2
    (2s-1)ω                                              -μ k^2-ω -μ Bρ_0ϕ_0 [(2s-1)^2Bρ_0/r̃-1]k^2
    -B^2ρ_0ϕ_0(2s-1)/r̃                                                  Bϕ_0                                               -k^2-r̃ ],


with k= k.
The eigenvalues of M can be shown to be always real which excludes oscillating patterns close to the threshold. We denote them  by σ_i, with σ_1<σ_2<σ_3. Solving M(k)=0 yields the modes for which temporal growth is marginal. In practice, we have M=-μ^2 k^2 Q(k^2), where Q(X)=X^2+q_1 X+ q_2 is degree 2 polynomial, with 

    q_1   =ω/μ+r̃-B^2 ρ_0 ϕ_0^2+B^3 ρ_0^2ϕ_0^2(2r̃-Bρ_0)(1-2s)^2/r̃^2,
    
    q_2   =ω/μ r̃^2(r̃^3-B^2ρ_0 r^2ϕ_0^2(1-2s)^2).

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







  * case (i): q_1≥0 and q_2≥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_+≡ 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. <ref>.  When ω is non zero, we sit in regime (ii) where q_2<0 is equivalent to 

    (2s-1)^2>1/Bϕ_0^2(1+Bρ_0/r)^2(1+r/Bρ_0),

which surprisingly does not depend upon the dynamical parameters. Physically, the instability comes from the frustrated field ϕ whose value at rest and without particles is 0, different from ϕ_h in the presence of particles with non-symmetric flipping rates. This regime is referred to as type II_s in <cit.>.

  * case (iii): q_1^2/4≥ q_2>0; patterns appear at finite wavelength (referred to as type I_s in <cit.>). We have X_-,X_+>0 and we set k_±≡ X_±^1/2. The eigenvalue σ_3(k) is positive for k∈[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 ω=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, ρ_0, ϕ_0 and s. By contrast, at fixed equilibrium parameters,  we transition from regime (i) to regime (iii) by changing ω or μ. 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. (<ref>), 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:

    σ_1   =-μ k^2,
    σ_2   =1/2(-r̃-ω-(1+μ)k^2-√(Λ)),
    σ_3   =1/2(-r̃-ω-(1+μ)k^2+√(Λ)),

with 

    Λ=[k^2(μ-1)+ω-r̃]^2+4μ k^2 ρ_0 B^2ϕ_0^2>0.

For the purpose of discussion we use ω=α+γ as the control parameter. Physically, we recall that for high flipping rates, the system remains homogeneous since particles locally efficiently mix, whereas for ω= 0, the system undergoes an equilibrium coarsening (see phase diagram Fig. <ref>).
Solving q_1^2-4q_2=0 yields a critical value of ω:

    ω_c=μ(Bϕ_0√(ρ_0)-√(r̃))^2,

below which the homogeneous system is no longer stable.
To study the system close to this transition, we write ω=ω_c-^2, where the distance to the threshold ^2>0 becomes our control parameter.
Since we sit in regime (iii), destabilization occurs at a mode k_c=k_±>0 when =0. Thus, when ω<ω_c, σ_3(k)⩾0 for k∈[k_-,k_+], with

    k_±=√(ρ_0 B^2 ϕ_0^2-r̃-Ω±√((ρ_0 B^2ϕ_0^2-r̃-Ω)^2-4r̃Ω))/√(2),

and where Ω=ω/μ.
The condition of existence of the k_± modes (namely that X_± are real) is given by

    Bϕ_0√(ρ_0)⩾√(r+Bρ_0)+√(ω/μ).

At ω=ω_c equality is achieved in Eq. (<ref>) and this allows us to infer the critical wavelength of patterns  λ_c≡ 2π/k_c= 2π(r̃ω_c/μ)^-1/4. This suggests that close to the threshold the patterns spatial periodicity λ_ p could be the combination λ_ p= 2π/k_ p=2π(r̃ω/μ)^-1/4.  

This prediction has been checked in simulations of a quasi 1D system of size 1000× 10 to force pattern formation along one direction, hence allowing us to achieve a good precision on the wavelength. We note on Fig. <ref> that the prediction on the pattern periodicity λ_ p applies beyond the pattern formation threshold.
Interestingly enough, λ_p  can be expressed as the geometric mean of the renormalized correlation length ξ̃ of the Gaussian field ϕ in presence of inclusions and of the typical diffusion length ℓ_d∼√(μ/ω) of a particle between two flips. The formula  λ_p∼(ξ̃ℓ_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 λ_ 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 ξ̃ of the field is t_d≈ξ̃^2/(2μ). On the other hand, the correlation time t_ϕ of the field over a scale ξ̃ is given by t_ϕ≈ 1/(2 r̃). Hence particles are fast with respect to the field when t_d≪ t_ϕ, or 1≪μ. In Fig. <ref>, 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.




 §.§ 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 <cit.>. 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 σ_3) along which we observe the temporal growth of the Fourier modes (see Fig. <ref>).

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

    P=[      1      0    0
    0   a(k) b(k)
    0      1      1 ],

with 

    a(k)=(1-μ)k^2+r̃-ω-√(Λ)/2Bϕ_0,
    
    b(k)=(1-μ)k^2+r̃-ω+√(Λ)/2Bϕ_0,

we define (U,V,W)^T=P^-1(ρ_1,ψ_1,ϕ_1)^T, with ρ_1, ψ_1 and ϕ_1 which are no longer infinitesimally small perturbations.
The new fields (U,V,W) now verify U∼ e^σ_1 t, V∼ e^σ_2 t and W∼ e^σ_3 t. 
However, we have to be careful: the existence of the conserved quantity ρ implies that the mode k=0 is marginal for the field U and has to be taken into account <cit.>. 
Indeed nonlinear terms couple modes k→ 0 with modes k∼ 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 σ_1, σ_2 and σ_3 around =0 and k∼ k_c. 
This would lead to neglect terms of the form V^qW^p and U^qW^p (for q⩾1 and p⩾0) since they exponentially go to 0 when k∼ 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

    ∂_t U   =σ_1 U+𝒩_1(U,W),
    ∂_t W   =σ_3 W+𝒩_2(U,W),

where 𝒩_1 and 𝒩_2 are nonlinear operators that couple U and W. To lowest order in , we find that 𝒩_1 contains terms ∼ O(W^2) and that 𝒩_2 contains terms ∼ O(UW). Thus, U will saturate to O(W^2), which renormalizes O(W^3) terms in 𝒩_2. Of course, this previous analysis holds for the case s=1/2, where equations are invariant upon the (ρ,ψ,ϕ)→(ρ,-ψ,-ϕ) symmetry. If s≠ 1/2, new terms appear in nonlinear equations (<ref>) and (<ref>). To lowest order, new terms in 𝒩_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 (ρ,ψ,ϕ,s)→(ρ,-ψ,-ϕ,1/2-s). These terms are directly responsible for the stability of hexagonal patterns <cit.> as confirmed in the numerical simulations (see Fig. <ref>). 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 <cit.>.


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 σ_3∼^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=^2 t; similarly in units of the large wavelength scale, we set X= 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 ω<ω_c. We expand the fields in a power series of the parameter . 
In the symmetric case s=1/2, the stable patterns are usually rolls and squares <cit.>.
To study their relative stability in two dimensions, using ψ_h=ϕ_h=0, the expansion for the fields reads

    ρ   =ρ_0+^2R(X,Y)+∑_n=1^∞^n ρ_n(x,y,X,Y),
    ψ   =∑_n=1^∞^n ψ_n(x,y,X,Y),
    ϕ   =∑_n=1^∞^n ϕ_n(x,y,X,Y),

with Y= 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 <cit.>). 
The functions ρ_n, ψ_n and ϕ_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, ∂_x→∂_x+∂_X, ∂_y→∂_y+∂_Y and ∂_t→^2∂_T. Next, we expand Eq. (<ref>), (<ref>), (<ref>) to successive orders to get a closed set of equations. In the canonical case of the Swift-Hohenberg equation <cit.>, the closed relation for the lowest order amplitude is obtained to order O(^3). In our case of existence of a conserved quantity, we have to extract field evolution up to order O(^4) to get a closed system of equations. We are going to proceed recursively to extract the evolution of the fields. To order O(), we find the following system:

    ℒ[ ρ_1; ψ_1; ϕ_1 ]=0,

with 

    ℒ=
    [          μ∇^2             0             0;             0     μ∇^2 -ω_c -μ Bρ_0ϕ_0∇^2;             0          Bϕ_0        ∇^2-r̃ ],

and where ∇^2=∂_x^2+∂_y^2. The solution of this system reads

    ρ_1(x,y,X,Y)=   0,
    ψ_1(x,y,X,Y)=   P_1(X,Y)e^i k_c x+Q_1(X,Y)e^i k_c y+c.c.,
    ϕ_1(x,y,X,Y)=   λ_1 ψ_1,

with λ_1=Bϕ_0/(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(^2), the system we arrive at is

    ℒ[ ρ_2; ψ_2; ϕ_2 ]=ζ^(2)(ρ_1,ψ_1,ϕ_1)

with ζ^(2)=(ζ^(2)_1,ζ^(2)_2,ζ^(2)_3)^T a vector which only depends on first order fields. The components of ζ^(2) read

    ζ^(2)_1=   μ B [ϕ_0(ψ_1∂_x^2ϕ_1+∂_xψ_1∂_xϕ_1)
       +ϕ_0(ψ_1∂_y^2ϕ_1+∂_yψ_1∂_yϕ_1)
       -ρ_0(∂_xϕ_1)^2-ρ_0ϕ_1∂_x^2ϕ_1
       -ρ_0(∂_yϕ_1)^2-ρ_0ϕ_1∂_y^2ϕ_1],
    ζ^(2)_2=   2μ[Bρ_0ϕ_0(∂_xX+∂_yY)ϕ_1-(∂_xX+∂_yY) ψ_1],
    ζ^(2)_3=   -2 (∂_xX+∂_yY) ϕ_1 .

The solution of Eq. (<ref>) reads:

    ρ_2=   λ_2 P_1^2(X,Y)e^2i k_c x 
    + λ_2 Q_1^2(X,Y)e^2i k_c y
       + 2λ_2 P_1(X,Y)Q_1(X,Y)e^i k_c x+ik_c y
       +2λ_2 P_1(X,Y)Q_1^*(X,Y)e^i k_c x-ik_c y+c.c.,
    
    ψ_2=    P_2(X,Y)e^i k_c x+Q_2(X,Y)e^i k_c y+c.c.,
    ϕ_2=   λ_1[P_2(X,Y)+2ik_cλ_1/r̃+k_c^2∂_X P_1(X,Y)]e^i k_c x
       +λ_1[Q_2(X,Y)+2ik_cλ_1/r̃+k_c^2∂_Y Q_1(X,Y)]e^i k_c y
       +c.c.,

with λ_2≡-Bλ_1(ρ_0λ_1-ϕ_0)/2=λ_1^2(k_c^2+r)/2 and where Q_1^* is the complex conjugate of Q_1. At O(^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 ψ_3 and ϕ_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(^3) we arrive at the following equations:

    ∂_T P_1=    a_1 P_1 + a_2 ∂_XXP_1
       - a_3 |P_1|^2 P_1-a_4 |Q_1|^2 P_1-a_5 R P_1 
    
    ∂_T Q_1=    a_1 Q_1 + a_2 ∂_YYQ_1
       - a_3 |Q_1|^2 Q_1-a_4 |P_1|^2 Q_1-a_5 R Q_1

with 

    a_1   =r̃/μ k_c^2+r̃,
    
    a_2   =4μ k_c^2 √(r̃)/B ϕ_0 (μ k_c^2 +r̃)√(ρ_0),
    
    a_3   =B^3 μ  k_c^2 ϕ_0^2 (B ϕ_0^2 ((k_c^2+r̃)^2-B^2 ρ_0^2)+2 (k_c^2+r̃)^2)/2 (k_c^2+r̃)^2 (μ  k_c^2 ρ_0B^2  ϕ_0^2+(k_c^2+r̃)^2)
    
    a_4   =2  μ  k_c^2 B^3ϕ_0^2 (B^2 ρ_0 ϕ_0^2 (k_c^2+ r)+(k_c^2+r̃)^2)/(k_c^2+r̃)^2 (μ  k_c^2 ρ_0B^2  ϕ_0^2 +(k_c^2+r̃)^2 ),
    
    a_5   =-μ  k_c^2 B^2ϕ_0^2 (k_c^2+ r)/μ  k_c^2 ρ_0 B^2  ϕ_0^2+(k_c^2+r̃)^2.

To order O(^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 ρ_4 equation. We obtain

    ∂_T R=   μ∇̃R+ κ_1(∂_X^2|P_1|^2+∂_Y^2|Q_1|^2)
       +κ_2(∂_Y^2|P_1|^2+∂_X^2|Q_1|^2),

with 

    κ_1   =μ B^2 ϕ_0^2(k_c^2-r)/(k_c^2+r̃)^2,
    κ_2   =-μ B^2 ϕ_0^2(k_c^2+r)/(k_c^2+r̃)^2.

We then perform a change of scale to fall back onto the canonical system found in <cit.>. Setting

    T→ T/a_1;
       X → X√(a_2/a_1)   ;   Y → Y√(a_2/a_1);
       P_1→ P_1 √(a_1/a_3)  ;   Q_1→ Q_1 √(a_1/a_3);
        R→ R  a_1/a_5,

we define

    g   =a_4/a_3>0,
    
    b_1   =μ/a_2>0,
    
    b_2   =(κ_1+κ_2)a_5/2a_2 a_3>0,
    
    b_3   =(κ_1-κ_2)a_5/2a_2 a_3<0,

and we finally obtain 

    ∂_T P_1=    P_1 + ∂_XXP_1
       -  |P_1|^2 P_1-g |Q_1|^2 P_1- R P_1 
    
    ∂_T Q_1=     Q_1 + ∂_YYQ_1
       -  |Q_1|^2 Q_1 - g|P_1|^2 Q_1- R Q_1
    
    ∂_T R=    b_1 ∇̃^2 R +b_2 ∇̃(|P_1|^2+|Q_1|^2)
       +b_3(∂_X^2-∂_Y^2)(|P_1|^2-|Q_1|^2),

with ∇̃≡∂_X^2+∂_Y^2.



 §.§ Roll and square stability





As we now deal with amplitude equations (<ref>), (<ref>) and (<ref>) in a canonical form, the results obtained by <cit.> 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 <ref>;

  – 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 <ref>.

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 <cit.> that the former then controls pattern formation. This explains why we observe the two-dimensional instability for rolls in Fig <ref>d. 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 <ref>).








§ 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, ϕ) 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 ϕ_k≡ϕ( r_k):

    w(S_k,ϕ) =   sω +η e^B ϕ_k ϕ_0≡ w_k^-     S_k=-1 
    
     (1-s)ω+η e^-B ϕ_k ϕ_0≡ w_k^+     S_k=+1 ,

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



 §.§ Mean-field analysis 

 
The mean-field equations are the same as Eq. (<ref>) – (<ref>) with α (resp. γ) changed into w^- (resp. w^+):

    ∂_tρ^+   =
    μ∇·[ρ^+∇∂ f_MF/∂ρ^+]+w^-ρ^–w^+ρ^+,
    
    ∂_tρ^-   =μ∇·[ρ^-∇∂ f_MF/∂ρ^-]
    -w^-ρ^-+w^+ρ^+,
    
    ∂_tϕ   = ∇^2 ϕ - rϕ- B ρ^+(ϕ-ϕ_0)- B ρ^-(ϕ+ϕ_0).

 First, we search for a homogeneous stationary solution of this system. The self-consistent equation now verified by ϕ_h is more involved. We find that a homogeneous solution has

    ρ_h   =ρ_0,
    ϕ_h   =G(ϕ_h)
    ψ_h   =r̃/Bϕ_0ϕ_h,

where the function G(ϕ_h) is given by

    G(ϕ_h)=Bρ_0/r̃(2s-1)ω+2ηsinh(Bϕ_hϕ_0)/ω+2ηcosh(Bϕ_hϕ_0)ϕ_0.

We show a graphical solution of G(ϕ_h)=ϕ_h in Fig. <ref>.

We remark on Fig. <ref>(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 ω. 
By contrast, when equilibrium flips are negligible, the homogeneous state is completely controlled by the fraction s (see Fig. <ref>b). 
The most interesting regime is for 2η/ω∼ 10^-1 where the self-consistent equation has up to five solutions for some parameters (see Fig. <ref>c), 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.










 §.§ Linear stability analysis



We now perform a linear stability analysis of equations (<ref>), (<ref>) and (<ref>), and study the stability for the different solutions of the self-consistent equation ϕ_h=G(ϕ_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 ϕ_h=0 is always a homogeneous and stationary solution. We perform LSA around this solution and the results are shown on Fig. <ref>(a). The results are interesting when compared to the existence of other solutions of the self-consistent equations (see Fig. <ref>(b)). In particular, the homogeneous state develops patterns before other solutions for ϕ_h appear (case ω=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 ϕ_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 ϕ (F is not a free energy since we are far from equilibrium). For a homogeneous density ρ_0, the evolution equation of the homogeneous field ϕ simply becomes

    ∂_tϕ=-∂ F/∂ϕ,

with
F(ϕ)=(r+Bρ_0)ϕ^2/2-ρ_0ln(ω+2ηcosh(Bϕ_0ϕ)   )
an even function of ϕ displayed on Fig. <ref>(c) for different values of ω. The global minimum of the function F corresponds indeed to the final state of the magnetization, namely ⟨ϕ⟩=0 or ⟨ϕ⟩>0, as observed in the PDE solution (see Fig. <ref>).

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.




§ 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 <cit.>, 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 <cit.>. This question can be asked for various collections of degrees of freedom; we choose to focus on a single particle of position X(t)≡ X_t interacting with the Gaussian field ϕ(x,t)≡ϕ_t(x). We restrict our derivation to a one dimensional system to keep the notation simple. The system evolves according to Eqs. (<ref>) and (<ref>) with T=Γ=1, and the spin flips from +1 to -1 (resp. from -1 to +1) with finite rate γ (resp. α). 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=Δ t. The probability of a noise history for the particle is given by 

    P[ξ|ξ_I]=   exp(-1/2∑_j=0^N-1Δξ_j^2 ).

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

    P[X|X_I]   =P[ξ|ξ_I]
    
    =   exp-1/4μ∑_j=0^N-1Δ t [ Ẋ_t_j+μ∂ H/∂ X_t_j[X_t_j,ϕ_t_j(X_t_j)]]^2
    ≃   exp-1/4μ∫_0^t_Fdτ[ Ẋ_τ+μ∂ H/∂ X_τ[X_τ,ϕ_τ(X_τ)]]^2,

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, ϕ(x,0)=ϕ_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 ξ̃(τ)=ξ(t_F-τ) and the reversed trajectory is then X̃(τ)=X(t_F-τ) such that X̃(0)=X(t_F)≡X̃_I. Entropy creation along a path is given by

    lnP[X(τ)|X_I]/P[X̃(τ)|X̃_I]
       =-1/4μ∫_0^t_F(Ẋ_τ+μ∂ H/∂ X_τ)^2-(Ẋ_τ-μ∂ H/∂ X_τ)^2 dτ
       =-∫_0^t_FẊ_τ∂ H/∂ X_τdτ.

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

    lnP[ϕ(x,τ)|ϕ_I(x)] / P[ϕ̃(x,τ)|ϕ̃_I(x)]
    =-∫_0^t_F dτ∫_x ∂_τϕ(x,τ)δ H/δϕ(x,τ).

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 Δ t then entropy creation for a flip history writes

    lnP[S(τ)|S_I]/P[S̃(τ)|S̃_I]=δ_t_Flnα/γ,

where δ_t_F=-1, 0 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

    H(t_F)-H(0)=   ∑_j=0^N-1  H[X_t_j+1,ϕ_t_j+1,S_t_j+1]
          -H[X_t_j,ϕ_t_j,S_t_j]
    
    
    =   ∑_j=0^N-1Δ t Ẋ_t_j∂ H/∂ X[X_t_j,ϕ_t_j,S_t_j+1]
         +Δ t ∫_xϕ̇_t_jδ H/δϕ(x)[X_t_j,ϕ_t_j,S_t_j+1]
         +H[X_t_j,ϕ_t_j,S_t_j+1]-H[X_t_j,ϕ_t_j,S_t_j]
         +o(Δ t)
    
    ≃   ∫_0^t_Fdτ[Ẋ_τ∂ H/∂ X_τ[X_τ,ϕ_τ,S_τ]
          +∫_xϕ̇_τδ H/δϕ(x)[X_τ,ϕ_τ,S_τ]]
       +∑_t_α[H[X_t_α^+,ϕ_t_α^+,S_t_α^+]
         -H[X_t_α^-,ϕ_t_α^-,S_t_α^-]],

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

    lnP[X_τ,ϕ_τ(x),S_τ|X_I,ϕ_I(x),S_I]/P[X̃_τ,ϕ̃_τ(x),S̃_τ|X̃_I, ϕ̃_I(x),S̃_I]
    
    =   -∫_0^t dτ(Ẋ_τ∂ H/∂ X+∫_x ∂_τϕ(x,τ) δ H/δϕ(x,τ))
       +δ_t_Flnα/γ
    
    =   -[H(t_F)-H(0)]+δ_t_Flnα/γ
       +∑_t_α[H[X_t_α^+,ϕ_t_α^+,S_t_α^+]-H[X_t_α^-,ϕ_t_α^-,S_t_α^-]]
    
    
    =   -[H(t_F)-H(0)]+δ_t_Flnα/γ
       +∑_t_α 2B ϕ_0 S_t_α^- ϕ_t_α(X_t_α).

Dividing this result by t_F and taking the limit t_F→∞ yields the entropy production rate σ. We immediately see that the first two terms in Eq. (<ref>) 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 

    σ   =lim_t_F→∞2B ϕ_0/t_F∑_0<t_α<t_F  S_t_α^- ϕ_t_α(X_t_α)
    
       =lim_t_F→∞2B ϕ_0/t_FN_t_F⟨ S_t_α^- ϕ_t_α(X_t_α)⟩,

with N_t_F the number of flips in [0,t_F]. For a variable that flips between two states at fixed rates α and γ, this number is given by

    N_t_F=2αγ/α+γ t_F,

when t_F→∞, and thus scales like O(t_F). We can further simplify the expression of entropy production:

    σ   =4B ϕ_0αγ/α+γ⟨ S_t_α^- ϕ_t_α(X_t_α)⟩
       =4 B ϕ_0 ω  s (1-s)⟨ S_t_α^- ϕ_t_α(X_t_α)⟩.

The average ⟨ S_t_α^- ϕ_t_α(X_t_α)⟩ 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(ϕ_0S_t^-). We can further say that the field can equilibrate between flips, thus the field ϕ is equal to ϕ_s.c. which satisfies the self-consistent equation ϕ_s.c.=Bϕ_0ρ_0 tanh(Bϕ_0ϕ_s.c.)/(r+Bρ_0) (see Sec. <ref>) in the bulk of each microphase. In this case the entropy production rate reads:

    σ_ω≪ω_c^max   =4 B ϕ_0ϕ_s.c.ω s(1-s).

In regime (ii) of fast flipping, patterns disappear, and S(t)  and  ϕ_t(X_t) are almost uncorrelated. We can actually predict that the entropy production rate saturates when ω→∞.  We notice that for N particles we have, by definition, ψ(x,t)= ∑_k=1^N S_k(t)δ(x-X_k(t)) and thus for one particle ψ(x,t)∼ S_t. We thus approximate ⟨ S_t_α^- ϕ_t_α(X_t_α)⟩ with its continuum description, namely ⟨ψ_t ϕ_t(X_t) ⟩. From the evolution equations of the fields, we compute ϕ×_tψ +ψ×_tϕ (where _tψ is given in Eq. (<ref>) and _tϕ is given in Eq. (<ref>)) to reconstruct a time derivative of a correlation, which is 0 in steady state. In the particular case of symmetric flips s=1/2, we have:

    0   =∂_t⟨ψϕ⟩
       =⟨ϕ(μ∇^2ψ+μ B ∇·[(ψϕ-ρϕ_0)∇ϕ]-ω ψ) 
              +ψ(∇^2 ϕ - rϕ- B ρϕ+ Bϕ_0ψ)⟩
       ω→∞=-ω⟨ψϕ⟩+Bϕ_0⟨ψ^2⟩,

which yields
⟨ψϕ⟩ω→∞→ Bϕ_0⟨ψ^2⟩/ω,
from which we infer the produced entropy for s=1/2:

    σ_ω→∞∝ B^2ϕ_0^2.

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

    σ_k(t)= 2B ϕ_0/t∑_0<t_α k<t  S^k_t_α k^- ϕ_t_α k(X^k_t_α k),

where the t_α k are the time of flips of particle k. In Fig. <ref> (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 ⟨ϕ⟩=0 and all particles spin up. Measuring entropy production for different values of ω, we recover that entropy production rate saturates when ω→∞,  and we find the value we predicted in Eq. (<ref>) for ω→ 0. These results are displayed in Fig. <ref> (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 σ=∫ d^2 rσ( r). Returning to a two-dimensional system, from Eq. (<ref>), we can identify such an entropy production density for a system with N particles:

    σ( r)=lim_t →∞1/t∑_k=1^N 2 Bϕ_0∑_t_α k<tS^k_t_α k^-ϕ_t_α k( r)δ( r- X^k_t_α k),

where the (t_α k)_α∈ℕ are the instants of flip of particle k, and where 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. <ref>, 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.









§ 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 <cit.>. 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 <cit.>. 
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 <cit.>.

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.




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.















































38
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Fisher and de Gennes(1978)]Fisher78
  author author M. E. Fisher and author P. G. de Gennes, @noop  journal journal C.
  R. Acad. Sci. Paris Ser. B volume 287, pages 207 (year 1978)NoStop
[Hertlein et al.(2008)Hertlein, Helden, Gambassi, Dietrich, and Bechinger]Hertlein08
  author author C. Hertlein, author L. Helden,
  author A. Gambassi, author S. Dietrich,  and author C. Bechinger, 10.1038/nature06443 journal journal Nature  volume 451, pages 172 (year
  2008)NoStop
[Ajdari et al.(1991)Ajdari,
  Peliti, and Prost]Ajdari91
  author author A. Ajdari, author L. Peliti,   and author J. Prost, 10.1103/PhysRevLett.66.1481 journal journal Phys. Rev. Lett. volume 66, pages 1481 (year 1991)NoStop
[Poulin et al.(1997)Poulin,
  Stark, Lubensky, and Weitz]poulin_1997
  author author P. Poulin, author H. Stark,
  author T. C. Lubensky,  and  author D. A. Weitz, 10.1126/science.275.5307.1770 journal journal Science volume 275, pages
  1770 (year 1997)NoStop
[Nicolson(1949)]Nicolson49
  author author M. M. Nicolson, @noop  journal journal
  Proc. Cambridge Philos. Soc. volume 45, pages 288 (year 1949)NoStop
[Hu and Bush(2005)]Hu05Nat
  author author D. L. Hu and author J. W. M. Bush, 10.1038/nature03995 journal
  journal Nature volume 437, pages 733 (year 2005)NoStop
[Goulian et al.(1993)Goulian, Bruinsma, and Pincus]Goulian93EPL
  author author M. Goulian, author R. Bruinsma,   and author P. Pincus, http://stacks.iop.org/0295-5075/22/i=2/a=012 journal
  journal EPL volume 22, pages 145 (year 1993)NoStop
[Dan et al.(1993)Dan,
  Pincus, and Safran]Dan93Lang
  author author N. Dan, author P. Pincus,  and  author S. A. Safran, 10.1021/la00035a005 journal journal
  Langmuir volume 9, pages 2768
  (year 1993)NoStop
[Dommersnes and Fournier(1999)]Dommersnes99EPJB
  author author P. Dommersnes and author J.-B. Fournier, 10.1007/s100510050968 journal
  journal Eur. Phys. J. B volume 12,  pages 9 (year 1999)NoStop
[Bitbol et al.(2012)Bitbol,
  Constantin, and Fournier]BitbolPlos12
  author author A.-F. Bitbol, author D. Constantin,
   and author J.-B. Fournier,  10.1371/journal.pone.0048306 journal
  journal PLoS ONE volume 7, pages e48306 (year 2012)NoStop
[Van Der Wel et al.(2016)Van
  Der Wel, Vahid, Šarić,
  Idema, Heinrich, and Kraft]vanderWel16
  author author C. Van
  Der Wel, author A. Vahid,
  author A. Šarić,
  author T. Idema, author D. Heinrich,  and author D. J. Kraft, 10.1038/srep32825 journal journal Scientific
  Reports volume 6 (year 2016),  10.1038/srep32825NoStop
[Brown et al.(2000)Brown,
  Smith, and Rennie]brown_2000
  author author A. B. D.  Brown, author C. G. Smith,  and author A. R. Rennie, 10.1103/PhysRevE.62.951 journal
  journal Phys. Rev. E volume 62,  pages 951 (year 2000)NoStop
[Noguchi and Fournier(2017)]noguchi_2017
  author author H. Noguchi and author J.-B. Fournier, 10.1039/C7SM00305F journal
  journal Soft Matter volume 13,  pages 4099 (year 2017)NoStop
[Miller and Bassler(2001)]miller_2001
  author author M. B. Miller and author B. L. Bassler, 10.1146/annurev.micro.55.1.165 journal journal Annual Review of Microbiology volume 55, pages 165 (year 2001),  note pMID: 11544353, http://arxiv.org/abs/https://doi.org/10.1146/annurev.micro.55.1.165
  https://doi.org/10.1146/annurev.micro.55.1.165 NoStop
[Stenhammar et al.(2017)Stenhammar, Nardini, Nash, Marenduzzo, and Morozov]stenhammar_2017
  author author J. Stenhammar, author C. Nardini,
  author R. W. Nash, author D. Marenduzzo,  and author A. Morozov, 10.1103/PhysRevLett.119.028005 journal journal
  Phys. Rev. Lett. volume 119, pages
  028005 (year 2017)NoStop
[Fribourg et al.(2014)Fribourg, Chami, Sorzano, Gubellini, Marabini, Marco, Jault, and Lévy]fribourg_2014
  author author P. F. Fribourg, author M. Chami,
  author C. O. S. Sorzano,
  author F. Gubellini, author R. Marabini, author
  S. Marco, author J.-M. Jault,  and author D. Lévy, 10.1016/j.jmb.2014.03.002
  journal journal Journal of Molecular Biology  volume 426, pages 2059 (year
  2014)NoStop
[Sumino et al.(2014)Sumino,
  Yamamoto, Iwamoto, Dewa,  and Oiki]sumino_2014
  author author A. Sumino, author D. Yamamoto,
  author M. Iwamoto, author T. Dewa,  and author
  S. Oiki, 10.1021/jz402491t journal journal The Journal of
  Physical Chemistry Letters volume 5, pages 578 (year 2014)NoStop
[Evans et al.(2013)Evans,
  Sun, Senyuk, Keller,
  Pergamenshchik, Lee, and Smalyukh]evans_2013
  author author J. S. Evans, author Y. Sun, author B. Senyuk, author
  P. Keller, author V. M. Pergamenshchik, author T. Lee,  and author I. I.  Smalyukh, 10.1103/PhysRevLett.110.187802
  journal journal Physical Review Letters  volume 110 (year 2013),  10.1103/PhysRevLett.110.187802NoStop
[Criante et al.(2013)Criante, Bracalente, Lucchetti,
  Simoni, and Brasselet]criante_2013
  author author L. Criante, author F. Bracalente,
  author L. Lucchetti, author F. Simoni,  and author
  E. Brasselet, 10.1039/c3sm50273b journal journal Soft Matter  volume 9, pages 5459 (year
  2013)NoStop
[Grawitter and Stark(2018)]grawitter_2018
  author author J. Grawitter and author H. Stark, 10.1039/C7SM02101A journal
  journal Soft Matter volume 14,  pages 1856 (year 2018)NoStop
[Chen(2004)]chen_2004
  author author H.-Y. Chen, 10.1103/PhysRevLett.92.168101 journal journal Physical Review Letters volume 92 (year 2004),  10.1103/PhysRevLett.92.168101NoStop
[Reigada et al.(2005a)Reigada, Buceta,  and Lindenberg]reigada_2005_1
  author author R. Reigada, author J. Buceta,   and author K. Lindenberg,  10.1103/PhysRevE.71.051906 journal journal Physical Review E volume 71 (year 2005a), 10.1103/PhysRevE.71.051906NoStop
[Reigada et al.(2005b)Reigada, Buceta,  and Lindenberg]reigada_2005_2
  author author R. Reigada, author J. Buceta,   and author K. Lindenberg,  10.1103/PhysRevE.72.051921 journal journal Physical Review E volume 72 (year 2005b), 10.1103/PhysRevE.72.051921NoStop
[Riecke(1992a)]riecke1992ginzburg_landau
  author author H. Riecke, https://doi.org/10.1016/0167-2789(92)90169-N
  journal journal Physica D: Nonlinear Phenomena  volume 61, pages 253  (year
  1992a)NoStop
[Matthews and Cox(2000)]matthews_2000
  author author P. C. Matthews and author S. M. Cox, 10.1088/0951-7715/13/4/317 journal
  journal Nonlinearity volume 13,  pages 1293 (year 2000), note arXiv:
  nlin/0006002NoStop
[Cox and Matthews(2003)]cox_2003
  author author S. Cox and author P. Matthews,  10.1016/S0167-2789(02)00733-9 journal
  journal Physica D: Nonlinear Phenomena volume 175, pages 196 (year
  2003)NoStop
[Winterbottom(2006)]winterbottom_2005
  author author D. M. Winterbottom, http://eprints.nottingham.ac.uk/10180/ title Pattern formation with a conservation law,  (year 2006)NoStop
[Helfrich(1973)]helfrich_1973
  author author W. Helfrich, @noop  journal journal Z.
  NaturForsch. volume C 28, pages 693
  (year 1973)NoStop
[West et al.(2009)West,
  Brown, and Schmid]west_2009
  author author B. West, author F. L. Brown,   and author F. Schmid, 10.1529/biophysj.108.138677 journal journal Biophysical Journal volume 96, pages 101 (year 2009)NoStop
[Krauth(2006)]krauth_2006
  author author W. Krauth, @noop  title Statistical Mechanics
  Algorithms and Computations (publisher Oxford University
  Press, address Oxford, year 2006)NoStop
[Dean(1996)]dean1996langevin
  author author D. S. Dean, @noop  journal journal Journal
  of Physics A: Mathematical and General volume 29,  pages L613 (year 1996)NoStop
[Cross and Hohenberg(1993)]cross_93
  author author M. C. Cross and author P. C. Hohenberg, 10.1103/RevModPhys.65.851 journal journal Reviews of Modern Physics volume 65, pages 851 (year
  1993)NoStop
[Swift and Hohenberg(1977)]swift_hohenberg_77
  author author J. Swift and author P. C. Hohenberg, 10.1103/PhysRevA.15.319 journal journal Physical Review A volume 15, pages 319 (year 1977)NoStop
[Coullet and Fauve(1985)]coullet_1985
  author author P. Coullet and author S. Fauve, 10.1103/PhysRevLett.55.2857 journal journal Physical Review Letters volume 55, pages 2857 (year
  1985)NoStop
[Riecke(1992b)]riecke_self-trapping_1992
  author author H. Riecke, 10.1103/PhysRevLett.68.301 journal journal Physical Review Letters volume 68, pages 301 (year
  1992b)NoStop
[Nardini et al.(2017)Nardini, Fodor, Tjhung, van
  Wijland, Tailleur, and Cates]nardini_2017
  author author C. Nardini, author E. Fodor,
  author E. Tjhung, author F. van Wijland, author
  J. Tailleur,  and author
  M. E. Cates, 10.1103/PhysRevX.7.021007 journal journal Phys.
  Rev. X volume 7, pages 021007
  (year 2017)NoStop
[Lebowitz and Spohn(1999)]Lebowitz-Spohn-1999
  author author J. L. Lebowitz and author H. Spohn, 10.1023/A:1004589714161 journal
  journal Journal of Statistical Physics volume 95, pages 333 (year 1999)NoStop
[Langer and Müller-Krumbhaar(1977)]langer_1977
  author author J. Langer and author J. Müller-Krumbhaar, 10.1016/0022-0248(77)90171-3
  journal journal Journal of Crystal Growth  volume 42, pages 11 (year
  1977)NoStop






