# Entropy production and its large deviations in an active lattice gas

Tal Agranov  DAMTP, Centre for Mathematical Sciences, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, United Kingdom  Michael E. Cates  DAMTP,
Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road,
Cambridge CB3 0WA, United Kingdom  Robert L. Jack  DAMTP, Centre for
Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge
CB3 0WA, United Kingdom  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom

###### Abstract

Active systems are characterized by a continuous production of entropy at
steady state. We study the statistics of entropy production within a lattice-
based model of interacting active particles that is capable of motility-
induced phase separation. Exploiting a recent formulation of the exact
fluctuating hydrodynamics for this model, we provide analytical results for
its entropy production statistics in both typical and atypical (biased)
regimes. This complements previous studies of the large deviation statistics
of entropy production in off-lattice active particle models that could only be
addressed numerically. Our analysis uncovers an unexpectedly intricate phase
diagram, with five different phases arising (under bias) within the parameter
regime where the unbiased system is in its homogeneous state. Notably, we find
the concurrence of first order and second order nonequilibrium phase
transition curves at a bias-induced tricritical point, a feature not yet
reported in previous studies of active systems.

##  I Introduction

The hallmark of active matter systems is the presence of a non-equilibrium
drive that acts at the single-particle level Marchetti  _et al._ (2013);
Bechinger  _et al._ (2016); Fodor and Cristina Marchetti (2018). This active
drive continuously dissipates heat into the environment, so that active
systems show a nonvanishing entropy production rate even at steady state.
Notably, these far-from-equilibrium steady states display a host of nontrivial
collective behavior that is not permitted in equilibrium systems. One of the
most striking and well studied of those is Motility Induced Phase Separation
(MIPS), where a system composed of self-propelled particles phase separates
into dense and dilute phases even when the inter-particle interactions are
purely repulsive Cates and Tailleur (2015).

It is now established that the entropy production rate (EPR) is a key
observable in the understanding of such collective phenomena and it has been
the subject of many recent studies Caprini  _et al._ (2019); Dabelow  _et al._
(2019); Ganguly and Chaudhuri (2013); Shankar and Marchetti (2018); Fodor  _et
al._ (2016); Mandal  _et al._ (2017); Li and Cates (2021); Caballero and Cates
(2020); Fodor  _et al._ (2022); O’Byrne  _et al._ (2022); Nardini  _et al._
(2017); Borthne  _et al._ (2020); Nardini  _et al._ (2017); Solon  _et al._
(2018a). This is partly because the EPR quantifies the actual level of
irreversibility in a system once the local rules for activity are combined
with interactions between particles. For instance if activity is modelled by
propulsive forces, but particles are jammed and unable to move, the active
forces do no work and the EPR vanishes. More generally it has proven useful
O’Byrne  _et al._ (2022); Fodor  _et al._ (2022) to define the informatic
entropy production rate (IEPR), defined by the log probability ratio of
forward and time reversed dynamics. In most settings, the dynamics of interest
is coarse-grained so that the IEPR quantifies the amount of irreversibility
that is detectable at coarse-grained level. Note that the IEPR is bounded
above by the microscopic EPR, with equality assured only if one avoids coarse-
graining entirely, and keeps track of all microscopic variables that are
present.

The IEPR enables us to quantitatively distinguish the emergent dynamics of
active systems from their passive counterparts. Its relation to collective
phenomena in active matter is two-fold. Firstly, measuring the IEPR was found
useful in identifying where activity comes to play in the collective behavior
of active systems, such as MIPS Fodor  _et al._ (2022). Secondly, a series of
recent studies found that tuning the value of IEPR in active systems, by
applying a dynamical bias, can cause phase transitions, leading to novel
states with a rich and interesting phenomenology, which may suggest new design
and control strategies for active matter Cagnetta  _et al._ (2017); Whitelam
_et al._ (2018); Tociu  _et al._ (2019); Gradenigo and Majumdar (2019); Nemoto
_et al._ (2019); Cagnetta and Mallmin (2020); Chiarantoni  _et al._ (2020);
Fodor  _et al._ (2020); Yan  _et al._ (2022); GrandPre  _et al._ (2021); Keta
_et al._ (2021). Such a bias allows one to explore sub-ensembles of dynamical
trajectories for which the IEPR is larger or smaller than usual. This
exploration is achieved by analysing the rate function for large deviations of
the IEPR; this can show singularities that are identified as non-equilibrium
phase transitions. The same biasing strategy can of course be used to look at
large deviations of other physical quantities such as currents. (This is done
for the specific model explored below in Agranov  _et al._ (2021) and Agranov
_et al._ (2022).)

Our present understanding of the nonequilibrium phase transitions that arise
on biasing IEPR is very far from complete. So far, these transitions of biased
active systems were uncovered primarily within numerical studies; there has
been no comprehensive analytical account. Numerically it is hard to resolve
the details (for instance the order of particular transitions) even for simple
models, and a complete mapping of the phase diagram is typically not possible.
Indeed, numerically studying large systems under a macroscopic bias,
especially near phase transitions, is nontrivial and requires advanced
sampling techniques Yan  _et al._ (2022). The general question of what to
expect for such a phase diagram therefore remained unanswered. Furthermore,
the relation of the phases and transitions of the biased system to those
appearing in the unbiased case such as MIPS, is still largely unknown.

For this reason, a fully-worked out analysis of large IEPR deviations for a
nontrivial model of active many-body dynamics is likely to be particularly
valuable. However, in almost all such models there is little prospect of
analytic progress: particles are strongly interacting and the system is far
from equilibrium, which hinders analytical characterisation of fluctuations.
Up until now, analytical progress was limited to either phenomenological
coarse-grained theories, or field-theoretical descriptions that employ mean-
field approximations and gradient expansions, see for example Peshkov  _et
al._ (2012); Bertin  _et al._ (2013); Wittkowski  _et al._ (2014); Nardini
_et al._ (2017); Solon  _et al._ (2018b, a); Tjhung  _et al._ (2018). While
these approaches proved most successful in accounting for the typical
behavior, they leave the full statistical description, and in particular, the
account of large macroscopic deviations which are at the center of the present
work, mostly unattended.

Another limitation of those approaches is that they do not admit a direct
comparison of the IEPR derived at the coarse grained or hydrodynamic scale
with the full EPR of the microscopic dynamics defined at particle level. To do
this one needs both quantities to be calculable within the same model. An
exact comparison of this type was so far limited to the case of noninteracting
particles only Pietzonka and Seifert (2018).

In this paper we provide the first analytically exact account of the IEPR
distribution in a well-chosen model of many interacting active particles of
the type capable of undergoing MIPS. We start by deriving exact expressions
for both the coarse-grained IEPR and the microscopic EPR evaluated over an
arbitrary system path. This offers the first exact comparison of the two in an
interacting active system. For reasons we explain, they are found to coincide,
creating an unexpectedly close link between the macroscopic IEPR and the
microscopic active work dissipated in this particular model Seifert (2012);
Pietzonka and Seifert (2018).

Next we turn to study the probability distribution of the IEPR at large times,
which obeys a large deviation principle and is described by a suitable rate
function. We provide a complete mapping of the dynamical phase diagram which
includes several previously reported phases Cagnetta  _et al._ (2017);
Whitelam  _et al._ (2018); Nemoto  _et al._ (2019); Chiarantoni  _et al._
(2020); Fodor  _et al._ (2020); Yan  _et al._ (2022); GrandPre  _et al._
(2021); Keta  _et al._ (2021). In particular, our analysis distinguishes two
kinds of stationary inhomogeneous states: phase separated states with sharp
interfaces, and states with smoothly modulated density. Arguments have been
presented for both such states Nemoto  _et al._ (2019); Fodor  _et al._
(2022), but previous analyses were not able to distinguish them. Our
analytically exact framework allows us to fully work out the intricate details
of the transitions to these phases, where we find a concurrence of first and
second order phase transition lines that meet at a nonequilibrium tricritical
point. Moreover, we uncover further novel phases that were not reported
before. We argue that these various phases should be generic in that they
should also arise if the statistics are biased with other observables of
interest, provided that these observables share the same symmetries as those
of the IEPR. We also establish the relation between the MIPS criticality in
the unbiased system and the transitions of the biased one. As was argued
before Nemoto  _et al._ (2019); Fodor  _et al._ (2022), when the system
approaches MIPS criticality, the onset of the transition to the smoothly
modulated phase happens at vanishingly small IEPR bias, so that the bias-
induced and motility-induced instabilities merge. A similar phenomena was
reported in the study of current fluctuations in Agranov  _et al._ (2022).

Our analytical breakthrough is possible thanks to recent progress in
establishing the exact fluctuating hydrodynamics of a lattice-based model for
repulsively interacting active particles in one dimension. For this active
lattice gas model, the noiseless coarse-grained hydrodynamics was established
first in Kourbane-Houssene  _et al._ (2018). Later, in Agranov  _et al._
(2021) and Agranov  _et al._ (2022) this was complemented with an exact large
deviation framework for macroscopic fluctuations. (We rely heavily on this
framework below.) The active lattice gas model shares many important features,
including MIPS Cates and Tailleur (2015), with previously studied active
particle models that are off-lattice and/or in higher dimensions. However, one
crucially important difference is that the microscopic rates are scaled with
system size. Consequently, in contrast to other standard models of active
particles, the particle magnetization density (the excess of right- over left-
moving particles) enters as a slow hydrodynamic field (a different scaling
would remove the slow magnetization). This slow field plays a major role in
allowing the analytic progress that we report below, essentially by ensuring
that the system remains mean-field-like, even in proximity to second-order
phase transitions. We believe that this is a small price to pay for obtaining
a set of exact results that can then serve as a benchmark for numerical
studies, and future approximate theories, involving models whose microscopic
motivation is more concrete and whose critical behavior is more complex. In
particular this choice of model allows an approximation-free analysis of the
phase diagram under bias (Fig. 4 below). This is already so complicated that
were it to emerge from a less exact theory, one would need to study its
robustness to the specific approximations made.

The rest of the paper is organized as follows: In Sec. II we review the active
lattice gas model, its corresponding fluctuating hydrodynamics description and
the associated MIPS. Then, in Sec. III we present the IEPR functional and its
expected value. In Sec. IV we present the large deviation framework that
accounts for the IEPR fluctuations and overview our main results, focusing for
simplicity on positive deviations. In Sec. V we derive the IEPR fluctuations
for the homogeneous phase. Sec. VI explains how this derivation helps expose
the other, inhomogeneous, phases whose derivation is detailed in Secs. VII and
VIII. Sec. IX concludes the mapping of the phase diagram and presents the
final expression for the large deviation function describing macroscopic IEPR
fluctuations. Here we also exploit a detailed fluctuation theorem that extends
all the results so far reported to negative IEPR values. Sec. X considers the
relation of the biased phase transitions to the MIPS criticality. Sec. XI
present a comparison to previous studies of IEPR in active systems. We
conclude and remark on future directions in Sec. XII.

##  II Active lattice gas model

###  II.1 Definition

The model we study was first introduced in Kourbane-Houssene  _et al._ (2018),
and we repeat its definition here. It consists of a one dimensional periodic
unit lattice of $L\gg 1$ sites where each site $i$ can either be occupied by a
$+$ particle ($\sigma_{i}^{+}=1$ and $\sigma_{i}^{-}=0$), a $-$ particle
($\sigma_{i}^{+}=0$ and $\sigma_{i}^{-}=1$), or be empty ($\sigma_{i}^{+}=0$
and $\sigma_{i}^{-}=0$). The number of particles in the system is
$N=\rho_{0}L$ where $\rho_{0}$ is the mean density.

The dynamics of the model is composed of three rules (see Fig. 1):

  1. (i)

Symmetric diffusion: A pair of neighboring sites exchange their state with
rate $D$.

  2. (ii)

Self propulsion: A $+$ ($-$) particle jumps to the right (left) neighboring
site with rate $\lambda/L$, provided the target site is empty.

  3. (iii)

Tumbling: A $+$ ($-$) particle converts into a $-$ ($+$) particle with rate
$\gamma/L^{2}$.

The scaling of the rates with $L$ ensures that in the hydrodynamic limit
($L\to\infty$ at fixed $\rho_{0}$) all processes occur on diffusive time
scales. Indeed, the time it takes for particles to traverse a macroscopic
system of size $L$ either via diffusive motion, or on account of self-
propulsion, scales as $L^{2}$ which is also the time scale for tumbling
events.

We identify two relevant length scales in this system, which are the
persistence length of active motion in the dilute limit,
$\hat{\ell}_{\mathrm{P}}=\lambda L/\gamma$, and the diffusive length
$\hat{\ell}_{\mathrm{D}}=L\sqrt{D/\gamma}$, which is the typical distance that
a particle diffuses, before it tumbles. One observes that both these length
scales are proportional to $L$ so they are macroscopic, due to the
$L$-dependence of the rates. These $L$-dependent rates may not be the obvious
choices for a physical model: one typically expects the tumbling time to be
faster, so that the length scales
$\hat{\ell}_{\mathrm{P}},\hat{\ell}_{\mathrm{D}}$ remain microscopic (of the
order of the particle size). However, the choice of $L$-dependent rates allow
for many exact results Kourbane-Houssene  _et al._ (2018) that would otherwise
not be achievable, including those described below. Interpretation of such
exact results can then be used to understand and interpret the results of more
realistic physical models.

.

Figure 1: Schematic representation of the microscopic dynamics of the active
lattice gas, as described in the text.

###  II.2 Hydrodynamic equations

To access hydrodynamic behaviour of this model, define a hydrodynamic co-
ordinate $x$ such that the position of site $i$ is $x=i/L$, which lies in the
unit interval $[0,1)$. The local density in the vicinity of this point is
obtained by a mesoscopic coarse-graining:

$\rho_{\pm}(x,t)=\frac{1}{2L^{\delta}}\sum_{|i-Lx|<L^{\delta}}\sigma_{i}^{\pm}\>,$
(1)

where the exponent $\delta$ can take any value in the range $0<\delta<1$. The
mean density is then

$\rho_{0}=\int_{0}^{1}\rho(x,t)dx\;.$ (2)

We write $\hat{t}$ for the time variable of the microscopic system: since the
system is diffusive, the corresponding hydrodynamic time variable is
$t=\hat{t}/L^{2}$. Define also the density and polarisation (or magnetisation)
fields:

$\displaystyle\rho$ $\displaystyle=\rho_{+}+\rho_{-}$ $\displaystyle m$
$\displaystyle=\rho_{+}-\rho_{-}\;.$ (3)

(The dependence of these fields on $(x,t)$ is left implicit, for a compact
notation.)

It was shown in Kourbane-Houssene  _et al._ (2018) that this model reaches
‘local equilibrium’ over mesoscopic time and length scales. By virtue of this
essential property, the density fields (1) obey hydrodynamic equations:

$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle-\partial_{x}J_{\rho}$ (4) $\displaystyle\partial_{t}m$
$\displaystyle=$ $\displaystyle-\partial_{x}J_{m}-2\gamma K,$ (5)

where $J_{\rho},J_{m}$ are the conservative fluxes and $K$ is the net rate of
change of the $\rho_{-}$ density due to tumbling reactions. (Note that $K$ is
used in other works to denote activity, which is not to be confused with the
net tumbling rate here.) The fluxes $J_{\rho,m}$ and the net tumbling rate $K$
fluctuate around their typical (most likely) values, which are:

$\displaystyle\bar{J}_{\rho}$ $\displaystyle=-D\partial_{x}\rho+\lambda
m(1-\rho)$ $\displaystyle\bar{J}_{m}$
$\displaystyle=-D\partial_{x}m+\lambda\rho(1-\rho)$ $\displaystyle\bar{K}$
$\displaystyle=m.$ (6)

(As above, the dependence of these mean fluxes and rates on the hydrodynamic
fields is left implicit, for compactness.)

The deterministic hydrodynamics given by these expressions was studied in
Kourbane-Houssene  _et al._ (2018) where a regime of MIPS was found. Steady
states of the model obey the time-independent equations

$\displaystyle 0$ $\displaystyle=$
$\displaystyle(\ell\partial_{x})\rho-\text{Pe}\,m(1-\rho)$ $\displaystyle 0$
$\displaystyle=$
$\displaystyle(\ell\partial_{x})^{2}m-\text{Pe}\,(\ell\partial_{x})\left[\rho(1-\rho)\right]-2m,$
(7)

where

$\displaystyle\ell$ $\displaystyle=\hat{\ell}_{\mathrm{D}}/L=\sqrt{D/\gamma}$
Pe $\displaystyle=\hat{\ell}_{\mathrm{P}}/\hat{\ell}_{\rm
D}=\lambda/\sqrt{D\gamma}\;.$ (8)

This $\ell$ is a rescaled version of the diffusive length, and the Péclet
number Pe is the ratio of the persistence length to the diffusive one. If Pe
and $\rho_{0}$ are sufficiently large, the hydrodynamic equations support non-
homogeneous steady states, corresponding to MIPS.

Two comments are in order. First, while this definition for the Péclet number
provides a measure for the self-propulsion strength, it is fundamentally
different from the common definition in studies of interacting active
particles, which compares the persistence length with the particle size. That
definition would be meaningless in this model, because the persistence length
${\hat{\ell}}_{\rm P}$ is macroscopic. Hence we choose the definition of Pe as
the ratio between the persistence length and the diffusive length, which are
both macroscopic.

Second, gradients in (7) always appear together with factors of $\ell$. The
result is that any inhomogeneous steady states will have spatial modulations
whose length scale are determined by this parameter (or equivalently by the
diffusive length $\hat{\ell}_{\rm D}$). In particular, MIPS appears in this
model Kourbane-Houssene  _et al._ (2018) as a coexistence of high- and low-
density phases, separated by narrow domain walls. In this work, we assume that
${\ell}\ll 1$ so that the interfacial width for MIPS is indeed narrow when
compared to the system size, although it is still macroscopic. (Anticipating
results below, note that biased ensembles can support phases with a smooth
density modulation: in these cases, narrow interfaces do not appear, even as
$\ell\to 0$.)

The binodal and spinodal curves for MIPS in this model have been calculated
(for $\ell\ll 1$) in Kourbane-Houssene  _et al._ (2018), and are shown in Fig.
2 for completeness. The critical point of the model
$(\rho_{c},\text{Pe}_{c})=(3/4,4)$ is located where the binodal and spinodal
curves meet.

.

Figure 2: The phase diagram for MIPS in the active lattice gas model. The
thick black curve is the binodal while the black dashed curve is the spinodal.
These meet at the critical point given by the red dot.

###  II.3 Hydrodynamic fluctuations

We now describe the fluctuating hydrodynamics of this model, as derived in
Agranov  _et al._ (2021) and Agranov  _et al._ (2022). For passive lattice
gasses, such a description was widely employed in the context of the
Macroscopic Fluctuation Theory (MFT), see Bertini  _et al._ (2015) for a
review. At the macroscopic level ($L\gg 1$), the system is overwhelmingly
likely to follow the deterministic hydrodynamic equations given above, but
significant deviations are still possible: they can be analysed by large
deviation theory. To obtain an explicit formula for the probability of these
rare events, we consider a trajectory (or history) of the system where the
hydrodynamic fields are measured, as are the corresponding fluxes: denote such
a trajectory as

${\cal
X}=\left\\{\rho(x,t),m(x,t),J_{\rho,m}(x,t),K(x,t)\right\\}_{x\in[0,1],t\in[0,T]}$
(9)

where $T$ is the trajectory duration (measured on the hydrodynamic scale). For
any trajectory, the measured fluxes and tumbling rates $J,K$ always obey
(4,5), but they differ in general from the deterministic values given in
(II.2). The probability of such a trajectory is then

$P({\cal X})\simeq{\rm e}^{-L{\cal A}({\cal X})}$ (10)

where the path action111Our notation for the action here is different from
previous works Agranov  _et al._ (2022) so as to avoid confusion with the sign
for the IEP in the following (17). is an integral of a Lagrangian density that
contains two terms:

$\mathcal{A}=\int_{0}^{T}dt\int_{0}^{1}\\!dx\,(\mathcal{L}_{J}+\mathcal{L}_{K}).$
(11)

The first contribution controls fluxes and is given by the Gaussian quadratic
form

$\mathcal{L}_{J}=\frac{1}{2}\begin{bmatrix}J_{\rho}-\bar{J}_{\rho}\\\
J_{m}-\bar{J}_{m}\end{bmatrix}^{\dagger}\mathbb{C}(\rho,m)^{-1}\begin{bmatrix}J_{\rho}-\bar{J}_{\rho}\\\
J_{m}-\bar{J}_{m}\end{bmatrix},$ (12)

with a symmetric correlation matrix

$\mathbb{C}(\rho,m)=D\begin{bmatrix}\sigma_{\rho}(\rho)&\sigma_{\rho,m}(\rho,m)\\\
\sigma_{\rho,m}(\rho,m)&\sigma_{m}(\rho,m)\end{bmatrix}$ (13)

where $\sigma_{\rho}(\rho)=2\rho(1-\rho)$, $\sigma_{m}(\rho,m)=2(\rho-m^{2})$,
and $\sigma_{\rho,m}(\rho,m)=2m(1-\rho)$. The contribution to the large
deviation functional controlling tumbling, $\mathcal{L}_{K}$, is non-Gaussian:

$\displaystyle\mathcal{L}_{K}=\gamma\left\\{\rho-\sqrt{K^{2}+(\rho^{2}-m^{2})}+K\ln\left[\frac{\sqrt{K^{2}+(\rho^{2}-m^{2})}+K}{(\rho+m)}\right]\right\\}.$
(14)

This form arises from the difference of the large deviation functionals for
the two Poisson processes ($\pm\to\mp$) that together determine the net rate
$K$; see Agranov  _et al._ (2022) for details.

We have neglected contributions to (11) from the initial condition of the
trajectory, these will be negligible in the following because we always
consider very long times $T$.

##  III entropy production functional

This work concentrates on fluctuations of the informatic entropy production in
the active lattice gas. For a given trajectory ${\cal X}$, the informatic
entropy production is defined Fodor  _et al._ (2022); O’Byrne  _et al._ (2022)
by comparing its probability with that of a corresponding time-reversed path
${\cal X}^{R}$, as222It has recently been emphasized that the proper
evaluation of informatic entropy production using a field theory requires a
careful definition of time discretization, say Itō or Stratonovich Cates  _et
al._ (2022). The scaling of fluctuation magnitude with $L^{-1/2}$ implies that
this choice is not important here.

$\Sigma({\cal X})=\ln[P({\cal X})/P({\cal X}^{R})].$ (15)

When defining ${\cal X}^{R}$, one has to specify the parity under time
reversal of the observable fields: different choices can reveal different
aspects of breaking of time reversal symmetry Fodor  _et al._ (2022). In the
present context, there is only a single choice that results in a nontrivial
expression for the entropy production: we must treat the fluxes and the
tumbling rate $K$ as odd under time reversal (this is natural since these are
velocity-like quantities). That is,

${\cal
X}^{R}=\left\\{\rho(x,T-t),m(x,T-t),-J_{\rho,m}(x,T-t),-K(x,T-t)\right\\}_{x\in[0,1],t\in[0,T]}$
(16)

which may be compared with (9). Notice that both $J$ and $K$ must be odd under
time reversal, to ensure consistency with the dynamical equations (4) and (5).

In non-equilibrium steady states, the entropy production $\Sigma(\cal X)$ is
proportional to the trajectory length $T$ so it is natural to define a time-
averaged IEPR as $\mathbb{S}_{T}({\cal X})=\Sigma({\cal X})/T$. Writing this
in terms of the action ${\cal A}$ and using the definition of ${\cal X}^{R}$
yields

$\mathbb{S}_{T}=\frac{L}{T}\int_{0}^{T}dt\int_{0}^{1}dx\left\\{2\begin{bmatrix}J_{\rho}\\\
J_{m}\end{bmatrix}^{\dagger}\mathbb{C}(\rho,m)^{-1}\begin{bmatrix}\bar{J}_{\rho}\\\
\bar{J}_{m}\end{bmatrix}+\gamma
K\ln\left(\frac{\rho+m}{\rho-m}\right)\right\\}.$ (17)

This expression can be significantly simplified by singling out a ‘non-
equilibrium’ contribution to (17). Indeed, notice that at vanishing bias
$\lambda=0$ the system becomes equilibrium. Correspondingly, we decompose the
expected currents (II.2) into their ‘equilibrium’ piece $-D\partial_{x}\rho$
and $-D\partial_{x}m$, and their non-equilibrium piece $\lambda m(1-\rho)$ and
$\lambda\rho(1-\rho)$ for the density and magnetization fluxes respectively.
For the equilibrium piece one finds that

$2\mathbb{C}(\rho,m)^{-1}\begin{bmatrix}-D\partial_{x}\rho\\\
-D\partial_{x}m\end{bmatrix}=\begin{bmatrix}-\partial_{x}\frac{\delta
F}{\delta\rho}\\\ -\partial_{x}\frac{\delta F}{\delta m}\end{bmatrix}$ (18)

with $F=L\int_{0}^{1}dx\mathcal{F}(\rho,m)$ being the free energy functional
of the equilibrium ($\lambda=0$) system with density

$\mathcal{F}(\rho,m)=\frac{\rho+m}{2}\ln\left(\frac{\rho+m}{2}\right)+\frac{\rho-m}{2}\ln\left(\frac{\rho-m}{2}\right)+(1-\rho)\ln\left(1-\rho\right).$
(19)

The expression Eq.(18) defines a vector of thermodynamic forces that drives
particles flow. The thermodynamic force that is conjugate to particles
flipping is derived from the same free energy

$\frac{1}{2}\ln\left(\frac{\rho+m}{\rho-m}\right)=\frac{\delta F}{\delta m}.$
(20)

Plugging these expressions into (17) we have

$\mathbb{S}_{T}=\frac{L}{T}\int_{0}^{T}dt\int_{0}^{1}dx\left\\{\begin{bmatrix}J_{\rho}\\\
J_{m}\end{bmatrix}^{\dagger}\begin{bmatrix}-\partial_{x}\frac{\delta
F}{\delta\rho}\\\ -\partial_{x}\frac{\delta F}{\delta m}\end{bmatrix}+2\gamma
K\frac{\delta F}{\delta
m}\right\\}+\frac{L\lambda}{T}\int_{0}^{T}dt\int_{0}^{1}dx\left\\{2\begin{bmatrix}J_{\rho}\\\
J_{m}\end{bmatrix}^{\dagger}\mathbb{C}(\rho,m)^{-1}\begin{bmatrix}m(1-\rho)\\\
\rho(1-\rho)\end{bmatrix}\right\\},$ (21)

Now using integration by parts and the dynamics (4) and (5) allows to
distinguish two contributions to the IEPR

$\mathbb{S}_{T}=-\frac{\Delta
F}{T}+\frac{L\lambda}{T}\int_{0}^{T}dt\int_{0}^{1}dx\left\\{2\begin{bmatrix}J_{\rho}\\\
J_{m}\end{bmatrix}^{\dagger}\mathbb{C}(\rho,m)^{-1}\begin{bmatrix}m(1-\rho)\\\
\rho(1-\rho)\end{bmatrix}\right\\},$ (22)

with $\Delta F=F(t=T)-F(t=0)$ the free energy difference between the final and
initial states of the system. As this difference is sub-extensive in time,
then at large times its contribution to (22) becomes negligible (for initial
and final conditions of bounded free energy $F$). From now on we will omit
this term. Finally, using the expression for $\mathbb{C}$ (13) one arrives at
the following simple expression for the IEPR

$\mathbb{S}_{T}=\frac{L\gamma\text{Pe}^{2}}{T}\int_{0}^{T}dt\int_{0}^{1}dx\alpha(x,t),$
(23)

where we defined the local rate of entropy production field

$\alpha(x,t)=\frac{1}{\lambda}\left[(1-\rho(x,t))J_{m}(x,t)+m(x,t)J_{\rho}(x,t)\right].$
(24)

The only contribution to this IEPR stems from the flux probability cost
$\mathcal{L}_{J}$, which alone encodes the irreversible, biased dynamics.
There is no explicit contribution from tumbling, but of course this reversible
process does influence the trajectories and hence the value of the IEPR.

It is interesting to compare this expression with the EPR for the underlying
microscopic model. A detailed analysis is given in Appendix A, which shows
that the macroscopic IEPR (24) coincides with the microscopic one. It is
unusual for the full microscopic EPR to be preserved under coarse-graining,
but this is a feature of the current model (which was, after all, constructed
to make the coarse-grained hydrodynamic equations more exact than usual). This
coincidence means that we can identify the coarse-grained IEPR in (23) with
the rate at which active work is expended at the particle level Seifert
(2012); Pietzonka and Seifert (2018).

###  III.1 The average IEPR

At macroscopic level, the mean IEPR can be computed from the most likely
trajectory. In homogeneous steady states, any such trajectory has
$\rho=\rho_{0}$, $J_{m}=\lambda\rho_{0}(1-\rho_{0})$ and $m=J_{\rho}=0$. The
mean IEPR is then

$\langle{\mathbb{S}}_{T}\rangle=L\gamma\text{Pe}^{2}\bar{a}(\rho_{0})$ (25)

where

$\bar{a}(\rho)=\rho(1-\rho)^{2}$ (26)

is a local rate of entropy production in a system at density $\rho$. To
interpret the non-trivial density dependence in (26), it is useful to consider
the microscopic dynamics.

In this model, entropy production arises from particles hopping to adjacent
vacant sites, as these are the only irreversible dynamical moves. The rate at
which the system performs such moves is related to the total current of
particles in their assigned ($\pm$) direction of motion, which is proportional
to $\rho_{0}(1-\rho_{0})$. Yet the total hopping rate into vacant sites is
only a fraction of the total current, as the latter also includes a
contribution from swaps between two occupied sites. With our chosen dynamical
rules, the latter are reversible and so do not contribute to the IEPR. As we
show in Appendix A, the relevant fraction of the total current is exactly
$(1-\rho_{0})$ which leads to the expression (26).

It will be important in the following that the function $\bar{a}(\rho_{0})$
loses convexity for $\rho_{0}<2/3$. This property will have far-reaching
consequences including a phase diagram of considerable complexity (arising, as
happens in equilibrium thermodynamics, from global minimization of a concave
thermodynamic potential). To the best of our knowledge, previous works, which
studied macroscopic fluctuations in related passive lattice gases, did not
consider observables with such a property, see e.g. Jack  _et al._ (2015);
Dolezal and Jack (2019); Lecomte  _et al._ (2012). In a forthcoming
publication Agranov  _et al._ we will further explore the phenomena related to
fluctuations of such observables in more general settings.

##  IV IEPR deviations: summary of main results

###  IV.1 Large deviations of the IEPR

Above we have defined the model and its IEPR. The central part of this work is
an analysis of large deviations of the IEPR. This section defines the problem
and gives an overview of the results. We work at the hydrodynamic level (which
is $L\gg 1$, with $L$ the number of lattice sites) and we further assume that
the trajectories are very long, $T\gg\gamma^{-1}$. Since $T$ is defined on the
hydrodynamic scale, this means that trajectories are long enough for typical
particles to tumble many times, and for a typical particle to explore the
whole system via diffusive hops.

For such long times, one expects a large-deviation form for the time-averaged
IEPR: that is, the the probability distribution for $\mathbb{S}_{T}$ behaves
as

$-\frac{1}{\gamma LT}\ln{\rm
Prob}\left(\frac{\mathbb{S}_{T}}{L\gamma\text{Pe}^{2}}\approx a\right)\simeq
I(a)$ (27)

where $I$ is the rate function. In this formula, the notation $\approx$
indicates that $\frac{\mathbb{S}_{T}}{L\gamma\text{Pe}^{2}}$ should lie in
some small interval around $a$, while $\simeq$ indicates asymptotic equality
in the joint limit of large $L$ and large $T$, see for example Touchette
(2009). This formula is valid when both $L$ and $T$ are large. For large $L$,
Eq. (10) already gives the probability of any individual hydrodynamic path, so
it is sufficient to minimise the action ${\cal A}$ over all paths that achieve
the desired IEPR, that is

$I(a)=\inf_{{\cal X}\colon\mathbb{S}_{T}({\cal
X})=aL\gamma\text{Pe}^{2}}\frac{\mathcal{A}({\cal X})}{\gamma T}\;,$ (28)

where the right hand side should be evaluated as $T\to\infty$. This means that
initial conditions can be neglected in the path probability formula (10). From
the definition of the IEPR it can be seen that if ${\cal X}$ solves the
minimisation problem for some $a>0$ then ${\cal X}^{R}$ solves the equivalent
problem for $a<0$, leading to a fluctuation theorem

$I(-a)=I(a)+\text{Pe}^{2}\,a$ (29)

see e.g. Seifert (2012); Jack (2020).

Solving this minimization problem fully describes the large deviations of
$\mathbb{S}_{T}$. The solutions will be computed in the following sections.
They depend on the mean density $\rho_{0}$ and on the parameters Pe and
$\ell$. As discussed in Sec. II.2, we shall focus on the behaviour for very
small $\ell$, for which we discuss the behaviour of $I(a)$ as a function of
$a,\text{Pe}$. As these variables are varied, the optimal trajectories can
change in character, leading to singularities in $I(a)$, which are interpreted
as dynamical phase transitions. We find five different dynamical phases in
this model, whose main characteristics are summarized below and in Table 1 and
Fig. 3.

  1. 1.

Homogeneous (H) phase. All fields are constant both in space and in time.
Moreover, the system is nonpolar with zero ‘magnetization’, $m=0$, and also
has vanishing mass flux $J_{\rho}=0$.

  2. 2.

Collective motion (CM) phase. All fields are still constant but the system
breaks reflection symmetry by developing a nonvanishing magnetization $m\neq
0$ and mass flux $J_{\rho}\neq 0$.

  3. 3.

Phase separated (PS) phase. The system breaks translational invariance with
the formation of sharply phase separate profiles. Reflection symmetry is
unbroken and there is no magnetization or current: $m,J_{\rho}=0$.

  4. 4.

Traveling band (TB) phase. Here a phase separated solution involves a
magnetized state $m,J_{\rho}\neq 0$ resulting in a sharply phase separated
traveling band solution. This state breaks translational and temporal
invariance and also reflection symmetry.

  5. 5.

Smoothly modulated (SM) phase. Here the system only breaks translational
invariance, but this time with a smoothly varying spatial profile.

Phase | Net current $J_{\rho}$ | Spatial structure | Temporal structure
---|---|---|---
Homogeneous, H | No | Homogeneous | Constant
Collective motion, CM | Yes | Homogeneous | Constant
Phase separated, PS | No | Sharp interface | Constant
Traveling band, TB | Yes | Sharp interface | Traveling wave
Smoothly modulated, SM | No | Smooth modulation | Constant
Table 1: The main characteristics of the different phases
---
Figure 3: States of the system at $\rho_{0}=0.75$ with increasing values of
$a$, the scaled IEPR. Here $\text{Pe}=2.5$ and $\ell=0.02$ in the panels of SM
and PS. The discontinuous profile in the panel TB corresponds to vanishing
$\ell=0^{+}$ where at finite $\ell$ the discontinuities are smoothed over a
scale $\ell$ as exemplified in the panel PS. The profile in the panel CM is
$\ell$-independent. The distinction between the two SM states is discussed in
Sec. VIII. The homogeneous phase that is found in the vicinity of $a=\bar{a}$
is omitted from the sequence.

Note that these phases (including those with non-zero magnetization) arise in
our model without any aligning interactions, for which typical steady-state
trajectories are homogeneous with zero mass flux. (The model also supports
MIPS steady states but we do not consider those parameter values in this work:
see Sec. XI for a discussion of that case.)

Fig. 3 shows the density profiles within the phases. The complete phase
diagram, as derived and explained in subsequent sections, is depicted in Fig.
4. The lines denote transitions between the different phases; these
transitions are accompanied by singularities of the rate function $I(a)$ as
detailed later in Fig. 11.

Figure 4: The nonequilibrium phase diagram for entropy production deviations
given by the minimization problem (28) at $\text{Pe}=2.5$ and vanishing
$\ell=0^{+}$ (A,B) and finite $\ell=0.05$ (C,D). (Panels B and D are
enlargements of the regions within the black dotted line in A and B
respectively.) The black solid line is the expected entropy production curve
(26). The black dots are tricritical points where a second order phase
transition line (solid orange) and first order (dashed orange) meet. The
miscibility gap that corresponds to the first order transition, not shown
here, is discussed in Sec. VIII, see also Fig 10. The critical lines in red,
orange, green, blue and purple are given by $a_{1-5}$ in Eqs. (36, 62, 43, 44)
and (45) respectively. The vertical blue segment is at $\rho_{0}=1/2$. The re-
entrant H phase arising for negative $a-\bar{a}$ at large density $\rho_{0}$
in (C,D) is stabilized by interfacial contributions and disappears as $\ell\to
0^{+}$ (A,B). Note that the first order (dashed orange) lines are associated
with miscibility gaps with tie-lines vertical (signifying time-like phase
separation). These gaps are too narrow to be shown here; see Fig. 10

###  IV.2 Finding the phases: overview

Since the derivation of these new phases and of the overall phase diagram is a
somewhat involved task, we give an overview of the subsequent calculations
here. The black line $a=\bar{a}(\rho)$ in the phase diagram (Fig. 4)
corresponds to the typical behaviour of the system; all other points
correspond to unusual values of the IEPR, and the corresponding optimal
trajectories are obtained by solving (28). Note that dynamical phase diagrams
are often plotted as a function of a biasing field, which parameterises the
distance of $a$ from $\bar{a}$ Fodor  _et al._ (2022); Jack (2020): in this
case the typical behaviour of the system would correspond to the zero-bias
axis of the phase diagram, instead of a non-trivial curve. Here, we plot the
phase diagram in terms of the IEPR itself. In the analogy with thermodynamics,
introducing the biasing field corresponds to a change of ensemble Touchette
(2009). In such phase diagrams, first-order phase transitions lead to regions
of phase coexistence in the phase diagram. In this work, we refer to such
regions as miscibility gaps, by analogy with the thermodynamics of fluid
mixtures. These miscibility gaps are very narrow in the system studied here so
they are not shown in Fig. 4 but they are included in Fig. 10, see Sec.
VIII.5.

Our derivation of the phase diagram starts with the H and CM phases in Sec. V.
In these cases all fields are constant both in space and time which means that
the minimization problem (28) becomes algebraic. We denote by
$I_{H}(\rho_{0},a)$ the rate function that corresponds to such homogeneous
solutions. In the H phase, the optimal trajectory in (28) has $m=0$ and
vanishing mass flux $J_{\rho}=0$. The resulting rate function is quadratic,
see Fig. 11 below. This quadratic form provides an upper bound for the rate
function everywhere in parameter space. The same upper bound was employed in a
series of recent works Pietzonka  _et al._ (2016); Horowitz and Gingrich
(2017) that established a non-equilibrium thermodynamic uncertainty relation
for current fluctuations. In this model, the bound is saturated only in the H
phase. In the other four phases, the minimisation problem is solved by
different kinds of trajectory, so $I(a)<I_{H}(\rho_{0},a)$, and dynamical
phase transitions occur.

The CM phase occurs when $a$ is larger than its typical value: here the system
spontaneously breaks the reflection symmetry, leading to a non-zero
magnetization $m\neq 0$ and mass flux $J_{\rho}\neq 0$. In Sec. V we identify
the mechanism behind this transition by looking at the microscopic dynamics. A
similar phase transition was reported in previous numerical studies of scalar
active systems Nemoto  _et al._ (2019); Fodor  _et al._ (2020); GrandPre  _et
al._ (2021); Keta  _et al._ (2021). Our analysis confirms that, within the
active lattice gas, this transition is described by a mean-field Landau
theory.

Inhomogeneous phases are discussed in Sec. VI. Analytical progress is possible
here by extending an analogy with equilibrium statistical mechanics that was
recently established in Agranov  _et al._ (2022) for the study of current
fluctuations in this model. Within this analogy, $I_{H}(\rho,a)$ takes the
role of the bulk term in a free energy functional of two conserved order
parameters. These are the spatially varying local entropy production rate
$a(x)$ and density $\rho(x)$ fields. As we show, the small length scale
$\ell\ll 1$ controls the analogue of the interfacial energy. Consequently, the
state of the system can be described by an analogy with the classical
equilibrium theory for phase separation where the rate function $I$ (27) is
given by a convex hull (common tangent plane) construction over $I_{H}$. This
construction is valid when the inhomogeneity involves sharp interfaces whose
width is controlled by $\ell$ and which contribute subdominantly relative to
the ‘bulk’ term $I_{H}$.

This leads to sharply phase-separated (PS) solutions that are found via a
common tangent construction over $I_{H}$, as discussed in Sec. VII. In this
phase the sharply separated profile does not break reflection symmetry
($m,J_{\rho}=0$) and is stationary. This phase was previously reported for
related systems when biased towards suppressed EPR values GrandPre  _et al._
(2021); Yan  _et al._ (2022); Fodor  _et al._ (2020); Chiarantoni  _et al._
(2020); Nemoto  _et al._ (2019); Cagnetta  _et al._ (2017). We show here that
this transition can also occur at enhanced IEPR biasing. We identify the
source of this phenomenon as the non-convexity of $I_{H}$ at enhanced IEPR
deviations for this model. We also report a second sharply phase-separated
state, for which the common tangent construction brings in a magnetized state
with $m,J_{\rho}\neq 0$, resulting in phase-separated profile that moves with
a constant velocity. This traveling band (TB) phase was not reported in
previous studies of other models, yet we expect it to be generic in the sense
that it it should appear whenever non-convexity of $I_{H}$ occurs alongside
reflection symmetry breaking ($m\neq 0$). Indeed, this phase was reported in a
recent study of the same model under bias for a different dynamical quantity
Agranov  _et al._ (2022).

Lastly we analyze the SM phase in Sec. VIII. Here the bulk term $I_{H}$ fails
to account for the state of the system due to a degeneracy in the common
tangent plane construction. This occurs because $I_{H}$ has the same value for
all unbiased states. (This could conceivably occur in the analagous problem of
a free energy function for a equilibrium system with two conserved order
parameters, but we have never seen an instance of it.) As a result, the state
of the biased system is determined by the sub-leading surface-tension-like
term in the rate functional. Solving the spatial minimization problem for this
functional, we find that the optimal solution is smoothly varying in space
rather than comprising bulk phases separated by sharp interfaces as in the PS
and TB states.

The SM phase was anticipated in a similar active context in Nemoto  _et al._
(2019) based on generic arguments. Here, by establishing its existence in the
active lattice gas model, we provide the first conclusive evidence for its
existence. Moreover, we find that the transition to this phase can be either
first or second order, depending on the value of the averaged density
$\rho_{0}$, see Fig. 4. Similar phase transitions were extensively studied in
the context of passive lattice gases Bodineau and Lagouge (2010); Lecomte  _et
al._ (2012); Jack  _et al._ (2015); Baek  _et al._ (2017, 2018); Dolezal and
Jack (2019), but for each given model were found to be of fixed order (first
or second). Here we find that the H$\to$SM transition occurs along both a
second order and a first order transition line, with these meeting at a
tricritical point. We trace back this phenomenon to the change of convexity in
(26).

In the next four sections we derive all the above-mentioned results. The
fluctuation theorem (29) means that it is sufficient to consider only $a\geq
0$ since the rate function with $a<0$ can be then be deduced, see Sec. IX.

##  V Homogeneous phases

We first consider solutions to the minimization problem (28) that are
homogeneous in space and time:

$\rho(x,t)=\rho_{0}\quad,\quad m(x,t)=m\quad,\quad
J_{\rho}(x,t)=J_{\rho}\quad,\quad J_{m}(x,t)=J_{m}.$ (30)

These will serve as a baseline for the following analysis of non-homogeneous
phases. Under these conditions, the minimization problem (28) is algebraic. In
the context of passive lattice gases the stationarity hypothesis is known as
the “additivity principle” Bodineau and Derrida (2004) and was shown to hold
in a variety of scenarios Bodineau and Derrida (2004); Hurtado and Garrido
(2009, 2010); Hurtado  _et al._ (2011); Prados  _et al._ (2011); Meerson  _et
al._ (2014); Meerson (2015); Agranov  _et al._ (2016); Agranov and Meerson
(2017, 2018); Agranov  _et al._ (2019). Its breaking signals a dynamical phase
transition, which has attracted significant interest in the literature, see
e.g. Bertini  _et al._ (2005); Bodineau and Derrida (2005); Bertini  _et al._
(2006); Shpielberg and Akkermans (2016); Zarfaty and Meerson (2016). As we
find below, however, a phase transition is possible in the active lattice gas
model even among simple time- and space-independent solutions.

We now turn to minimize (28) subject to (30). Since $\rho_{0}$ is a parameter,
the quantities to be optimized are $m,J_{\rho},J_{m},K$. The dynamical
equation (5) requires the net tumbling rate $K=0$. Second, constraining the
IEPR in (23) to obey $\mathbb{S}_{T}=aL\gamma\text{Pe}^{2}$ requires

$(1-\rho_{0})J_{m}+J_{\rho}m=a\lambda.$ (31)

Plugging these expressions into (11) and minimizing with respect to $J_{\rho}$
provides an additional algebraic relation between the fluxes, and one finally
arrives at

$J_{\rho}=\frac{ma\lambda}{\rho_{0}(1-\rho_{0})+m^{2}}\quad,\quad
J_{m}=\frac{\rho_{0}a\lambda}{\rho_{0}(1-\rho_{0})+m^{2}}.$ (32)

With these relations, the minimization problem (28) reduces to a single
minimization over the magnetization $m$:

$I_{H}(\rho_{0},a;\text{Pe})=\min_{m}\tilde{I}(m),$ (33)

with

$\tilde{I}(m)=\rho_{0}-\sqrt{\rho_{0}^{2}-m^{2}}+\text{Pe}^{2}\,\frac{\left\\{a-(1-\rho_{0})\left[\rho_{0}(1-\rho_{0})+m^{2}\right]\right\\}^{2}}{4(1-\rho_{0})\left[\rho_{0}(1-\rho_{0})+m^{2}\right]}$
(34)

Anticipating a second-order phase transition to the CM phase, we expand in
powers of $m$ to obtain

$\tilde{I}(m)=\text{Pe}^{2}\frac{(a-\bar{a})^{2}}{4\bar{a}}-\text{Pe}^{2}\,(1-\rho_{0})\frac{a^{2}-a_{1}(\rho_{0},\text{Pe})^{2}}{4\bar{a}^{2}}m^{2}+\frac{\bar{a}^{2}+2\text{Pe}^{2}\,\rho_{0}^{2}a^{2}}{8\rho_{0}^{3}\bar{a}^{2}}m^{4}+\mathcal{O}(m^{6})$
(35)

where $\bar{a}=\bar{a}(\rho_{0})$ is given by (26) and

$a_{1}(\rho_{0},\text{Pe})=\bar{a}(\rho_{0})\sqrt{\frac{2+\text{Pe}^{2}\,\rho_{0}(1-\rho_{0})}{\text{Pe}^{2}\,\rho_{0}(1-\rho_{0})}}.$
(36)

Interpreting Eq. (35) as a Landau free energy, one sees that $m$ acquires a
non-zero value whenever $a$ exceeds the critical value
$a_{1}(\rho_{0},\text{Pe})$. Clearly
$a_{1}(\rho_{0},\text{Pe})>\bar{a}(\rho_{0})$, so this symmetry breaking
transition occurs when biasing the system towards increased IEPR. This is a
second order phase transition: the optimal value $m$ (and the corresponding
mass flux $J_{\rho}$ (32)) grows continuously from zero on increasing $a$
through $a_{1}$. It follows that near the critical point

$m,J_{\rho}\sim\left[a-a_{1}(\rho_{0},\text{Pe})\right]^{1/2}$ (37)

which is the usual critical exponent for such transitions.

Inserting these results into (33), one obtains

$\displaystyle I_{H}(\rho_{0},a)=$
$\displaystyle\frac{\text{Pe}^{2}}{4}\frac{\left[a-\bar{a}(\rho_{0})\right]^{2}}{\bar{a}(\rho_{0})},$
$a\leq a_{1}(\rho_{0},\text{Pe})$,  $\displaystyle I_{H}(\rho_{0},a)=$
$\displaystyle\rho_{0}-\sqrt{\rho_{0}^{2}-m}+\frac{\text{Pe}^{2}}{4}\frac{\left[a-\bar{a}(\rho_{0})-(1-\rho_{0})m^{2}\right]^{2}}{\bar{a}(\rho_{0})+(1-\rho_{0})m^{2}},$
$a>a_{1}(\rho_{0},\text{Pe})$, (38)

where the relation $m=m(a;\rho_{0},\text{Pe})$ is given implicitly by

$a=\left[\bar{a}(\rho_{0})+(1-\rho_{0})m^{2}\right]\sqrt{1+\frac{2}{\text{Pe}^{2}}\frac{1}{(1-\rho_{0})\sqrt{\rho_{0}^{2}-m^{2}}}}.$
(39)

This $I_{H}$ is the variational estimate for the rate function, under the
assumption that the deviations in the IEPR occur via homogeneous trajectories.

As usual in the Landau-like scenario, the expression (38) for the minimized
rate function (analogous to a free energy) has a discontinuous second
derivative at $a=a_{1}$. Note that the quadratic form in (38), valid for
$a<a_{1}$, coincides with a universal bound on current fluctuations Pietzonka
_et al._ (2016); Horowitz and Gingrich (2017). This bound is saturated in the
H phase, but not in the CM (collective motion) phase, which is stationary and
homogeneous, but has nonzero magnetization $m$.

The mechanism behind this transition can be illuminated by looking at the
microscopic picture. As mentioned above in Sec. III and detailed in Appendix
A, the entropy production arises from the total particle current into adjacent
vacant sites. This current is only a fraction of the total current which,
acccording to the dynamical rules in the chosen model, includes a contribution
also from reversible swaps between adjacent sites (Fig. 1). By increasing the
imbalance $m$ between right-moving ($+$) and left-moving ($-$) particles, the
part of the particle current $J_{\rho}$ carried by reversible swaps is
reduced. At the extreme $|m|=\rho_{0}$, where all particles are of the same
type, reversible swaps contribute nothing to the current, maximizing the rate
of entropy production alongside the current itself.

This mechanism is closely related to the one behind a similar transition that
was reported in a system of active Brownian particles with repulsive pairwise
interactions Nemoto  _et al._ (2019). That model undergoes a transition to a
collective motion phase, when considering sufficiently large values of the
active work (which is effectively the microscopic EPR). In this CM state,
particles align their orientation and move collectively in a common direction.
This alignment decreases the rate of collisions between the particles, thereby
enhancing the rate at which particles perform irreversible motion (through
self-propulsion), and increasing the EPR. In the model we study here, particle
collisions are less effective, since particles can always swap positions by
diffusive dynamical moves (Fig. 1). Nevertheless, encounters between
oppositely oriented particles still decrease the rate of irreversible moves,
similar to the active Brownian particles. Thus, by aligning the particles’
orientation and reducing the rate of those encounters, the IEPR is enhanced.

##  VI Concavities and inhomogeneous phases

###  VI.1 Squared gradient approximation (small $\ell$)

The expression (38) in the previous Section was obtained by optimizing (28)
under the assumption of homogeneous trajectories. As such, it is a valid bound
on the rate function $I$, but there are cases where the global minimiser of
(28) is an inhomogeneous trajectory. Studying the expression (38) over the
entire $(\rho_{0},a)$ parameter space will allow us to understand why the
homogeneous trajectory becomes sub-optimal in some regions, and expose the
inhomogeneous phases that take over.

We first assume that the system is stationary with zero mass flux
$J_{\rho}=0$, so the fields vary only in space, not time (this is the
additivity principle Bodineau and Derrida (2004)). One can make progress in
the ensuing minimization problem by exploiting the small length scale $\ell$.
This length scale controls the gradient terms appearing in the action
functional (11). We show in Appendix B that the magnetization $m$ can be
optimized analytically, yielding a two-field minimization problem over the
density $\rho(x)$ and the flux $J_{m}(x)$. It is convenient to reparameterise
$J_{m}$ in terms of the local rate of entropy production (24)
$\alpha(x)=J_{m}(x)(1-\rho(x))/\lambda$.

At leading order in a gradient expansion, the minimization problem (28)
becomes

$I(a)\simeq\min_{\alpha(x),\rho(x)}\int_{0}^{1}dx\left\\{I_{H}(\rho(x),\alpha(x))+\ell^{2}\begin{bmatrix}\partial_{x}\rho\\\
\partial_{x}\alpha\end{bmatrix}^{\dagger}\mathbb{A}\begin{bmatrix}\partial_{x}\rho\\\
\partial_{x}\alpha\end{bmatrix}\right\\},$ (40)

subject to the constraints

$a=\int_{0}^{1}\alpha(x)dx\quad,\quad\rho_{0}=\int_{0}^{1}\rho(x)dx.$ (41)

with the matrix $\mathbb{A}(\rho,\alpha)$ given in Eq.(B.2.1) of Appendix B.

The bulk term $I_{H}$ could have been anticipated by the following argument.
Neglecting the gradient terms in the action functional altogether, we have a
minimization problem that is locally algebraic and coincides with that of the
previous section V. This means that the local contribution to the action
functional should match that for homogeneous profiles (38), yet now these
profiles can vary in space. Meanwhile the derivation of the sub-leading term
in (40) necessitates retaining the gradient terms in the functional (11) and
is derived in Appendix B.

###  VI.2 Optimization in states with sharp interfaces: analogy with
thermodynamics

To perform the minimization in (40), we again take inspiration from
equilibrium statistical mechanics: $I_{H}(\rho,\alpha)$ can be interpreted as
an analogue of a bulk free energy density as a function of two order
parameters $\alpha$ and $\rho$ that are spatially conserved due to (41). This
is analogous to globally minimizing the free energy of a ternary fluid
mixture, which is a well-studied problem in chemical physics (see, e.g.,
Andelman  _et al._ (1987)). So long as the interfacial energy term in (40) is
negligible to leading order, then whenever $I_{H}$ becomes non-convex, the
homogeneous solution becomes sub-optimal across a region that encloses the
non-convexity (but generally extends beyond it into locally convex regions).
The optimal solution then comprises a sharply phase-separated state whose
coexisting phases have $(\rho_{0},a)$ values found by a common tangent
construction over $I_{H}$. The resulting sum of ‘bulk free energies’
represents the lower convex envelope of $I_{H}$. (In standard thermodynamic
terminology, this envelope is referred to as the convex hull of $I_{H}$, we
adopt this convention throughout the following.)

For two conserved densities, as here, there can generically be either two
phases or three phases in coexistence, but the latter does not arise in our
model. Within any region of two-phase coexistence, each globally unstable
density point $(\rho_{0},a)$ lies on a unique ‘tie-line’, which is one of a
family of straight lines each connecting two stable points on the perimeter of
that region. Upon global minimization, the uniform state at $(\rho_{0},a)$
phase-splits into the two phases that are end-points of its tie line, with
phase volumes determined by (41). The interfacial energy in such a phase-
separated state scales as $\ell$ and is negligible for $\ell\ll 1$. Notice
however, that for such a sharply phase separated solution, the gradient
expansion (40) breaks down in the region of the sharp interface where higher
order gradient terms contribute. Nonetheless, the interfacial contribution is
sub-extensive so that the minimized rate function $I$ consists, to leading
order in $\ell$, of the convex hull of $I_{H}$. Wherever $I_{H}$ does not
coincide with its convex hull, inhomogeneous trajectories are preferred to
homogeneous ones in (28).

Figure 5: (A) The function $I_{H}$ (38) at $\text{Pe}=3.5$. The black curve
marks the minimum $I_{H}=0$ obtained for
$a=\bar{a}(\rho_{0})=\rho_{0}(1-\rho_{0})^{2}$. This is the locus of unbiased
states each of which is homogeneous. The red curve marks the critical line
(36) between H and CM phases. It is denoted by a dashed line within the region
where $I_{H}$ becomes non-convex, and the optimal solution is no longer
homogeneous. (B) and the blow up (C) show the regions of non-convexity of
$I_{H}$, see (38). Local non-convexity
$\text{Det}\left[\text{Hess}(I_{H})\right]<0$ in orange, and global non-
convexity in purple. See the Appendix C for explicit expressions and details
of their derivation. The purple curve is given by (45) and denotes the
boundary of the TB phase.

The problem is thereby reduced to finding the required convex hull. Notably,
$I_{H}(\rho_{0},a)$, is everywhere convex when viewed as a function of the
single variable $a$. However, viewed as a two dimensional function over the
parameter space $(\rho_{0},a)$ it is not everywhere convex. Fig. 5(A) shows a
plot of $I_{H}$ (38). To analyse the convexity of this function, we separate
two ways in which it can be non-convex. For some values of $(\rho_{0},a)$, the
function $I_{H}$ is locally non-convex (its Hessian has at least one negative
eigenvalue). For other values, the Hessian may be positive definite but
$I_{H}$ does not coincide with its convex hull: we refer to such regions as
globally non-convex. All these regions can be found analytically, see Appendix
C; they are plotted in Fig. 5(B,C).

Although these arguments assumed stationarity, they can also be easily adapted
to the travelling band (TB) phase in which $J_{\rho}\neq 0$ and a phase-
separated state translates at fixed speed through the system. As a result, the
leading order term that controls the rate function $I$ is controlled solely by
$I_{H}$ in the TB phase as well, with the same construction of the convex hull
via the common tangent plane. In this case, the tangent construction then
describes the two phases which comprise the traveling band state. This is
proven in Appendix B. Thus both the PS and TB regions are accessible by this
quasi-thermodynamic approach to the minimization in (40), essentially because
in each case the minimiser constists of two phases separated by a narrow
domain wall whose contribution to $I_{H}$ is sub-leading at small $\ell$. The
quantitative results of this procedure are detailed in Section VII below.

###  VI.3 The problem of SM states

The analysis so far has considered inhomogeneous states with sharp interfaces.
However, the system studied here also supports SM (smoothly modulated) phases,
where this approach fails. This is attributable to a degeneracy that is
generic in the large deviation context, although it would be quite exceptional
in a thermodynamic setting. Namely, the function $I_{H}(\rho_{0},a)$ in (38)
attains its minimum $I_{H}=0$ everywhere along the curve
$a=\bar{a}(\rho_{0})$, see Fig. 5. (In the thermodynamic context, one expects
instead to find minima of the bulk free energy at one or more isolated
points.) The curve where $I_{H}=0$ represents the set of unbiased states, one
for each $\rho_{0}$, all of which have zero rate function. As a result, at any
point in the two-dimensional parameter space $(\rho_{0},a)$ that is engulfed
by the convex hull of the curve $a=\bar{a}(\rho_{0})$, the rate function drops
onto the plane $I=0$ (modulo interfacial contributions). The convex hull of
the curve $a=\bar{a}(\rho_{0})$, refers to the convex hull of the two-
dimensional domain $0<a<\bar{a}(\rho_{0})$, in the $(\rho_{0},a)$ plane: this
is the green shaded region in Fig. 4(A,B). It is given by

$a<\begin{cases}\bar{a}(\rho_{0}),&\rho_{0}\leq(1/2)\\\
a_{3}(\rho_{0})&\rho_{0}>(1/2)\end{cases}$ (42)

with the line

$a_{3}(\rho_{0})=\frac{(1-\rho_{0})}{4}$ (43)

which is tangent to the curve $\bar{a}$ at the point of contact at
$\rho_{0}=1/2$. This makes the tangent construction degenerate, as any point
on that plane is now shared by infinitely many possible tie lines, each
representing a different pair of bulk phases in coexistence. More importantly,
in all of these phases the bulk contribution $I_{H}$ vanishes, so that the
interfacial tension terms now dominate.

To find the optimal solution in the presence of the degeneracy we therefore
retain the interfacial tension term in (40). We show in Section VIII that the
solution to this minimization problem is given by a smoothly modulated
solution that is not characterized by a sharp interface, making the minimized
rate function scale as $\mathcal{O}(\ell^{2})$.

##  VII Results for states with sharp interfaces: PS and TB

Having set out the general principles for solving the minimization problem
(28) in the presence of inhomogeneous states, we now analyse the PS and the TB
states in detail. (Recall that these states contain sharp interfaces between
bulk phases. They differ qualitatively in that the TB phase has a nonvanishing
mass flux and is nonstationary.) Both states are found by the common tangent
construction over $I_{H}$. Note that both these states are restricted to
densities above half filling, $\rho_{0}\geq 1/2$. The resulting high-density
phase has $\rho(x)=\rho_{h}$ and the low-density one has $\rho(x)=\rho_{l}$;
the corresponding local EPRs are $a_{h},a_{l}$.

Figure 6: A blow up of the ‘pinch point’ singularity at $(\rho_{0}=1,a=0)$ of
the homogeneous rate function (38)

It is explained in Appendix C that the common tangent construction is
dominated by the point $(\rho_{0}=1,a=0)$ where the function $I_{H}$ has a
‘pinch point’ singularity: $I_{H}\sim a^{2}/(1-\rho_{0})^{2}$, see Fig. 6.
Indeed, at close packing $\rho_{0}=1$, the rate function diverges for any non
vanishing IEPR $a\neq 0$, as all accessible dynamical moves are now
reversible. At the same time, the rate function vanishes at $a=0$ since this
is the expected value $\bar{a}$ when $\rho_{0}=1$. This point forms part of
the convex hull of $I_{H}$ and in fact represents the high-density phase for
all coexisting states. (Note that the usual requirement for local convexity at
a point of common tangency is bypassed because this point lies at the very
edge of the allowed $(\rho_{0},a)$ domain; ‘tangency’ here does not constrain
the derivatives.) Thus the high-density phase in both the PS and TB phases is
always the state of close packing which has vanishing IEPR
$(\rho_{h}=1,a_{h}=0)$.

On the other hand, the low-density phase $(\rho_{l},a_{l})$ is given by
different expressions in the PS and TB cases, as follows. For the phase
separated (PS) phase, the tangent construction involves the parabolic branch
of the rate function $I_{H}$ (38). The resulting low density phase is given by
half filling ($\rho_{l}=1/2$) but with an IEPR that varies with the global
density: $a_{l}=a/2(1-\rho_{0})$.

Finding the low-density phase in the travelling band (TB) phase is more
involved. The main result (Appendix C) is that the low-density state develops
a spontaneous magnetisation when $a>a_{4}(\rho_{0},\text{Pe})$ with

$a_{4}(\rho_{0},\text{Pe})=\frac{1}{4}(1-\rho_{0})\sqrt{1+\frac{8}{\text{Pe}^{2}}}$
(44)

The resulting symmetry-broken state is a TB: it corresponds to a low-density
CM region, while the high-density region remains symmetric ($m=0$) with
$\rho_{h}=1$. The properties of the low-density phase are
$(\rho_{l},m_{l},J_{\rho,l},J_{m,l})$; they can be computed by constructing
the relevant tie-line on the $(\rho_{0},a)$ plane as we detail in Appendix
C.2. Since the high density phase is invariant, this is simpler than in most
equilibrium problems. We first construct a line connecting the high density
phase at $(\rho_{h}=1,a_{h}=0)$ to the point of interest $(\rho_{0},a)$ set by
the chosen global mean values for the density and IEPR. Now extend this line
until it intersects the critical curve

$a_{5}(\rho_{0},\text{Pe})=\frac{\rho_{0}+\sqrt{2\rho_{0}-1}}{2}\sqrt{(1-\rho_{0})^{2}+\frac{2\sqrt{2}(1-\rho_{0})}{\text{Pe}^{2}\,\sqrt{\rho_{0}-\sqrt{2\rho_{0}-1}}}},$
(45)

which marks the boundary between the TB and CM phases. This curve is derived
in Appendix C.2. The intersection point marks the low-density phase,
$(\rho_{l},a_{l})$, whose values then determine $m_{l}$ via (39), and the
fluxes via (32) (see Fig. 7). As in the PS state (or any equilibrium-like
phase coexistence), the relative amounts of the two coexisting phases is set
by global conservation such that as the point of interest approaches one end
of the tie-line, the amount of the other phase vanishes linearly. Upon
crossing the curve $a_{5}(\rho_{0})$ set by (45) this causes a jump in the
second derivative of the rate function (27).

Figure 7: Geometrical construction for the high and low density phases in the
TB region. Here $\text{Pe}=3.5$.

Until now we did not discuss the dynamical equations (4) and (5) that must be
satisfied by any phase-separated solution. These are trivially met in the bulk
of the phases, and also at the interfaces in the PS state where there is no
mass flux in either coexisting phase. However, in the TB state, a stationary
solution would violate the conservative equation (4) since the high density
phase has zero mass flux, while the low density phase has a nonzero mass flux.
To resolve this issue, the entire phase separated solution must propagate with
velocity

$V=-\frac{J_{\rho,l}}{1-\rho_{l}},$ (46)

as indicated in Fig. 3. Meanwhile the dynamical constraint (5) can always be
accommodated with the help of the additional tumbling rate degree of freedom
$K$ which is nonvanishing across the sharp interface. This results in a
negligible $\mathcal{O}(\ell)$ contribution to the rate function that is
negligible here.

##  VIII Results for the smoothly modulated (SM) phase

The last region of the phase diagram to be analysed is the SM state, which
occurs in the green shaded regions of Fig. 4(A,B). As discussed in Sec. VI.3,
this state occurs in regions where the common tangent construction is
degenerate, which means that the behaviour as $\ell\to 0$ cannot be deduced by
considering bulk phases separated by sharp interfaces. At the level of the
minimisation problem (40), this can be seen by considering profiles such that
the local entropy production is (almost) typical for the local density

$\alpha(x)=\bar{a}(\rho(x))+o(\ell)\;.$ (47)

Inserting any such profile in (40) yields $I(a)=O(\ell^{2})$, where
$a=\int_{0}^{1}\alpha(x)dx$. That is, the “interfacial” (gradient) terms in
(40) determine which density profile minimises the action.

We therefore assume (47), so the minimization problem (40) reduces to a single
field minimization

$I(a)\simeq\ell^{2}\inf_{\rho(x)}\int_{0}^{1}\\!dx\,M(\rho)(\partial_{x}\rho)^{2},$
(48)

where $M(\rho)$ is derived in Appendix B, as

$M(\rho)=\frac{[2+\text{Pe}^{2}\,(1-\rho)(1-2\rho)]^{2}}{8\rho(1-\rho)[2+\text{Pe}^{2}(1-\rho)]}\;,$
(49)

and the constraints (41) become

$\rho_{0}=\int_{0}^{1}dx\,\rho(x)\qquad,\qquad
a=\int_{0}^{1}dx\,\bar{a}(\rho(x)).$ (50)

It can be verified that $M(\rho)\geq 0$. Notice that the spatial averages (50)
over any density profile $\rho_{0}(x)$ are bounded within the convex hull of
the curve $\bar{a}(\rho_{0})$ which comprises the SM region (42).

The next step is to analyse the single-field minimisation problem (48). This
can be performed numerically (indeed, that method was used to generate the SM
profiles in Fig. 3). However, we will show that analytic progress is also
possible, and allows phase transitions to be characterised. The integral
constraint $a=\int\alpha(x)dx$ makes the analysis more challenging so we
employ a standard tactic in large deviation theory Touchette (2009) which
corresponds to a change of ensemble. This is equivalent in this case to using
a Lagrange multiplier to enforce the constraint.

###  VIII.1 Scaled cumulant generating function

To deal with the constraint, we consider a scaled cumulant generating function
(CGF) Jack (2020); Touchette (2009) for a suitably rescaled version of the
IEPR $\mathbb{S}_{T}$,

$\psi(\Lambda)=\frac{1}{\gamma LT}\ln\left\langle
e^{T\Lambda\mathbb{S}_{T}/(\gamma\text{Pe}^{2})}\right\rangle,$ (51)

where the average is over the system dynamics. Since individual paths obey the
large deviation principle (10), the average can be computed by the saddle
point method, which gives

$-\psi(\Lambda)=\frac{1}{\gamma T}\inf_{\cal
X}\left[{\mathcal{A}}(\mathcal{X})-\Lambda\mathbb{S}(\mathcal{X})T/(\gamma\text{Pe}^{2}L)\right]$
(52)

Comparing with (28), we have exchanged a constrained minimisation over ${\cal
X}$ for a “biased minimisation” that includes the Lagrange multiplier
$\Lambda$. If one now optimises over trajectories ${\cal X}$ with
$\mathbb{S}_{T}({\cal X})=aL\gamma\text{Pe}^{2}$, one sees from (28) that

$\psi(\Lambda)=\sup_{a}\left[\Lambda a-I(a)\right]\;,$ (53)

which we recognise as a Legendre-Fenchel transform. (In the context of large
deviation theory, this is Varadhan’s lemma Touchette (2009), which relates the
CGF of $a$ to its rate function.) In the homogeneous state where $a<a_{1}$ one
has the Legendre transform of $I_{H}$

$\psi_{H}(\Lambda)=\left(\frac{\Lambda^{2}}{\text{Pe}^{2}}+\Lambda\right)\bar{a}(\rho_{0})\;.$
(54)

Computation of $\psi(\Lambda)$ also allows reconstruction of the rate function
by inverse Legendre transform as

$I(a)=\sup_{\Lambda}\left[\Lambda a-\psi(\Lambda)\right]\;.$ (55)

If the derivative $\psi^{\prime}(\Lambda)$ has discontinuities, some care is
required with this transformation Touchette (2009). Such situations correspond
to first-order dynamical phase transitions. In the system considered here,
(55) does hold even in the presence of such phase transitions; this leads to a
common-tangent construction for $I(a)$, associated with time-like phase
separation as we find in the subsequent Sec.VIII.4.

###  VIII.2 Stability of homogeneous state : Landau-like theory

We now consider the minimisation problem (52) and we analyse the stability of
the homogeneous state. The aim is to characterize the system state in the SM
phase and the associated phase transition. In contrast to previous Sections
where these were computed in the limit $\ell\to 0$, here these are determined
by the sub-leading $\mathcal{O}(\ell^{2})$ corrections. As a side result we
will derive also the boundaries of the (green) SM regions in Fig. 4 up to
$\mathcal{O}(\ell^{2})$.

To this end, we construct a variational ansatz for (52). It is sufficient to
consider trajectories $\cal X$ that are homogeneous in time (consistent with
the additivity principle). For parameters in the non-convex region given by
(42), we also take the local entropy production rate to satisfy (47). Then the
action ${\cal A}$ in (52) can be approximated as in (48) by ${\cal
A}(\mathcal{X})=\ell^{2}\int dxM(\rho)(\partial_{x}\rho)^{2}$. Since the
action is $O(\ell^{2})$ it is convenient to define
$\tilde{\Lambda}=\Lambda/\ell^{2}$ and
$\tilde{\psi}(\tilde{\Lambda})=\psi(\ell^{2}\tilde{\Lambda})/\ell^{2}$. Eq.
(52) finally becomes

$\displaystyle-\tilde{\psi}(\tilde{\Lambda})=\inf_{\rho(x)}\int_{0}^{1}dx\left[M(\rho)(\partial_{x}\rho)^{2}-\tilde{\Lambda}\bar{a}(\rho(x))\right].$
(56)

where the minimisation is subject (as usual) to
$\rho_{0}=\int_{0}^{1}dx\rho(x)$. Now consider a variational ansatz for the
density profile, as:

$\rho(x)=\rho_{0}+A\cos(2\pi x)+B\cos(4\pi x)\;.$ (57)

This ansatz is similar to that of Bodineau and Derrida (2007); Lecomte  _et
al._ (2012); Baek  _et al._ (2018); Dolezal and Jack (2019). It is sufficient
in this context because the instabilities of the homogeneous state occur in
the largest wavelength term (with coefficient $A$), but it will be necessary
to consider also the second harmonic (with coefficient $B$), to accurately
describe the behaviour of $A$ beyond leading order.

Injecting (57) into (56), the functional minimisation over $\rho$ reduces to a
minimisation over $A$ and $B$. Details of the computation are given in
appendix D. For homogeneous profiles ($A=B=0$), one reproduces the linear term
in the homogeneous CGF $\psi_{H}$ (54), as the quadratic part is sub-leading
$\mathcal{O}(\ell^{4})$ in the scaling regime $\Lambda\sim\ell^{2}$. Assuming
now that $A$ and $B$ are small, one finds that $B=O(A^{2})$ and this parameter
may be minimised analytically. Collecting terms up to quartic order in $A$,
(56) becomes

$\tilde{\psi}(\tilde{\Lambda})=\tilde{\psi}_{H}(\tilde{\Lambda}_{c})-\min_{A}\Psi(A)\,.$
(58)

where

$\Psi(A)=(\tilde{\Lambda}_{c}-\tilde{\Lambda})\bar{a}^{\prime\prime}(\rho_{0})\frac{A^{2}}{4}+\frac{\beta(\rho_{0},\text{Pe})}{2}A^{4}+\mathcal{O}(A^{6})$
(59)

is a Landau-like free energy with
$\tilde{\Lambda}_{c}=8\pi^{2}M(\rho_{0})/\bar{a}^{\prime\prime}(\rho_{0})$. A
formula for $\beta$ is given in Eq.(119) of Appendix D.

The boundaries of the SM state are obtained by computing values of
$(\rho_{0},\text{Pe},\tilde{\Lambda})$ for which $\Psi$ is minimised by a
value $A\neq 0$, corresponding to an inhomgeneous density. Since $\Psi$ is
known to $\mathcal{O}(A^{4})$ we can perform this calculation for small $A$.
The behaviour of this Landau-like theory depends on the signs of
$\beta(\rho_{0},\text{Pe})$ and $\bar{a}^{\prime\prime}(\rho_{0})$, as well as
that of $(\tilde{\Lambda}-\tilde{\Lambda}_{c})$. We separate two cases
according to the sign of $\beta$.

###  VIII.3 Linear instabilities and second-order phase transitions
($\beta>0$)

Assume first that $\beta>0$. In this case, analysing the stability of the
homogeneous state only requires that we consider (58) up to quadratic order in
$A$. From (58), the system will be inhomogeneous ($A\neq 0$ and
$\tilde{\psi}>\tilde{\psi}_{H}$) if

$(\tilde{\Lambda}-\tilde{\Lambda}_{c})\bar{a}^{\prime\prime}(\rho_{0})<0.$
(60)

From (26) then

$\bar{a}^{\prime\prime}(\rho_{0})=6\rho_{0}-4$ (61)

which changes sign at $\rho=2/3$. Note also that $\tilde{\Lambda}_{c}$ has the
same sign as $\bar{a}^{\prime\prime}(\rho_{0})$. Using these facts with (54),
the result is that the homogeneous phase has two regions of linear
instability, which must correspond to $A\neq 0$, see also Appendix D.1. The
first occurs for $\rho_{0}<2/3$, it is
$a<a_{2}(\rho_{0},\text{Pe},\ell)<\bar{a}(\rho_{0})$ with

$a_{2}(\rho_{0},\text{Pe},\ell)=\left[1+\ell^{2}\frac{8\pi^{2}M(\rho_{0})}{\text{Pe}^{2}\,(3\rho_{0}-2)}\right]\bar{a}(\rho_{0}).$
(62)

The second linearly unstable region occurs for $\rho_{0}>2/3$ and corresponds
to $a>a_{2}(\rho_{0},\text{Pe},\ell)>\bar{a}(\rho_{0})$, where $a_{2}$ is
given by the same formula (but note this function has two branches). See Fig.
8 and Appendix D.1 for further details. In both cases, the minimum of the
Landau-like free energy occurs at $A=\bar{A}$ with

$\bar{A}=\sqrt{\frac{(\tilde{\Lambda}_{c}-\tilde{\Lambda})\bar{a}^{\prime\prime}(\rho_{0})}{4\beta(\rho_{0},\text{Pe})}}.$
(63)

This is a continuous transition, with
$\bar{A}\sim|\tilde{\Lambda}_{c}-\tilde{\Lambda}|^{1/2}$.

To put these results in context, we note from Lecomte  _et al._ (2012); Baek
_et al._ (2018); Dolezal and Jack (2019); Jack  _et al._ (2015) that for
functional minimisations of the form (56), one expects inhomogeneous states to
form for either positive or negative values of the biasing parameter
$\Lambda$, according to the sign of $a^{\prime\prime}$. In contrast to
previous studies, the function $a^{\prime\prime}$ considered here changes
sign, leading to two separate regions of instability. We will see in the
following that the situation is in fact more complex than this linear analysis
can reveal.

Figure 8: Illustration of the boundaries of linear instability dented by cyan
line. They are given by the curve $a_{2}(\rho_{0})$ (62). As $\ell\to 0$ these
tend to the boundaries of local non-convexity of $I_{H}$ denoted by orange
shading. The purple shading denote regions where $I_{H}$ is globally non-
convex yet locally convex, see also Fig. 5 (B) and (C). The black curve is the
expected IEPR curve (26). Here $\text{Pe}=2.5$ and $\ell=0.1$.

As anticipated above, one also sees that for $\ell\to 0$, the critical curve
(62) tends to the boundary of the region where $I_{H}$ is locally convex,
which is given by $(2/3-\rho_{0})(a-\bar{a}(\rho_{0}))<0$, see Fig. 8. This
arises from (62) by virtue of the fact that the $\ell^{2}$ piece has
coefficient $\pm\infty$ for $\rho_{0}\to 2/3$, shown as a vertical orange line
in Fig. 8, so that as $\ell\to 0$, the stability boundary is the union of this
vertical line with the unbiased IEPR curve $\bar{a}(\rho_{0})$. In
consequence, at $\rho_{0}<2/3$ the solution is unstable for vanishingly small
suppressions of the entropy production,
$a<a_{2}=\bar{a}-\mathcal{O}(\ell^{2})$, while at $\rho_{0}>2/3$ it is
unstable for vanishingly small enhancements,
$a>a_{2}=\bar{a}+\mathcal{O}(\ell^{2})$.

###  VIII.4 Beyond linear stability ($\beta<0$): First-order phase transition
and timelike phase coexistence

To go beyond linear stability analysis, we return to the Landau-like free
energy (59). If the coefficient $\beta<0$ [but assuming that the theory is
ultimately stabilised by terms at $O(A^{6})$ or higher], this minimisation
problem can support first-order phase transitions. That is, while the
homogeneous state may be a local minimum of $\Psi$, there may be additional
minima corresponding to SM states. On varying the parameters, the global
minimum can then change discontinuously from homogeneous to SM, in direct
analogy with first-order thermodynamic phase transitions. Such transitions
occur while the homogeneous state is locally stable: they cause the SM states
to appear for even smaller IEPR deviations than (62) would predict, see also
Bodineau and Derrida (2007).

To see this, recall (63). For $\beta>0$ we have a continuous phase transition,
as noted above. For $\beta<0$, higher-order terms in $\Psi$ must be
considered, in order that the minimum can be found. However, the order of the
phase transition can be deduced without computing these terms. For
convenience, consider the case $\rho_{0}>2/3$ so that
$\bar{a}^{\prime\prime}(\rho_{0})>0$. Suppose that $\beta<0$ and
$\tilde{\Lambda}$ is increased from zero towards $\tilde{\Lambda}_{c}$. For
$\tilde{\Lambda}=0$ the free energy $\Psi$ has its global minimum at zero. At
$\tilde{\Lambda}=\tilde{\Lambda}_{c}$, this minimum becomes unstable: the
negativity of $\beta$ shows that the global minimum of the Landau-like free
energy must be already at $A\neq 0$. That is, the global minimum of $\Psi$
must become non-zero for some $\tilde{\Lambda}_{1}<\tilde{\Lambda}_{c}$, at
which point the local minimum at zero is still locally stable. This
corresponds to a first-order (discontinuous) phase transition with a jump
discontinuity in the magnitude $A$. For $\rho_{0}<2/3$, a similar analysis
holds, except that $\tilde{\Lambda}_{c}<0$ so one must consider the behaviour
on reducing $\tilde{\Lambda}$ towards $\tilde{\Lambda}_{c}$. Within the Landau
theory, this argument shows that points with
$(\Lambda_{c},\beta)=(\tilde{\Lambda}_{c},0)$ are tricritical.333Technically,
they might be critical points of even higher order, since we have not checked
the positivity of the $A^{6}$ terms in (59), but this does not affect whether
the transition is continuous or discontinuous, which is our main concern here.
We find that for $\text{Pe}<\text{Pe}^{*}\simeq 2.641$ the Landau theory
supports two such points $\rho_{c1,2}$, one at $\rho_{c1}<2/3$ and one at
$\rho_{c2}>2/3$. For $\text{Pe}>\text{Pe}^{*}$ additional points appear on
both sides of $\rho_{0}=2/3$, see also Sec.X. Fig. 9, shows the jump
discontinuity in the spatial modulation amplitude
$\Delta\rho=\rho_{\text{max}}-\rho_{\text{min}}$ which grows smoothly when
advancing along the first order transition line away from the tricritical
point $\rho_{c1}$.

Figure 9: The spatial modulation amplitude
$\Delta\rho=(\rho_{\text{max}}-\rho_{\text{min}})$, at different values of
averaged density $\rho_{0}$ in purple points. these are obtained from
numerical solution of (56). The black dot marks the critical value for linear
instability. The black curves are the theoretical prediction close to
criticality $\Delta\rho=2A$ with $A$ given by (63). The point $\rho_{0}=0.58$
is just below the tricritical value $\rho_{c1}=0.600\dots$, while
$\rho_{0}=0.61$ is just above it. For densities $\rho_{0}<\rho_{c1}$ the
transition is continues and second order. For densities $\rho_{0}>\rho_{c1}$
the transition is first order and occurs before the critical value for linear
instability $|\Lambda_{1}|<|\Lambda_{c}|$. Here $\text{Pe}=2$ and $\ell=0.02$.

Figure 10: (A) Dynamical phase diagram for representative parameters
$\text{Pe}=1$ and $\ell=0.02$. Phases are labelled as in Fig. 4. Orange solid
lines: second order transitions between H and SM, given by (62). Dashed orange
and magenta lines: boundaries of the miscibility gap associated with
discontinuous transitions between H and SM. Black circles: tricritical points
that separate first-order and second-order transitions. Black dashed lines:
tie-lines connecting points at equal density, within the miscibility gap. (B)
Illustration of the miscibility gap for $\text{Pe}=2$, $\rho_{0}=0.74$ and
$\ell=0.02$. Solid black line: rate function $I(a)$. Dashed blue line: the
solution of the minimisation problem (48), which assumes trajectories are
homogeneous in time, which shows a kink at the point where the minimiser of
that problem changes from H to SM. The true rate function is obtained by (55)
as the convex hull of the dashed blue line: the dominant trajectories in this
case exhibit time-like phase separation. Orange and magenta dots mark the
edges of the miscibility gap. As $\ell\to 0$, the rate function is
$O(\ell^{2})$ throughout this range, the orange dot approaches $\bar{a}$,
while the magenta dot approaches a finite value of $a$, so the miscibility gap
stays open.

The resulting dynamical phase diagram is shown in Fig. 10(A). We focus on the
behaviour for $a<\bar{a}$ (below the black line in that Figure), which
includes the tricritical point at $\rho_{0}<2/3$. The orange line indicates
the linear instability of the H state, which is one boundary of the SM phase.
This line terminates at the tricritical point at $\rho_{c1}\simeq 0.55$ (black
dot), where $\beta=0$. Beyond this point, the transition between H and SM
states is first-order, corresponding to a jump discontinuity in $a$ as
$\tilde{\Lambda}$ is reduced from zero. The two ends of the jump are indicated
by orange and magenta dashed lines. For values of $a$ below the magenta dashed
line, the system is SM.

The first order transition means that the rate function defined by (48) is not
convex when viewed as a function of the single variable $a$ but instead has a
kink singularity along the $a$ axis, see Fig. 10 (B). Similarly, the CGF
defined by (56) has a kink singularity along the $\tilde{\Lambda}$ axis.
Correspondingly, the relation (53) is not invertible Jack (2020); Touchette
(2009). Acting with a Legendre-Fenchel transform on $\tilde{\psi}$ we will now
arrive at the convex envelop of (48) obtained by a common tangent construction
over the ‘miscibility gap’ in $a$. This results in the final expression for
the true rate function (27), which does not equal the expression (48), see
Fig. 10 (B). This construction corresponds to a timelike phase coexistence
between the magenta and orange dashed lines: the optimal path in (28) spends
part of its time period in the H state and the remainder in the SM state Jack
(2020). These time periods are separated by instantons (domain walls in time)
whose contribution to the action is negligible Bertini  _et al._ (2006); Baek
_et al._ (2017, 2018) in the limit $T\to\infty$. In the thermodynamic context,
the domain of phase coexistence in the phase diagram is called a miscibility
gap.

All these phenomena – tricritical point, miscibility gap, and time-like phase
coexistence – are duplicated for $a>\bar{a}$, although the fine detail is not
graphically resolved in Fig. 10(A). In this regime, it is evident from Fig.
10(A) that the transition from homogeneous to inhomogeneous profiles may
proceed directly to PS state, skipping the SM state altogether (possibly
subject to a small miscibility gap). However, when $\ell$ becomes small, the
orange lines (both solid and dashed) approach the black curve,
$a=\bar{a}(\rho_{0})$, and the first instability of the H state is always to
an SM state. On the other hand, the miscibility gap remains finite even as
$\ell\to 0$. Indeed, its endpoints correspond to the minima of $\Psi$ (59)
which are $\ell$ independent, as is the gap between these endpoints. For
further discussion see Appendix D.4.

Finally, we note that timelike phase coexistence is only possible when
$I(\rho_{0},a)$ is concave in $a$ at fixed $\rho_{0}$, so that the two phases
that correspond to the ends of the miscibility gap can share the same
$\rho_{0}$. This is why we did not consider such solutions when discussing the
convex hull of $I_{H}$.

###  VIII.5 Character of the H $\to$ SM and SM $\to$ PS transitions

The preceding arguments based on (56) and (58) are sufficient to derive the
boundaries of the SM region in Fig. 10(A), as advertised above. We give a
brief summary of the underlying physics.

For $a=\bar{a}$ (typical values of the IEPR) then the system is homogeneous
(H). For any $a<\bar{a}$, the convex hull of $I_{H}$ is zero, so one expects
inhomogeneous states to appear as $\ell\to 0$, but the structure of these
states cannot be predicted by the quasi-thermodynamic (common tangent)
approach used in Sec. VII. Instead, the physics in this regime is entirely
under the control of the gradient terms in (40). The result is that the H
state survives in a narrow region with $a$ close to $\bar{a}$, beyond which
the system becomes SM. This may take place either by a linear instability
(continuous transition) or by a finite jump in the density profile, from a
constant density to one with a finite modulation amplitude. In this latter
case, the phase diagram includes a miscibility gap, where SM and H phases
coexist within a single trajectory (in different time periods). The boundary
between the continuous and discontinous transitions is a tricritical point.

The distance of the H$\to$SM transition from the line $a=\bar{a}$ is
determined by the balance of two costs in the minimisation problem (40). The H
state has a cost from the reduction in IEPR at fixed $\rho$, coming from the
bulk term $I_{H}$, while the SM state has an interfacial cost proportional to
$\ell^{2}$. At moderate densities, the costs balance for
$\bar{a}-a=O(\ell^{2})$, leading to a narrow band of the H phase between
$a=\bar{a}(\rho)$ and the SM state in Fig. 10. However, if the density is
large enough that $\bar{a}(\rho_{0})=O(\ell)$, this band of the H phase can
extend all the way to $a=0$. This is a non-perturbative effect of $\ell$ on
the phase boundary; the range of densities over which it can occur must vanish
as $\ell\to 0$. [In that limit, Fig. 4 already shows that the SM state extends
throughout the range $0<a<\bar{a}(\rho_{0})$.]

For $a>\bar{a}$ and $\rho>1/2$, one observes again an H$\to$SM phase boundary
on which a tricritical point separates regions with first-order and second-
order phase transitions. In this case, the H state occurs when
$a-\bar{a}=O(\ell^{2})$. In addition, the system supports PS states for
$a>a_{3}(\rho_{0})$, as defined in (43). On approaching this line, the SM
solution develops increasingly sharp gradients in order to accommodate the
constraints (41) and the assumption (47) simultaneously: the SM profiles start
to resemble coexisting bulk phases whose densities are approximately
$\rho_{l,h}=1/2,1$. At the SM to PS transition, the interfaces sharpen (their
widths become $O(\ell)$), while the coexisting densities are (exactly)
$\rho_{l}=1/2$ and $\rho_{h}=1$, consistent with common tangent construction
of Sec. VII. Note that the SM and PS states are distinguished by whether their
interfacial widths are sharp $O(\ell)$ or smooth $O(1)$: this distinction is
only meaningful in the limit $\ell\to 0$. For finite $\ell$, there is a smooth
crossover between these states.

In fact, a similar SM$\to$PS transition also occurs upon approaching the $a=0$
line from the the SM phase with $a>0$. Here too, the inhomogeneous SM solution
must develop increasingly sharp gradients until becoming sharply phase
separated with co-existing densities $\rho_{l,h}=0,1$. However, this PS state
only occurs when $a=0$ exactly. In the next Section we address what lies
beyond the $a=0$ line, in the region at negative $a$.

##  IX The rate function and its extension to the negative half plane

The preceding sections provide a full description of the nonequilibrium phase
diagram for positive values of the IEPR, $a>0$. This half plane contains
regions of both positive and negative deviation from the unbiased value
$\bar{a}(\rho_{0})$. Here we consider the resulting minimized rate function.
We also extend both the rate function and the phase diagram to $a<0$.

###  IX.1 The small $\ell$ limit

Although the phase diagram for $a>0$ is complicated, even in the limit
$\ell\to 0^{+}$ (Fig. 4 (A,B)) the resulting rate function $I(\rho_{0},a)$
becomes simple in this limit: it everywhere approaches the convex hull of
$I_{H}$ obyeing (38). For $\rho_{0}<1/2$, this hull is given by

$\displaystyle I_{<}(a)=$ $\displaystyle 0,$ $0\leq a\leq\bar{a}$,
$\displaystyle I_{<}(a)=$ $\displaystyle I_{H}(\rho_{0},a),$ $\bar{a}\leq a$,
(64)

and for $\rho_{0}>1/2$:

$\displaystyle I_{>}(a)=$ $\displaystyle 0,$ $0\leq a\leq a_{3}$,
$\displaystyle I_{>}(a)=$
$\displaystyle\frac{\text{Pe}^{2}}{(1-\rho_{0})}(a-a_{3})^{2},$ $a_{3}\leq
a\leq a_{4}$,  $\displaystyle I_{>}(a)=$
$\displaystyle\frac{1-\rho_{0}}{1-\rho_{l}}I_{H}(\rho_{l},a_{l}),$ $a_{4}\leq
a\leq a_{5}$,  $\displaystyle I_{>}(a)=$ $\displaystyle I_{H}(\rho_{0},a),$
$a_{5}\leq a$. (65)

with $\rho_{,}a_{l}$ obtained by the geometrical construction that is shown in
Fig. 7, see Sec.VII. These expressions for $I_{H}$ are plotted in Fig. 11.

---
Figure 11: The large deviation rate function $I(a)$ (27) in thick black line
for $\rho_{0}=0.3$ (A) where it is given by (64), and for $\rho_{0}=0.74$ (B)
and (C) where it is given by (65). The dashed black lines are the quadratic
form in (38) which serves as an upper bound on the rate function $I$
everywhere in parameter space. The colored points mark the critical values
that correspond to the critical phase transition curves in Fig. 4. In all
panels $\text{Pe}=3.5$. All the results shown are for $\ell=0^{+}$.

These results can be extended to $a<0$ using the detailed fluctuation theorem
(29). The corresponding trajectories (which solve (28)) are given by simple
time reversal (16) of the corresponding points with $a>0$. As a result, all
the critical lines $a_{1-5}$ and the various phases that they separate are
simply reflected to the negative half plane. Note that the additional linear
term in (29) does not affect the construction of phase equilibria within this
half plane (it has no effect on concavities there) and moreover, because its
coefficient is positive, it cannot create additional concavities linking the
two half planes either. The full minimization problem for the rate function is
thus solved independently in each half plane, and the phase diagram for $a<0$
is simply the reflection, in the $a=0$ axis, of that found already for $a>0$ .

Examples showing the resulting rate function (for $\ell\to 0$) are presented
in Fig. 12: these are composites of direct and mirrored segments for positive
and negative $a$. This composite function exhibits a kink singularity at
$a=0$. Indeed in the SM region with $\bar{a}(\rho_{0})>a>0$ we have that
$I=0$, while its reflection to the negative half plane has
$I=\text{Pe}^{2}|a|$. The two surfaces meet along the line $a=0$, with
discontinuous first derivative.

The fluctuation theorem (29) applies for all $\ell$, so the dynamical phase
diagrams shown in Fig. 4 (C,D) are again mirrored in the lower half plane. The
identities of the phases are the same and the corresponding trajectories are
related by time reversal.

---
Figure 12: The fluctuation theorem construction (29) extending our results to
the negative half plane $a<0$. All critical values are reflected accordingly.
The large deviation function $I(a)$ (27) in thick black line for
$\rho_{0}=0.3$ (A) and for $\rho_{0}=0.74$ (B). In both panels
$\text{Pe}=3.5$. All the results shown are for $\ell=0^{+}$. In panel (A) the
rate function in solid line follows the dashed parabola (the quadratic from in
(38)) in the two H phases (between the red and green points), where in panel
(B) there is no H phase and the rate function differs from the parabola
everywhere.

###  IX.2 Finite $\ell$ corrections

We have discussed the role of finite $\ell$ for these phases and transitions
at various points in the preceding sections, but summarize them here for ease
of reference. These effects are quite subtle and responsible for all the
differences between the diagrams in panels (A,B) of Fig. 4 and those in (C,D)
of the same figure, some of which are nonperturbative, such as the re-entrant
H phase in the lower right corner. The reason for these relatively large
effects is the degeneracy whereby the interfacial terms dominate the action
everywhere within the two-dimensional convex hull of the $\bar{a}(\rho_{0})$
curve. In most other regions of the phase diagram, the parameter $\ell$ has
(at most) a perturbative effect on the phase boundaries. For example, the
transition between H and CM states occurs at $a=a_{1}(\rho_{0},\text{Pe})$, as
given in (36). Since these states both have homogeneous density, this phase
boundary is not affected by $\ell$.

On the other hand, the parameter $\ell$ does affect the phase boundaries of
the SM state: as discussed in Sec. VIII.4 and Fig. 10, increasing $\ell$
diminishes the SM regions of the phase diagram, making way for the re-entrant
H phase as just described. The SM phase is diminished by $\ell$ because its
associated interfacial cost (proportional to $\ell^{2}$) may become greater
than the cost of a homogeneous state with mildly reduced IEPR. The tricritical
points on the boundary of the SM phase are present in the limit $\ell\to 0$,
while their positions depend perturbatively on $\ell$. The widths of the
miscibility gaps are finite as $\ell\to 0$, and in general they depend non-
perturbatively on $\ell$.

Turning to the SM $\to$ PS transition (which occurs for $a=a_{3}(\rho_{0})$,
as given by (43)): the distinction between the SM and PS states is whether
there are sharp interfaces between bulk phases or a smooth modulation of the
density. This is a clear distinction only in the limit $\ell\to 0$, since
sharp interfaces are never observed for $\ell>0$. Hence, the SM$\to$PS
transition becomes a smooth crossover when $\ell$ is finite. A similar
rounding occurs of the kink that appears at $a=0$, when the rate functions for
the upper and lower half plane get glued together; the discontinuous
derivative described in Sec. IX.1 is smoothed across a boundary layer of width
$\ell$.

Lastly, the domains of the PS and TB phases are diminished due to a finite
$\mathcal{O}(\ell)$ cost in the action for sharply phase separated solutions.
This means that the transitions at $a=a_{4}$ (44) and $a=a_{5}$ (45) are
perturbatively distorted on increasing $\ell$ from zero, leaving the nature of
the transition unchanged.

##  X Tricriticality near the MIPS critical point

All the work reported above addresses the large deviation structure for IEPR
fluctuations in regimes of the active lattice gas model where the unbiased
dynamics results in a homogeneous (H) phase. On the one hand, as found above,
the biased H state loses linear stability beyond the critical IEPR level set
by $a_{2}$ (62). On the other hand, the same model has a MIPS critical point
at $(\rho_{0}=0.75,\text{Pe}=4)$, where the unbiased homogeneous state loses
linear stability towards an inhomogeneous (in this case PS) state. That is,
the MIPS transition in the unbiased system and the transition from H to SM
under bias share the same sort of criticality. Notice that this is not true
for any of the other phase transitions shown in Fig. 4, which do not arise
from the loss of linear stability. A similar property was reported in the
study of current statistics of the same model Agranov  _et al._ (2022).
However, the transition in that work was into a traveling band state which
breaks reflection symmetry, in contrast to the SM state here which does not.

The MIPS criticality can then be interpreted as the point where the biased
system is critical at exactly zero bias. Indeed, the difference
$a_{2}-\bar{a}$ vanishes when approaching the MIPS critical point. This can be
verified from the expression (62) for $a_{2}$, because the mass term $M$ (49)
entering this expression vanishes at the MIPS critical point. This behavior
was anticipated without proof in a previous study Nemoto  _et al._ (2019); the
current work provides the first definitive substantiation of such a scenario.

Figure 13: The difference $(a_{2}-\bar{a})/\ell^{2}$ (62) tends to zero when
approaching the MIPS criticality at $(\rho_{0}=0.75,\text{Pe}=4)$ denoted by
the red dot. Here the solid line denotes the second order transition curve
$a_{2}$ (62), while the dashed line denotes its extensions beyond the
tricritical points (the black dots) where the transition is first order and is
not captured by the $a_{2}$ curve. Two tricritical points collide at the MIPS
critical point of the unbiased model. The leftmost of them is the one analyzed
before in Sec. VIII.4 and shown in Fig. 4 and Fig. 10. Above
$\text{Pe}>\text{Pe}^{*}\simeq 2.64$ two additional tricritical points appear
with the same phenomenology as discussed before. Of these two points, the one
with a lower $\rho_{0}$ is the one of the two points that collide at the MIPS
criticality.

Intriguingly, we find that at the MIPS criticality condition two tricritical
points collide as shown in Fig. 13 and verified analytically in Appendix D.5.
One of these points is the one we analyzed in Sec. VIII.4 above. The second
one is an additional tricritical point that emerges at higher values of Pe
than we considered thus far, see Fig. 13. In terms of the biased phase
diagrams, this collision occurs only at one specific Peclet number (Pe $=4$),
the one for which the unbiased MIPS is critical. In this setting it is a
higher order critical point of some kind.

Note however, that this does not imply that the MIPS critical point itself is
any kind of tricritical or higher order critical point of the unbiased system.
Indeed, an ordinary critical point on an unbiased phase diagram might become
first order at infinitesimal bias, but this does not make the unbiased model
tricritical.

##  XI Discussion : Comparison with numerical studies of off-lattice models

We have studied the large deviation theory for fluctuations in the informatic
entropy production rate (IEPR) of an active lattice gas model. Remarkably, the
model admits an exact analytical derivation of the IEPR large deviation
function (64) and (65), and the corresponding phases, see Figs. 3 and 4. It
displays a rich phenomenology with an intricate phase diagram unlike any
reported before in studies of related passive lattice gases, or other active
systems. The IEPR fluctuations are characterized by five different phases that
correspond to the breaking of different symmetries in the system. All these
exact derivations are possible in the limit of small diffusive length scale
$\ell\ll 1$ which is where the system also admits MIPS. (Note however that we
addressed only parameter values where the unbiased system is homogeneous
rather than showing MIPS.) Yet finite $\ell$ corrections leave most aspects of
the qualitative picture unchanged, albeit with some stronger effects in
regions of the phase diagram where the rate function vanishes in the absence
of these corrections. These can in turn be analysed in detail to leading order
in $\ell$.

The nature of the dynamical phases and the associated phase transitions have
been discussed extensively in previous sections, and we will not repeat them
here. However, it is useful in what follows to compare the results for this
one-dimensional lattice-based model with off-lattice two-dimensional models
for which large deviations of the IEPR have been estimated numerically
GrandPre  _et al._ (2021); Keta  _et al._ (2021); Nemoto  _et al._ (2019); Yan
_et al._ (2022).

###  XI.1 Fast vs slow magnetisation fields

When comparing off-lattice models with the active lattice gas, it is important
to recall that the tumbling rate in the lattice model is slower than particle
hopping by a factor of $L^{2}$. In terms of the hydrodynamic description
(4,5), this means that both the density and the magnetisation enter on an
equal footing: the magnetisation is a slow field in the hydrodynamic sense. By
contrast, the off-lattice models are defined such that the time to move a
diameter is of the same order (in system size) as the time for decorrelation
of their direction of propulsion. This means that the corresponding
magnetisation fields for the off-lattice systems relax quickly; they are not
hydrodynamic.

This difference is important for several reasons: it helps to make the lattice
model tractable for analytical calculations, and it leads to large deviation
probabilities with different dependences on the system size. However, the main
phenomenology of the models’ steady states is similar, including MIPS. To be
precise about the scaling of the large deviation probabilities, it is useful
to return to the microscopic time variable $\hat{t}$, which we recall is
defined such that a time period $\hat{t}=1$ is long enough that a typical
particle will diffuse by one lattice spacing. Expressed in this way, the
probability to observe IEPR $a$ scales in the off-lattice models as

$-\ln P(a)\sim\hat{T}L^{d-2}I(a)$ (66)

where the dimension $d=1$, we used (27), and the trajectory length in
microscopic time units is $\hat{T}=TL^{2}$.

In off-lattice systems, there are two possibilities for the scaling of large
deviations: either (66), or alternatively

$-\ln P(a)\sim\hat{T}L^{d}I_{m}(a)\,$ (67)

for some rate function $I_{m}$. This latter form corresponds to rare events
with much smaller probability. The scaling (66) occurs when the large
deviation occurs by a hydrodynamic mechanism, which means that only the slow
hydrodynamic fields behave in a non-typical way during the rare event. For the
active lattice gas model, all the large deviations that we have considered
occur by hydrodynamic mechanisms and obey (66). On the other hand, large
deviations that require non-typical behaviour in fast (non-hydrodynamic)
fields are expected to follow (67). In the off-lattice models of Keta  _et
al._ (2021); Nemoto  _et al._ (2019), some large deviations do occur by such
mechanisms, because the magnetisation is now a fast field.

When considering biased ensembles, it is important that the bias $\Lambda$ was
chosen in (51) under the assumption that large deviations scale as in (66). In
contrast, the corresponding bias $s$ in Keta  _et al._ (2021); Nemoto  _et
al._ (2019) was defined under the assumption that (67) holds. To access
hydrodynamic scaling regimes one should define a corresponding quantity
$\Lambda=sL^{2}$. For a unified treatment of lattice and off-lattice models we
could choose to consider large deviations of the IEPR by taking the limit
$\hat{T}\to\infty$ at fixed $L$, and then consider the behaviour of the
corresponding rate function as $L\to\infty$. This is a different limiting
procedure to that described in Sec. IV.1, but it leads to the same results.

###  XI.2 Dynamical phase transitions: inhomogeneous states and collective
motion

In off-lattice models, Ref. Nemoto  _et al._ (2019) identified two dynamical
phase transitions. In the present context, those results correspond to
inhomogeneous states for $a<\bar{a}$, and states with collective motion that
occur for $a\geq a_{c}>\bar{a}$, where $a_{c}$ is a critical IEPR. Fig. 4
shows that this is similar (at least qualitatively) to the behaviour of the
active lattice gas for low density, $\rho_{0}<1/2$. This is the connection
that we explore in the following. (To our knowledge, the off-lattice models do
not support complex phase behaviour like that of the active lattice gas for
$\rho_{0}>1/2$, this point is discussed further below.)

Based on the discussion above, the expected scaling for large deviations in
the off-lattice models is (66) for the inhomogeneous states (which occurs by a
hydrodynamic mechanism), and (67) for the collective motion state (which
requires the participation of the fast magnetisation field Keta  _et al._
(2021)). This is consistent with the results of Nemoto  _et al._ (2019): the
inhomogeneous states were analysed in that work by assuming (66) with an
$L$-dependent rate function $I_{m}$, which is one possible way to encapsulate
both (66) and (67) in a single framework.

Since the magnetisation of the active lattice gas is a slow field, the
transition to CM behaviour is hydrodynamic in this model. This is different
from the off-lattice models. However, it was already argued in Sec. V that the
microscopic origin of the transition is similar in both cases: the IEPR is
suppressed by “collisions” between particles with different orientations, and
CM states reduce the possibilities for such collisions.

---
Figure 14: (A) and the blow up panel (B) are schematic representation of the
scaling behavior for the rate function $I$ for $\rho_{0}<0.5$ where the
governing phases are H, SM and CM. The dashed curve is the upper bound given
by the quadratic form in (V). The dotted straight line $\text{Pe}^{2}|a|$
marks the $\ell\to 0^{+}$ asymptotic of $I$ that is obtained by the detailed
fluctuation theorem (29). The corresponding kink singularity at $a=0$ is
smoothed over a narrow region around $a=0$ of width $\ell$ which is
characterized by a PS phase.

The relationship between the inhomogeneous states of the active lattice gas
and the off-lattice models is more subtle. The schematic Fig. 14 summarizes
the scaling behaviour of the rate function $I(a)$ when $\rho_{0}<1/2$ and
$\ell$ is small but finite. Different parts of the rate function have
different scalings with $\ell$: it is $O(\ell^{2})$ in the SM region and
$O(\ell)$ in the PS state near $a=0$. For $a>\bar{a}$ it is $O(1)$, including
for the CM state. Interestingly, the arguments in Keta  _et al._ (2021);
Nemoto  _et al._ (2019) about the behaviour of the ($L$-dependent) rate
function $I_{m}$ for the off-lattice models yield a similar picture: one has
$I_{m}(a)=O(1)$ for $a>\bar{a}$, there is an inhomogeneous (SM) regime for
reduced IEPR where the rate function $I_{m}(a)=O(L^{-2})$, which is equivalent
to a hydrodynamic scaling (66). There also is a phase separated (and
dynamically arrested) regime with sharp interfaces, where the IEPR $a\approx
0$ and $I_{m}(a)=O(L^{-1})$.

In other words, the behaviour of the active lattice gas and the off-lattice
models appear similar if one identifies the dimensionless parameter $\ell$ in
the lattice picture with the ratio $\sigma_{0}/L$, where $\sigma_{0}$ is the
particle diameter, in the off-lattice models. Recalling that
$\ell=\hat{\ell}_{D}/L$ where $\hat{\ell}_{D}$ is the diffusive length that
sets the natural scale for density gradients, this correspondence is natural.
The particle diameter $\sigma_{0}$ sets the natural scale for density
gradients in the off-lattice models, it remains fixed in the large-system
limit $L\to\infty$, so one generically expects phase separation into states
with sharp interfaces. The corresponding length scale $\hat{\ell}_{D}$ in the
active lattice gas is taken to be proportional to the system size $L$, in
order to facilitate the analytic computation, but one later assumes that
$\ell=\hat{\ell}_{D}/L\ll 1$, in order to restore the possibility of sharp
interfaces.

For this reason, the physical arguments in this work about the scalings shown
in Fig. 14 can be similarly applied to the off-lattice models, leading to the
same scalings, but with $\ell$ replaced by $\sigma_{0}/L$. This is indeed
consistent with Keta  _et al._ (2021); Nemoto  _et al._ (2019). It shows (as
expected) that the cost in probability of generating inhomogeneous density
profiles is controlled by the natural scale for density gradients, while the
transition to CM phase is independent of this parameter, leading to a rate
function that is $O(1)$ in the relevant parameter.

Recall that in the Introduction we described the $L$-dependence chosen for the
tumbling rate in the active lattice gas as ‘unnatural’ in that it causes the
magnetisation to become a slow field. However, this discussion shows that the
resulting phase behaviour is nonetheless ‘natural’, given a careful
interpretation of the control parameters. Specifically, many features of the
phase behaviour then echoe what is known about off-lattice models for which
magnetisation relaxes much more quickly. This reinforces our view that the
advantages of the active lattice gas model, in terms of what can be calculated
exactly, outweigh any disadvantage in terms of its ‘unnaturalness’.

##  XII Outlook and conclusions

We end with a few remarks on possible future directions. As noted above, the
active lattice gas behaves similarly to the off-lattice models for densities
$\rho_{0}<1/2$, but the lattice model has complex behaviour for higher
density, which has not been seen off-lattice Nemoto  _et al._ (2019); GrandPre
_et al._ (2021); Keta  _et al._ (2021); Nemoto  _et al._ (2019); Yan  _et al._
(2022). Much of the complexity of the active lattice gas behavior can be
traced back to the nontrivial density dependence of the expected IEPR
$\bar{a}(\rho_{0})$ (26), whose curvature changes sign at $\rho_{0}=2/3$. It
would be interesting to investigate off-lattice models where the mean IEPR
shows a similar property, which might well lead to novel phases such as the TB
reported here, as well as to tricritical points and other interesting
behaviour. (More specifically, the TB state operates via a combination of
phase separation with collective motion; it would be natural to see such
combinations in a system where $\bar{a}^{\prime\prime}$ can change sign.)

With this in mind, we recall from Sec. III.1 that the function
$\bar{a}(\rho_{0})$ can also be derived directly from the microscopic
dynamics, so the microscopic EPR coincides with the macroscopic IEPR. This is
because, as shown in Appendix A, the particles in system trajectories that
contribute to IEPR traverse macroscopic length scales over macroscopic time
scales, and are thus exactly captured by the coarse-grained fluctuating
hydrodynamics description. Consequently, the coarse-grained IEPR also fixes
the rate at which work is expended in the system by propulsive forces. This
equivalence is in accord with Ref. Kaiser  _et al._ (2018) which points at a
relation between the MFT large deviation functional, and the one corresponding
to the underlying Markov process.

All our analysis in this work assumed that the steady state of the active
lattice gas was homogeneous, and MIPS has not occurred. Nevertheless, we
anticipate that to leading order at small $\ell$, most features of the problem
go through to the MIPS phase. In particular, the rate function $I$ should
still be given by the convex hull of $I_{H}$. In all phases except the SM one,
we have considered possible sharply phase-separated solutions where the
optimal one among those is singled out via a convex hull construction. The
MIPS state can be regarded as one such possibility, and hence we have not
excluded it as part of the minimization procedure. The transition to the SM
phase is more subtle as it is dominated by sub-leading corrections at small
$\ell$. In this regime, the optimal profile and the nature of this phase
transition might be affected by the possibility of MIPS: this could certainly
be investigated further.

To conclude, in this paper we have calculated the large deviation properties
of the entropy production rate in an active gas model. The special properties
of this model (in which the orientational dynamics share the hydrodynamic
timescales of diffusion and self-propulsion) allow such calculations to be
done exactly. Despite these special features, the result show many
similarities with what is known numerically from models in which the
orientational dynamics represent a fast field rather than a hydrodynamic one.
The resulting phase diagram shows considerable complexity and features phase
separated, travelling band, and smoothly modulated states. It also shows
regions of time-like phase separation in which the system spends part of its
trajectory in a modulated state and part of it in a homogeneous one. Many
interesting features of the phase diagram, including but not limited to
tricritical behaviour, stem from the fact that the biasing variable (in this
case the entropy production rate) is not everywhere a convex function of
density. We expect similar features in other systems subject to bias with that
property; this will be explored further in Agranov  _et al._ (2022).

###### Acknowledgements.

We thank Jakub Dolezal for sharing with us his code for functional
minimization and Sunghan Ro for helping prepare Fig.1. TA is funded by
Blavatnik Postdoctoral Fellowship Programme. MEC is funded by the Royal
Society. Work funded in part by the European Research Council under the
Horizon 2020 Programme, ERC Grant Agreement No, 740269.

##  Appendix A The microscopic EPR

Our aim in this section is to compute the microscopic EPR and then to express
it in terms of hydrodynamic fields.

###  A.1 Preliminary: defining microscopic flows

The microscopic calculation of the EPR involves counting all microscopic
dynamical moves that contribute to the entropy production. To this aim we now
define three counting (or flow) variables:

  * •

Total flow: $N_{+}^{T}$ is the difference between the total numbers of right
and left hops performed by $+$ particles up to time $T$. The corresponding
quantity for $-$ particles is $N_{-}^{T}$.

  * •

Flow into vacant sites: $N_{+,0}^{T}$ is the difference between the total
numbers of right and left hops into vacant sites performed by the $+$
particles up to time $T$. Similarly $N_{-,0}^{T}$ for $-$ particles.

  * •

Flow due to swaps: $N_{+,-}^{T}$ is the difference between the total numbers
of right-swaps and left-swaps between a $+$ and a $-$ particles up to time
$T$, where a right-swap is a swap move in which the $+$ particle moves to the
right.

These three variables are trivially related by

$N_{+}^{T}=N_{+,0}^{T}+N_{+,-}^{T}\quad,\quad
N_{-}^{T}=N_{-,0}^{T}-N_{+,-}^{T}.$ (68)

###  A.2 Computing the microscopic EPR

Out of the possible dynamical moves of the system, the only one which is not
reversible and entails an entropy gain is when a particle swaps positions with
an adjacent vacant site. For a $+$ particle hoping to a right neighboring
vacant site this move contributes an entropy gain of

$\Delta\mathbb{S}=\log\left(1+\frac{\lambda}{DL}\right)=\frac{\lambda}{DL}+\mathcal{O}(L^{-2}).$
(69)

The reverse move (a $+$ particle swapping position back with a left vacant
site) contributes the opposite amount $-\Delta\mathbb{S}$. Thus, the total
entropy gain up to some time $T$ due to the $+$ particles dynamics is given by
the flow into vacant sites $N_{+,0}^{T}$ (and similarly for $-$ particles).
All in all, the total entropy gain up to time $T$ is

$T\mathbb{S}_{T}=(N_{+,0}^{T}-N_{-,0}^{T})\Delta\mathbb{S}.$ (70)

We must now express the flows into vacant sites $N_{+(-),0}^{T}$ in terms of
hydrodynamic fields. We will do so by first expressing the total flows
$N_{+(-)}$ and total swaps $N_{+,-}$ in terms of hydrodynamic fields, and then
using (68) to derive $\mathbb{S}_{T}$.

###  A.3 Expressing the flows in terms of hydrodynamic fields

The hydrodynamic currents $J_{+}(x,t)$ and $J_{-}(x,t)$ are directly related
to the microscopic dynamics: for a short time $\Delta t$ (measured on the
hydrodynamic scale), $L\Delta tJ_{+}(x,t)$ is the number difference between
right and left hops of $+$ particles across a single edge of the lattice, at
position $x$. Summing over all edges on the lattice and integrating over time
gives the total flow

$N_{+}^{T}=L^{2}\int_{0}^{T}dt\int_{0}^{1}dxJ_{+}(x,t)\;,$ (71)

A similar expression holds for $-$ particles.

The next step is to relate the number of swaps $N_{+,-}$ to the currents. In
the absence of tumbling, the existence of counterpropagating currents
determines how many swaps must take place: at leading order in $L$ one has

$N_{+,-}^{T}=L^{2}\int_{0}^{T}dt\int_{0}^{1}dx\left[J_{+}(x,t)\rho_{-}(x,t)-J_{-}(x,t)\rho_{+}(x,t)\right],$
(72)

We will derive this result below, showing at the same time that it remains
valid in the presence of tumbling. (Inutuitively, the number of swap events
between each tumble is large enough that the result without tumbling still
applies.)

The scaling of the microscopic rates with $L$ (see Sec. II) means that we can
define a mesoscopic time scale $\Delta\hat{T}=L^{1+\delta}$ with $0<\delta<1$
(here $\Delta\hat{T}$ is measured in microscopic time units). Over such a time
scale, the typical displacement of a particular particle during
$\Delta\hat{T}$ will scale as $\sim\lambda L^{\delta}$ on account of the
biased hops. This defines a mesoscopic length scale corresponding to a large
number of lattice sites, so a large number of swaps must occur in this time.
(Still, this scale is much smaller then the system size $L$.) It is on such
mesoscopic length compartments that the local values of the hydrodynamic
fields are defined, see Eq.(1). In the following we will assume the large $L$
limit and express local sample means in terms of the appropriate hydrodynamic
fields.

It is important that a very small fraction of particles will tumble within the
mesoscopic time scale. In fact, the probability that a finite fraction of
particles will tumble is even beyond the large deviation scaling for the
tumbling rate probability cost (14), and hence can be neglected altogether,
see also Agranov  _et al._ (2022). Thus, for what follows, we can neglect
tumbling during $\Delta\hat{T}$.

Hence we can define an individual displacement of a $+$ particle $i$, which we
denote by $n_{+}^{i}$ (here the superscript $\Delta\hat{T}$ is suppressed for
ease of notation). We relate this quantity to hydrodynamic fields. Consider a
mesoscopic compartment $\Delta x$ centered at $x$ ($x\in[0,1)$ being a
macroscopic variable). Now define the sample mean of particle’s displacement
during $\Delta\hat{T}$ over all particles in such a compartment $n_{+}(x,t)$.
Then the contribution to the total flow $N_{+}^{\Delta\hat{T}}$ coming from
that compartment is expressed in terms of the local hydrodynamic fields via

$\sum_{i\in\Delta x}n_{+}^{i}=L\Delta x\rho_{+}(x,t)n_{+}(x,t)=L\Delta
xJ_{+}(x,t)\frac{\Delta\hat{T}}{L}.$ (73)

Notice that here we use microscopic time. This is why there is an $L^{2}$
scaling difference in the right hand side, compared to (71). From (73) we have
that

$\rho_{+}(x,t)n_{+}(x,t)=J_{+}(x,t)\frac{\Delta\hat{T}}{L}.$ (74)

Next, we can also relate the individual flows $n_{+}^{i}$ to swaps. Focus on a
specific particle $i$ and denote by $\Delta_{i,j}$ its relative displacement
from a specific $-$ particle labeled $j$ at time $\hat{t}=0$. Then the
contribution to the total swaps $N_{+,-}^{\Delta\hat{T}}$ which is coming only
from the particle $i$ is $N_{+,-}^{i}$, given by

$N_{+,-}^{i}=\sum_{j}\Big{\lfloor}\frac{n_{+}^{i}-n_{-}^{j}}{\Delta_{i,j}}\Big{\rfloor},$
(75)

where $\lfloor\quad\rfloor$ is the floor function. Now, since
$n_{+}^{i},n_{-}^{j}\sim\lambda L^{\delta}$ are mesoscopic, all relevant $j$
particles that contribute to this sum are within the mesoscopic compartment
$\Delta_{i,j}<\Delta x\sim L^{\delta}$ around particle $i$. Contribution from
particles outside of this compartment is extremely rare and is beyond the
large deviation scaling we consider here. In particular, we can assume the
particles did not make a complete revolution around the system.

On the other hand, as this compartment contains a large (mesoscopic) number of
$-$ particles, then their individual flows $n_{-}^{j}$ are self averaging and
can be replaced by the sample mean $n_{-}(x,t)$. Thus, we have that the sum in
(75) has contributions from all $-$ particles at distance
$n_{+}^{i}-n_{-}(x,t)$ from particle $i$; the total number of such particles
is $(n_{+}^{i}-n_{-}(x,t))\rho_{-}(x,t)$, and so

$N_{+,-}^{i}=\left[n_{+}^{i}-n_{-}(x,t)\right]\rho_{-}(x,t)\;.$ (76)

Summing over all $+$ particles in the system we have that

$N_{+,-}^{\Delta\hat{T}}=\sum_{i}N_{+,-}^{i}=L\int_{0}^{1}dx\left[n_{+}(x,t)-n_{-}(x,t)\right]\rho_{-}(x,t)\rho_{+}(x,t).$
(77)

Using the relation (74), switching to the hydrodynamic time scale, and
integrating over time, we finally arrive at (72).

###  A.4 Expressing $\mathbb{S}_{T}$ in terms of the hydrodynamic fields

We are left with substituting the two results (71) and (72) into (68) to
arrive at the expression for the flow into vacant sites in terms of
hydrodynamic fields to obtain

$N_{+(-),0}=L^{2}\int_{0}^{T}dt\int_{0}^{1}dx\left[J_{+(-)}(x,t)(1-\rho_{-(+)}(x,t))+J_{-(+)}(x,t)\rho_{+(-)}(x,t)\right].$
(78)

Finally, using this expression in (70) we arrive the at the microscopic EPR
expressed in terms of hydrodynamic fields

$\displaystyle\mathbb{S}_{T}$ $\displaystyle=$
$\displaystyle\frac{L\gamma\text{Pe}^{2}}{T}\,\int_{0}^{T}dt\int_{0}^{1}dx\left[(1-2\rho_{-}(x,t))\frac{J_{+}(x,t)}{\lambda}-(1-2\rho_{+}(x,t))\frac{J_{-}(x,t)}{\lambda}\right]$
(79) $\displaystyle=$
$\displaystyle\frac{L\gamma\text{Pe}^{2}}{T}\,\int_{0}^{T}dt\int_{0}^{1}dx\left[(1-\rho(x,t))\frac{J_{m}(x,t)}{\lambda}+m(x,t)\frac{J_{\rho}(x,t)}{\lambda}\right],$
(80)

in exact agreement with the expression (23) that was obtained from the coarse
grained fluctuating hydrodynamics. In particular, the expected microscopic EPR
coincides with the macroscopic IEPR (26).

###  A.5 Insights from the microscopic derivation

The microscopic derivation helps explain two features of the IEPR. First, the
density dependence of the expected IEPR $\bar{a}(\rho_{0})$ (26) can be
understood. Indeed, for the expected IEPR one has
$\rho_{+}=\rho_{-}=\rho_{0}/2$, $J_{+}=-J_{-}=\text{const}$ and
$N_{+}^{T}=-N_{-}^{T}=L^{2}J_{+}$. Plugging these into (72) one finds that the
total swaps are given by $N_{+,-}^{T}=N_{+}^{T}\rho_{0}$. Using this in (68)
we find that the flow into vacant sites is a fraction $(1-\rho_{0})$ of the
total flow

$N_{+,0}^{T}=(1-\rho_{0})N_{+}^{T}\quad;\quad
N_{-,0}^{T}=(1-\rho_{0})N_{-}^{T}.$ (81)

As we detail in Sec. III, the factor $(1-\rho_{0})$ helps explain the density
dependence in (26).

Secondly, the transition to collective motion is also explained by the
microscopic picture. By increasing the plus minus imbalance $m$, the number of
reversible swaps is reduced, while leaving the sum of total currents
unchanged. By decreasing the portion of reversible moves out of the total
dynamical moves that comprise the integrated current, the entropy production
is increased.

Indeed, assuming stationary and homogeneous density profiles and using the
expressions (II.2) for the expected fluxes, one finds that the expected total
swap (72) is given by

$N_{+,-}^{T}=\frac{L^{2}\lambda T(1-\rho_{0})}{2}(\rho_{0}^{2}-m^{2}).$ (82)

This expression monotonically decreases with $|m|$ until it vanishes for the
completely aligned case $|m|=\rho_{0}$. At the same time, the sum of total
flows in the direction of the self propulsion is independent of $m$

$N_{+}^{T}-N_{-}^{T}=L^{2}T\lambda\rho_{0}(1-\rho_{0})$ (83)

##  Appendix B Deriving the square gradient expansion (40)

This appendix derives the expansion (40) from the minimization problem (28),
which gives an expression for the rate function $I$ up to quadratic order in
$\ell$. As a preliminary step, we first derive a similar expression which is
valid at zeroth order in $\ell$, namely

$\lim_{\ell\to
0}I(a)=\min_{\alpha(x),\rho(x)}\int_{0}^{1}dxI_{H}\big{(}\rho(x),\alpha(x)\big{)},$
(84)

where the minimisation is subject to the usual constraints
$a=\int_{0}^{1}\alpha(x)dx$, and $\rho_{0}=\int_{0}^{1}\rho(x)dx$. (Recall
that the rate function $I$ depends implicitly on $\rho_{0},\ell$ as well as
$a$.) Building on this analysis, the full result (40) is then derived in
Appendix B.2. Since (84) does not involve any gradient terms, the infimum is
given by the convex hull of $I_{H}$. We will show that this is true for all
$\rho_{0},a$, even if the minimizer of (28) is non-stationary, as in the TB
phase.

###  B.1 Derivation of (84)

Let us first write the Lagrangian ${\cal L}_{J}$ in terms of dimensionless
quantities, rescaling the currents $J_{\rho}$ and $J_{m}$ by factors of
$\lambda$:

$\mathcal{L}_{J}=\frac{\gamma}{2}\begin{bmatrix}\text{Pe}\,\frac{J_{\rho}}{\lambda}-\text{Pe}\,m(1-\rho)+\ell\partial_{x}\rho\\\
\text{Pe}\,\frac{J_{m}}{\lambda}-\text{Pe}\,\rho(1-\rho)+\ell\partial_{x}m\end{bmatrix}^{\dagger}D\mathbb{C}(\rho,m)^{-1}\begin{bmatrix}\text{Pe}\,\frac{J_{\rho}}{\lambda}-\text{Pe}\,m(1-\rho)+\ell\partial_{x}\rho\\\
\text{Pe}\,\frac{J_{m}}{\lambda}-\text{Pe}\,\rho(1-\rho)+\ell\partial_{x}m\end{bmatrix}.$
(85)

Since all gradient terms appear with a prefactor $\ell$, they disappear as
$\ell\to 0$ for the computation of the action (28). [This is still true even
when the optimal fields are phase separated with a domain wall of width
$\sim\ell$. The gradient terms are $\mathcal{O}(1)$ in that case but the
region of integration where they contribute is $\mathcal{O}(\ell)$, and hence
their overall contribution to the functional is $\mathcal{O}(\ell)$, and
therefore negligible.]

Next, write the dynamical equations (4) and (5) as

$\displaystyle\gamma^{-1}\partial_{t}\rho$ $\displaystyle=$
$\displaystyle-\text{Pe}\,\ell\partial_{x}(J_{\rho}/\lambda)$ (86)
$\displaystyle\gamma^{-1}\partial_{t}m$ $\displaystyle=$
$\displaystyle-\text{Pe}\,\ell\partial_{x}(J_{m}/\lambda)-2K.$ (87)

As $\ell\to 0$, Eq. (86) shows that either $\rho$ must be stationary, or the
gradient of $J_{\rho}$ must diverge. The latter scenario can occur at sharp
interfaces, leading to a TB state where the density is stationary within the
bulk phases, with currents at the interfaces which cause the band to travel
(see Sec. VII).

The same property holds for the optimal magnetization history $m$. That is, it
must also be either stationary or take the form of a traveling band. Here this
does not follow immediately from (87), and the proof is slightly more
involved. One way of arriving at it is employing the Hamiltonian MFT
formulation of the problem and using it to prove stationarity of $m$. Then
from Eq.(87), one has $K=0$ everywhere except (possibly) in the vicinity of
sharp interfaces. Recall, $K=0$ is the same value that one finds for the
homogeneous solution, see Sec. V.

Setting $K=0$, and neglecting the gradient terms in (85), the action
functional (11) is of the form considered in Sec. V, except that the fields
$\rho,m,J_{\rho},J_{m}$ entering the action functional are not constants.