Stochastic Hydrodynamics of Complex Fluids: Discretisation and Entropy Production

Many complex fluids can be described by continuum hydrodynamic field equations, to which noise must be added in order to capture thermal fluctuations. In almost all cases, the resulting coarse-grained stochastic partial differential equations carry a short-scale cutoff, which is also reflected in numerical discretisation schemes. We draw together our recent findings concerning the construction of such schemes and the interpretation of their continuum limits, focusing, for simplicity, on models with a purely diffusive scalar field, such as ‘Model B’ which describes phase separation in binary fluid mixtures. We address the requirement that the steady-state entropy production rate (EPR) must vanish for any stochastic hydrodynamic model in a thermal equilibrium. Only if this is achieved can the given discretisation scheme be relied upon to correctly calculate the nonvanishing EPR for ‘active field theories’ in which new terms are deliberately added to the fluctuating hydrodynamic equations that break detailed balance. To compute the correct probabilities of forward and time-reversed paths (whose ratio determines the EPR), we must make a careful treatment of so-called ‘spurious drift’ and other closely related terms that depend on the discretisation scheme. We show that such subtleties can arise not only in the temporal discretisation (as is well documented for stochastic ODEs with multiplicative noise) but also from spatial discretisation, even when noise is additive, as most active field theories assume. We then review how such noise can become multiplicative via off-diagonal couplings to additional fields that thermodynamically encode the underlying chemical processes responsible for activity. In this case, the spurious drift terms need careful accounting, not just to evaluate correctly the EPR but also to numerically implement the Langevin dynamics itself.


Introduction
Numerous complex fluid systems can be described by continuum equations formulated at the hydrodynamic level. This reflects the fact that their important structure and dynamics arises at a mesoscopic scale not a molecular one. Examples include theories of flowing liquid crystals described by vector or tensor order parameters [1,2], and those of partially miscible binary fluid mixtures, described by a conserved scalar composition variable [3]. The latter can undergo phase separation via a combination of diffusive motion and fluid flow, for which the canonical model is called Model H in the classification of Hohenberg and Halperin [4]. An important special case of Model H, in which the fluid velocity is set to zero so that phase separation proceeds by diffusion only, is called Model B. The latter describes various physical processes in complex fluids, such as Ostwald ripening of emulsion droplets, where the coupling between diffusion and fluid flow is unimportant.
These hydrodynamic-level descriptions are often first encountered as deterministic equations of motion. This is sometimes sufficient, for example, in Ostwald ripening of emulsions where large droplets grow at the expense of small ones via deterministic diffusive fluxes. However, there are many other processes in binary fluids (and also liquid crystals), ranging from droplet nucleation to dynamics near critical points, where the stochasticity of the continuum models must be retained so as to maintain a faithful description of thermal fluctuations. Note that this is even true of single-phase fluids whose true quiescent state involves a Boltzmann distribution for the velocity field v(r), not the state of zero velocity predicted by the Navier-Stokes equation in the absence of forcing. As first shown by Landau and Lifshitz, this is fixed by adding a fluctuating thermal stress to the Navier-Stokes equation [5]. The resulting thermal fluctuations in the fluid then impart Brownian motion to any colloidal particle suspended in it, without the need for a separate Langevin force on the colloid.
In the hydrodynamic modelling of complex fluids, it is therefore important to be able to handle thermal noise terms correctly, both at a conceptual level in the continuum and when creating discrete implementations of the continuum equations for use in computer simulation studies. The first of these tasks poses technical challenges of surprising complexity, which can only be resolved by studying the discretisation issue. The reason for this is simple: adding noise converts the PDEs of deterministic complex fluid models into Stochastic PDEs (SPDEs), which, in general, have no mathematical meaning without some sort of cutoff at short scales. (In a few favourable cases, meaning has been restored directly at the continuum level by a procedure that effectively constructs the renormalization group and the continuum limit simultaneously [6].) In terms of physical modelling, the existence of a cutoff is unproblematic: continuum descriptions, such as the Beris-Edwards equations for liquid crystals or Models H and B for binary fluids, only hold at scales larger than the molecular one. Mathematically, however, once noise is included, the cutoff can infiltrate the continuum models in unexpected ways. For example, we will find below that trying to work directly in the continuum limit gives in the equations under study undefined mathematical objects, such as δ(0)-the Dirac delta-function evaluated at zero argument. This is symptomatic of a quantity that diverges as the cutoff becomes small. Moreover, we know from equilibrium statistical physics that a particular quantity of interest may or may not depend on the cutoff according to details of the model. For example, if a scalar-order parameter field has Gaussian fluctuations at wavenumber q, |φ q | 2 = G −1 (q), then the corresponding real-space variance |φ(r)| 2 either remains finite or blows up with the cutoff according to the convergence at high q of G −1 (q)dq. This real-space variance is a legitimate object of enquiry. However, hydrodynamic descriptions such as Model H and B effectively expand G as a low-order polynomial in q on the basis that the high q behaviour is not important. For this reason, it is unwise to assume that the continuum limit of the models studied by physicists always make sense.
Turning from that conceptual issue to the more practical one of numerically discretising the hydrodynamic equations of a thermal complex fluid, there emerges a crucial requirement for the treatment of noise that creates further surprising traps for the unwary. This is the requirement that the discretised equations respect the principle of detailed balance. Put differently, if one sets up a numerical model for a complex fluid and calculates its entropy production rate (EPR) in a steady state of thermal equilibrium, the EPR should vanish. We will see below that there are various different ways in which numerical analyses can fail this test.
One setting in which the issue of entropy production comes to the fore is in the study of active field theories [7]. These are stochastic hydrodynamic models intended to describe active complex fluids whose microscopic components are driven by an internal power supply. Examples of such active fluids include suspensions of motile bacteria and of autophoretic colloidal particles with asymmetric surface chemistry that catalyses a chemical reaction, creating chemical gradients that drive the colloids forward. The study of active matter has exploded into a field whose detailed discussion would take us far beyond the topic of this paper; see [8]. For the present purposes, we can regard active field theories as extensions of the stochastic hydrodynamic equations for complex fluids in which detailed balance is deliberately broken by the inclusion of new terms that do this, usually at the lowest possible order in the expansion in order parameter fields and their gradients.
A strategy we have recently developed in studying such active field theories is to quantify their mesoscopic irreversibility by calculating the steady-state EPR directly at the level of the fluctuating order parameter field dynamics [9][10][11][12]. This quantity is best-called the informatic EPR or IEPR [13]: it makes no attempt to capture all the microscopic irreversibility or heat flows associated with the particle motions underlying the coarse-grained, hydrodynamic SPDEs. Instead, the IEPR is computed informatically from forward and reverse path probabilities using the tools of stochastic thermodynamics [14] applied to the SPDEs themselves. These tools have also found applications in active matter systems such as biochemical signalling [15], mechanosensory processes [16] and bacterial motion [17]. We have further extended these ideas and embedded a large class of active field theories in a thermodynamically consistent setting that accounts for their driving mechanism, in which case, the irreversibility of the enlarged system capture the actual rate of heat production. In our studies of active field theories, we have found interesting physics to be laid bare when one considers the way the IEPR (and the heat rate) depends on the spatial configuration of the system and also the way different contributions to it (e.g., bulk or interfacial) scale with the noise level. To address these issues by computer simulation, it is clearly crucial to have a numerical implementation in which the calculated entropy production arises solely by virtue of the active, detailed-balance-breaking terms, unpolluted by any failure of the numerical discretisation scheme to respect detailed balance even in thermal equilibrium.
Accordingly, in our recent studies of active field theories, we have been forced to carefully consider the conceptual and discretisation issues for the stochastic hydrodynamics of complex fluids generally. We have found that, beyond a few important contributions such as [18,19], these issues are not widely discussed in the literature accessible to physicists-especially not in relation to entropy production and its numerical evaluation. Thus far, our own results on these topics have been presented only incidentally, if at all, in technical appendices and side remarks in papers on how active hydrodynamic models actually behave [9][10][11][12]. We attempt here a coherent perspective on these issues. For simplicity, our main focus is on Model B and its active counterparts, in which the sole order parameter is a scalar field and the only dynamics is diffusive. Indeed, between here and the concluding section, we say nothing of the wider class of complex fluid models containing vector and tensor order parameters (for liquid crystals) or even a coupling to fluid flow (for a scalar field, Model H). We emphasize, however, that the conceptual and discretisation issues addressed here apply, in varying degrees, to all these other cases.
The rest of this paper is structured as follows. To set the stage, Section 2 reviews the questions of discretisation and spurious drift for a single particle Langevin equation with multiplicative noise, discussing also the Fokker-Planck equation, path integrals, and entropy production in this simplified setting before addressing the stochastic calculus for finitely many degrees of freedom. This establishes a core set of ideas that are utilized subsequently for the continuum case. In Section 3, we turn to continuous fields, focusing on the case of (active) Model B where the noise is additive rather than multiplicative, and show how the spatial discretisation must be carefully handled to avoid erroneous evaluation of the (informatic) entropy production. We focus on finite difference schemes, as opposed to spectral ones, for spatial discretisation because, besides being widely used, this approach offers the most direct way to illuminate problems with the continuum limit. In Section 4, we consider how to embed an active field theory within an enlarged description that is thermodynamically consistent in the sense that it accounts for heat flow (caused, in this instance, by chemical reactions that drive the system microscopically) at the level of linear irreversible thermodynamics. We review how this generically leads to multiplicative noise even where none was previously present and describe the further conceptual and discretisation problems arising from this. Finally, in Section 5, we offer some brief concluding remarks.

Stochastic Thermodynamics of Particles
In this section, we establish some basic concepts concerning stochastic differential equations and thermal motion, starting in the context of a single particle and then turning to the case with several degrees of freedom.

Langevin Equation
Let us consider a colloidal particle suspended inside a viscous solvent in one dimension. The solvent acts as a heat bath for the particle with temperature T, and the particle is assumed to be in thermal equilibrium with the heat bath at all times. Let us denote x(t) the stochastic trajectory of the centre of mass of the particle. The equation of motion for the particle is then given by the overdamped Langevin equation: where U(x) is an external potential (provided, e.g., by an optical trap), and η(t) is a Gaussian white noise with zero mean η(t) = 0 and (Dirac) delta-function correlation Note that x(t) is a stationary process since the potential U does not explicitly depend on time. In Equation (1), we neglect the inertia of the particle, which is valid if the Reynolds number is much smaller than unity; Γ(x) is the mobility or the inverse of the friction coefficient. (For a spherical particle, Γ = 1/(6πηR), where η is the viscosity of the solvent and R is the radius of the particle.) In this example, we also allow the mobility Γ(x) to vary locally in space. D(x) in Equation (1) is the diffusion coefficient or the noise strength. The noise strength D(x), the mobility Γ(x), and the solvent temperature T are all related via the Stokes-Einstein relation, which is a direct consequence of the fluctuation-dissipation theorem (FDT): D(x) = Γ(x)T (note that we work in units of k B = 1). Since the mobility, and hence the diffusion constant, vary locally in space, the noise in Equation (1) is multiplicative. Since the noise is multiplicative, the Langevin equation as written in Equation (1) is ambiguous, unless we specify how we discretise the dynamics in time. Depending on how we do so, we may need to include the spurious drift term ν a (x) in Equation (1) to recover Boltzmann distribution in the steady state. The 'spurious drift' terminology is conventional but may be confusing: the term ν a (x) arises in effect because the noise, depending on the discretisation scheme used, might or might not still have zero average. (For instance, if the noise variance increases with x and is evaluated mid-step, then random steps in the positive x direction are larger than those towards negative x.) Finally, to simplify the notation, we shall also define:

Discretised Langevin Equation
Let us discretise the time into t n = t 0 + n∆t, where n = 0, 1, 2, . . . , N. The Dirac delta function correlation in the continuous noise η(t)η(t ) = δ(t − t ) can be regularized into a Kronecker delta η m η n = δ mn /∆t. The discretised Langevin equation is then given by: where {ξ 0 , ξ 1 , . . . ξ N−1 } are a set of independent Gaussian random variables with zero mean, ξ n = 0, and Kronecker delta correlation ξ m ξ n = δ mn . In (4), a ∈ [0, 1] is the discretisation parameter, which tells us when during the timestep we should evaluate the particle position x for the purpose of sampling the noise (whose variance is, we recall, xdependent). Thus, a = 0 corresponds to the Itô choice (initial postion), a = 1/2 corresponds to Stratonovich (midpoint position), and a = 1 corresponds to anti-Itô (final position). Now, using the mean value theorem x n+a = x n + a∆x n , we write Equation (4) as: In order to derive the Fokker-Planck equation below, we first need to compute the first and second moment of ∆x n :

Fokker-Planck Equation
Let us denote P(x, t|x 0 , t 0 )dx to be the probability of finding the particle at [x, x + dx] at time t, given that it was at x 0 at time t 0 , where t 0 < t. The time evolution of this probability density function is given by Kramers-Moyal expansion (see [20] for derivation): where ∆x(t) = x(t + ∆t) − x(t). Substituting Equations (6) and (7) into the equations above and taking the limit ∆t → 0, we obtain: We can also write this as a continuity equation ∂P/∂t = −∂J/∂x, where the probability current is given by: For an equilibrium system, which is the case in our example, the probability current should be equal to [18,21]: Together with FDT, D(x) = Γ(x)T, the probability current from Equation (11) will guarantee Boltzmann distribution in the steady state: P(x, t → ∞) ∝ e −U(x)/T . Comparing Equation (11) to Equation (10), we thus require the spurious drift to be In the case of additive noise, where D and Γ are constant, the spurious drift is always zero. In the case of multiplicative noise, where D(x) and Γ(x) vary in space, the spurious drift can be made to vanish only by choosing the anti-Itô discretisation a = 1. Generally speaking, numerical strategies are simplest for Itô (a = 0), whose update statistics depend only on the state at the start of the timestep in which the update is to occur. On the other hand, the Stratonovich discretisation (a = 1/2) has some desirable properties in relation to temporal reversibility, see Section 2.5. Moreover, as we will see in Equation (38) below, setting a = 1 does not eliminate all spurious drift terms in higher dimensions, where such terms remain a generally unavoidable feature.

Path Integral Formalism
The Fokker-Planck equation in Equation (9) is usually rather difficult to solve when generalising to higher dimensions. In many situations (e.g., when calculating the entropy production rate), it is often easier to work with the path probability.

Transition Probability
Suppose that our particle is initially at x n at time t n . For a given noise realization ξ n , the position of the particle x n+1 in the next timestep t n+1 is given by the discretised Langevin Equation (4), in which ξ n is a Gaussian random variable with probability density function: We can then substitute Equation (4) into Equation (13) to obtain the probability of finding the particle at [x n+1 , x n+1 + dx n+1 ] at time t n+1 , given that it was at x n at the previous timestep t n : Note that the Jacobian |dξ n /dx n+1 | is inserted when we change the random variable from ξ n to x n+1 . To find the Jacobian, we first express ξ n as a function of x n+1 from Equation (4) to obtain: where we have used the mean value theorem x n+a = x n + a(x n+1 − x n ) again. Taking a derivative with respect to x n+1 , we then obtain: We want to exponentiate the terms inside the curly bracket; however, we note that x n+1 − x n ∼ √ ∆t, so we cannot conduct this directly. Instead, we shall use the following Taylor expansion [19], for some constant C ∼ ∆t 0 : Here, we have approximated ∆x n ∆x n g 2 (x n+a )∆t, which is valid for small ∆t, c.f. Equation (7) and [19]. The Jacobian can then be exponentiated as follows: where f , g, f , and g are evaluated at x n+a . Substituting Equation (19) into Equation (14), we obtain: where f , g, f , and g are again evaluated at x n+a . Finally after completing the square, we obtain the transition probability: Below, we will often use a shorthand notation whereby the O(∆t 3/2 ) term is implicit in expressions such as this.

Path Integral
Suppose that, initially, the particle is at x 0 at time t 0 . What is the probability that we find the particle at [x N , x N + dx N ] at time t N ? This probability can be written as (Chapman-Kolmogorov equation): Substituting the transition probability from Equation (21) into the equation above, we obtain: where A{x n } is called the dynamical action (Onsager-Machlup action) [18,19]: (24) and N {x n } is some normalization prefactor, which is constant for additive noise: In the limit ∆t → 0, Equation (23) becomes a path integral, i.e., we sum over all possible trajectories {x n |n = 0, 1, 2, . . . , N}, each with a weight or path probability P {x n } = N {x n } e −A{x n } . Here, f (x) = −Γ(x)U (x) + (1 − a)D (x) as in Equations (1) and (12). Note that the expressions for N and A in Equations (24) and (25) are generic for any processes. When we invoke FDT below, D(x) = Γ(x)T, we then assume {x n } to be a stationary process. Furthermore, note that for additive noise, where D and Γ are constant, the dynamical action still depends on the discretisation parameter a, even though the discretised Langevin dynamics does not depend on a anymore. This is important when calculating entropy production via the path probabilities, as we consider next.

Entropy Production
Consider now a single stochastic trajectory of our overdamped particle, {x n |n = 0, 1, 2 . . . N}, which was generated by the discretised Langevin Equation (4). A foundational result of stochastic thermodynamics [14,22] is that the total heat dissipated from the particle to the environment, as it moves along this single trajectory, obeys: The second equality states that the heat dissipation determines the increase in the entropy of the medium or heat bath supplying the noise, ∆S m . In (26), P {x n } is the probability of obtaining this particular trajectory {x n |n = 0, 1, 2, . . . N} and P R {x n } is the probability of observing the time-reversed trajectory {x R n = x N−n |n = 0, 1, 2, . . . N} (see Figure 1) under the same Langevin dynamics (4). For example, the chosen trajectory might be a particle going from high to low energy, in which case the time-reversed trajectory is much less probable to be seen under the same forward Langevin dynamics (4). Note that to keep the same notation as in the previous literature [12,13,23], P R {x n } := P {x R n } in Equation (26). Furthermore, note that if the specific trajectory {x n } results from a specific protocol, such as changing one of the parameters inside the potential energy U(x) (which we do not consider in this paper), in the time-reversed trajectory, we also have to reverse the direction of this protocol [14,16,24,25].

Evaluation via Discretised Action
For the forward trajectory {x n }, the path probability obeys P {x n } = N {x n } e −A{x n } , where N {x n } and A{x n } are given in Equations (24) and (25), respectively. For the timereversed trajectory {x R n }, the path probability is given by P {x R n } = N {x R n } e −A{x R n } , where N and A are still the same expressions, given in Equations (24) and (25), except that we replace the arguments by {x R n }. Let us first calculate the normalization prefactor N {x R n } for the reversed trajectory {x R n }: In the second equality above, we have substituted x R n = x N−n . Comparing Equation (25) to Equation (27), the normalization factor N differs between forward and backward discretised paths and does not cancel in Equation (26), unless we choose the Stratonovich discretisation, a = 1/2. This is a compelling reason to choose Stratonovich when calculating the heat dissipation or entropy production via Equation (26), and we do so hereafter. In Stratonovich, the action for the reversed trajectory {x R n } is then given by: The heat dissipated is simply the difference between the backward and the forward action: Finally, we apply FDT where ∆U = U(t N ) − U(t 0 ), and since this is a Stratonovich integral, we have used the standard chain rule in the last equality. Thus, we recover the first law of thermodynamics. The Stratonovich integral over U (x) dx in Equation (31) can, if desired, be converted to an Itô integral by setting U (x n+ 1 2 ) = U (x n ) + 1 2 U (x n )∆x n to obtain: Finally, substituting Equation (31) back to Equation (26), we may also show that detailed balance is obeyed: as is indeed required for any system in thermal equilibrium.

Non-Equilibrium Steady State
We may generalize the above result to non-equilibrium steady states. For example, we may imagine applying a non-conservative force F(x) on the particle so that the Langevin equation now reads: where D(x) = Γ(x)T. In the case of a periodic potential U(x), a constant external force F may give rise to a steady-state current, which indicates a non-equilibrium steady state and thus breaks detailed balance. The EPR in the steady-state ensemble is found as follows: where P [x(t )] is the path probability for some forward trajectory {x(t )|t ∈ [0, t]} and P R [x(t )] is the path probability for the same trajectory going backwards in time {x(t − t )|t ∈ [0, t]}. The angle bracket indicates ensemble averaging or the average over different noise realizations {η(t )|t ∈ [0, t]}. Note that, since the entropy content of the system is unchanging in the steady state, all the entropy produced within it ends up in the medium or heat bath so that the EPR, which is the rate of change of entropy in the bathṠ, equates to the dissipation rateQ within a factor T. The notation in (35) is chosen to connect with subsequent sections and with the previous literature; note, however, that in [9,10,12], the un-accented symbol S is used to denotes the entropy production rate, which is calledṠ in this paper. Following the same derivation as above, we can show that the steady-state heat production rate is:Q Note that dU/dt is zero in the steady state, on average. Thus, the rate of heat dissipation is equal to the average rate of work conducted by the external force F, again consistently with the first law.
The above results are given in thermodynamic language which ultimately rests on the first law (conservation of energy). However, Equation (34) is not generically thermodynamically consistent [26,27]. For example, one could interpret Equation (34) as describing an active particle, such as a swimming microorganism, for which the term is a spatially varying propulsive velocity. The x dependence of V might then have no connection with energetics (that is, F(x) is no longer a mechanical force), reflecting instead a tendency to swim in the positive or negative direction depending on external stimuli, such as an imposed gradient in nutrient levels (for instance, ΓF ∝ ∂ x H(x) with H a food concentration) [28]. In such cases, there is no first law behind Equation (34), and we cannot associate ln(P /P R ) in Equation (35) with energy dissipation or heat production. We can nonetheless define an informatic entropy production rate, or IEPR, via the first equality only in Equation (35). It is this IEPR that we will generalize in Section 3 as a tool for quantifying the irreversibility of active field theories. Thereafter, in Section 4, we will restore a link with thermodynamics and the first law, under specific assumptions concerning the near-equilibrium character of the microscopic dynamics responsible for activity.

Stochastic Calculus for d > 1 Degrees of Freedom
Let us consider the general Langevin equation for d > 1 degrees of freedom in a system with detailed balance. We denote the coordinates to be x i , where i = 1, 2, . . . , d. (This could describe either one particle in d dimensions, or N > 1 particles in d/N dimensions.) The Langevin equation for {x i (t)} is given by where {η i (t)} are Gaussian white noises with zero mean η i (t) = 0 and delta-correlations . The deterministic part f i and the noise prefactor g ij can be written as [18]: Here, U({x i }) is the potential energy, and a ∈ [0, 1] is the time discretisation parameter as before. The second and the third terms in Equation (38) constitute the spurious drift ν i , whose presence ensures the Boltzmann distribution in the steady state: (38) Γ ij and D ij are the mobility and the diffusion matrix, respectively, which must satisfy FDT: Onsager symmetry requires Γ ij and D ij to be symmetric with respect to i ↔ j and semi-positive definite (to check this, one can insist − dU/dt to be semi-positive definite). Hence, g ij can also be chosen to be symmetric without loss of generality.

Conversion from Stratonovich to Itô Integral
Now suppose we discretise the time t into t n = t 0 + n∆t, where n = 0, 1, . . . , N.
. . , N}, and the discretised Langevin equation reads: where {ξ n i } is a set of independent Gaussian random variables with zero mean, ξ n i = 0, and Kronecker delta-correlations, ξ m i ξ n j = δ ij δ mn . Of particular interest below are the Itô (a = 0) and Stratonovich (a = 1/2) discretisations. Let usconsider the following two integrals Here, h i is a general function of {x i (t)}. Note that ξ n j √ ∆t on the right-hand side of (42,43) is also called the Wiener process t n +∆t t n η(t) dt. In I S ij , the Stratonovich integral, h i is evaluated at the mid-points {x n+ 1 2 i }, whereas in the Itô integral, I I ij , h i is evaluated at the start-points {x n i } of each time increment. To connect the two integrals, we expand x (42) to give: Finally we can approximate ξ n l ξ n j δ lj (which is valid in the limit ∆t → 0 [19]) to obtain where the conversion term I S→I ij is just the (noiseless) Riemann integral

Dynamical Action
Following a similar derivation for the case d = 1 given above, the path probability for some discretised trajectory {x n i |n = 0, 1, . . . , N} is given by P {x n i } ∝ e −A{x n i } , where the action is [18] where f i , g ij , D ij , and their derivatives are evaluated at {x n+a i }, and D −1 is the inverse matrix of D, with matrix elements D −1 ij . The transition probability from {x 0 i } at time t 0 to {x N i } at time t N can then be written as a path integral in the limit of ∆t → 0. For future reference, we shall also write: where A conv contains all terms which depend on a explicitly. For instance, for additive noise, where g ij , Γ ij , and D ij are constant, the a-explicit term is As already described for d = 1 in Section 2.5, when calculating the EPR, the preferred choice for the time discretisation is a = 1/2 (Stratonovich) so that the pre-exponential product in Equation (48) is the same for any forward and backward pair of paths. With this choice of a = 1/2 (only), A conv is identical for the pair and therefore cancels when the difference of their actions is taken to give the EPR.

Scalar Active Field Theories with Additive Noise
We now turn our focus to field-theoretical models. These require discretisation in space as well as time. We will see that the analysis of time-reversibility for fluctuating hydrodynamics brings additional difficulties with respect to finite dimensional systems.
In what follows, we show that these difficulties can be resolved by carefully choosing the spatial discretisation scheme, as well as the temporal one.
Throughout this section, we address the fluctuating hydrodynamics for a single conserved scalar field, governed by diffusive (Model B-like) dynamics. This describes a system that undergoes phase-separation. We allow for activity but insist that the steadystate EPR must vanish when active terms are switched off. The various considerations set out here generalize in varying degrees to more complex models of the kinds mentioned in the Introduction.
The dynamics of a diffusive conserved scalar order parameter φ(r, t) is governed bẏ where J d is a deterministic current and Λ a spatio-temporal Gaussian white noise current satisfying Here, T is the temperature and Γ is the collective mobility. In principle, Γ = Γ[φ], but we now take it to be constant so that the noise is additive [18,29]. This gives vast technical simplifications that we freely exploit below, with almost no modification to the physics of interest, namely phase separation. For passive systems en route to equilibrium, the deterministic part of the current takes the form This is Model B [4,30]. The chemical potential µ E derives from a free energy F [φ], which is conveniently chosen of the φ 4 -type with a 4 and κ(φ) strictly positive. Phase separation then arises, at mean-field level, whenever a 2 < 0.
Extensions of Model B have recently played a crucial role in understanding phase separation in active systems. In the simplest setting [10,[31][32][33], these theories only retain the evolution of the density field φ, while hydrodynamic [34,35] or polar [36,37] fields can be added if the phenomenology requires. The top-down construction of these field theories, via conservation laws and symmetry arguments, closely retraces the path leading to Model B for passive phase separation [7]. However, locally broken time-reversal symmetry implies that new non-linear terms are allowed. The ensuing minimal theory, Active Model B+ [9,10], includes all terms that break detailed balance up to order O(∇ 4 φ 2 ) in a gradient expansion of the dynamics ofφ [9,10]. It is defined by replacing J in Equation (51) by which contains two activity parameters, λ and ζ, which are independent in more than one dimension. Model B is recovered at vanishing activity (λ = ζ = 0) [4]. Note that we retain constant noise amplitude; such noise need not be thermal in origin in an active system, although it can be in some interesting near-equilibrium cases as will be addressed in Section 4. This model could be further complemented by a coloured noise, a feature that has been recently considered [38,39]. Note also that the decomposition of µ into its equilibrium and nonequilibrium parts is not unique. Since the only defining property of µ A is that it does not derive from a free energy, an arbitrary equilibrium contribution can be moved into it from µ E . For simplicity, we set Γ = 1 without loss of generality and also set ζ = 0. In addition, we will now choose the (positive) square gradient coefficient to be a constant, κ(φ) = κ, following [4,30]. This simplified model was introduced in [31] and is known in the literature as Active Model B (AMB); as just described, it is a special case of AMB+ but sufficient for our present purposes.
In analogy with the finite-dimensional case discussed in Section 2.6.2, the action of AMB can be written as where A conv depends on the scheme employed for the time-discretisation. (Note that the inverse Laplacian in Equation (56) is well defined as a Coulomb integral in either an infinite or periodic domain.) At first sight, it is straightforward to generalise the expression for A conv that was given for finite-dimensional systems in Equation (50) as where s ∈ [0, t], here and below, is a time variable. Importantly, however, no mathematical sense can be given to Equation (57) without an explicit choice of spatial discretisation. Indeed, if we try to retain continuous spatial variables, from Equations (55) and (57), we obtain Here, the presence of δ(0) (the Dirac delta evaluated at zero argument) does not allow a continuum interpretation even in the distributional sense. The problem arises from the fact that Equation (57) contains a functional derivative at point r of a function (∇ · J d ) evaluated at the same spatial location r. As we shall see in Section 3.2, a proper interpretation can be given only after discretising the dynamics in space. We will then find that A conv not only diverges as the continuum limit is taken (resulting in the δ(0) terms), but that it depends on the spatial discretisation scheme used.

Informatic Entropy Production
It is straightforward to notice that A conv is symmetric in time; thus, although it reweights paths in a configuration-dependent manner, it does not contribute to the steadystate IEPR [9], which readṡ where the integral over time is performed within the Stratonovich scheme and the average is taken with respect to noise realizations. For active systems,Ṡ ≥ 0, with equality only if, at the coarse grained scale of the field φ(r, t), the emergent dynamics is reversible. It is perfectly possible, in principle [40], that reversible dynamics do emerge after coarse graining even though the microscopic processes powering the dynamics of φ are very irreversible. However, the generic case in active matter is, of course, to have irreversible dynamics at the mesoscopic scale described by φ(r, t), and hence, have positive IEPR in Equation (59).
Recall that in contrast with the case of a forced thermal particle considered in Equation (35), but just as for the single active particle considered in Equation (34), the informatic entropy production rateṠ given by Equation (59) cannot be interpreted as the ratio between the heat produced and the temperature. There are several reasons for this. Firstly, in a general active setting, even the passive-looking terms in the model (those entering µ E ) need have no connection with interparticle forces: like the active terms, they could emerge from purely behavioural rules among swimming microorganisms, say. Thus, there is no first law, and no direct connection with heat. Second, even in a system where these connections can be made and a first law established, to capture the full heat production of the system, one must consider all microscopic degrees of freedom, not just the coarse-grained fields. However, for systems whose activity can be viewed as a small departure from thermal equilibrium, there is a middle path in which one can embed an active field theory within a larger model whose thermodynamics is consistent at the level of the degrees of freedom actually retained. This approach was developed in [12] and will be reviewed in Section 4.
Meanwhile, as explored in [9][10][11]13], the IEPR has emerged as a useful tool for quantifying the extent to which the behaviour of active complex fluids at hydrodynamic level (as described by φ and/or additional order parameters such as fluid velocity, nematic or polar order, etc.) is irreversible. We give an example of such calculations, which can only be performed numerically and therefore requires further consideration of discretisation, in Section 3.3 below.

Spatial Discretisation
We now discuss spatial discretisation strategies for AMB. The reason is two-fold. First, as we have seen above, we are unable to give a precise mathematical meaning to the action A[φ] of a fluctuating hydrodynamic theory working directly at the continuum level; it is natural to expect, and we confirm this here, that the issue can be solved by discretising the dynamics in space. Second, to numerically integrate any field theory, it is necessary to employ some form of spatial discretisation. A desirable feature of the discretisation used, which becomes crucial if one is interested in measuringṠ, is that the ensuing discrete system respects time-reversal symmetry if the field theory one intends to approximate does. We thus describe here how to perform spatial discretisation of AMB such that detailed balance is always recovered in the equilibrium limit for AMB (λ → 0). For simplicity, we focus on the one-dimensional AMB of finite width L such that x ∈ [0, L] with periodic boundary conditions; extending these results to higher dimensions is straightforward. (Note also that in one dimension, the ζ and λ nonlinearities in Equation (55) are not independent, so we include AMB+ up to the parameter shift λ → λ − ζ/2.) We discretise x into N lattice points with equal lattice spacing ∆x so that N∆x = L, and the density field as φ( Representing the discrete gradient and Laplacian operators as the discretised dynamics reads with η i (t)η j (t ) = δ ij δ(t − t ). Given the spatial reflection symmetry of the underlying model (x → −x), a natural choice is to use midpoint spatial discretisation for the gradient operator, which corresponds to the choice A ij = (δ i+1,j − δ i−1,j )/(2∆x), and hence, B ij = (−δ i+2,j + 2δ ij − δ i−2,j )/(2∆x) 2 .
In the passive limit λ = 0, µ i = (1/∆x)∂F /∂φ i so that Notably, to ensure that the model respects time-reversibility in the passive limit, we are not free in the choice of the discrete gradient and Laplacian operators. Indeed, Equation (62) respects detailed balance only if AA T = A T A = B [20,29], corresponding to ∇ 2 = ∇ · ∇ at the discrete level. Happily, the mid-point spatial discretisation indeed satisfies this condition, and so time is reversible as required.
A separate discretisation issue is to make sense of A conv for AMB, which we found to be divergent if computed directly in the continuum limit. From Equations (50) and (62), we obtain As expected from Equation (58), these terms are divergent as ∆x → 0. Interestingly, A conv not only depends on the choice of the time-discretisation encoded in a ∈ [0, 1] but also on the choice of the spatial discretisation encoded in the matrices A and B. Still, with the Stratonovich choice (a = 1/2), we have that A conv − A R conv = 0. This shows that, even for active fields, A conv does not contribute to the IEPR, which we consider next.

Computing the IEPR
Evaluating the informatic entropy production rateṠ in numerical simulations of fluctuating hydrodynamics exposes a subtlety which is once again related to the precise spatial discretisation used. When simulating the dynamics numerically, it is often preferable to employ Itô's prescription, so that the update at a given timestep depends only on prior data (thus avoiding use of predictor-corrector or other iterative procedures). However, for reasons given in Section 2.5 above,Ṡ is reliably accessible only within the Stratonovich framework. Following standard stochastic calculus rules as recalled in Section 2.6.1 for finite-dimensional systems, one might be tempted to transform the Stratonovich integral definingṠ into an Itô integral that in turn can be computed using trajectories obtained directly from integrating the Itô-discretised time dynamics. Subtleties, however, arise when pursuing this path for stochastic PDEs, which can be fully clarified only by also discretising the spatial dynamics as we do here.
We again consider the case of AMB, for which the IEPR is given by Equation (59). Working directly at the continuum level, let us first try to transform the Stratonovich integral appearing in Equation (59) into an Itô integral by generalising to the infinite dimensional case the conversion term that we have given in Equation (46) for finite dimensions. We obtainṠ where in which ∇ r {1,2} denotes the gradient operator with respect to r {1,2} , and α, β are spatial coordinates. Given that the correction toṠ is given by an integral over space of I S→I (r, r), and that the latter is a divergence, one might speculate that there is no correction due to the Stratonovich to Itô transformation (at least for periodic boundary conditions). However, taking r 1 = r 2 in Equation (65), as required to evaluate Equation (64), produces an undefined δ(0) divergence. We, therefore, consider the entropy production rate of the fully discretised dynamics (61) and perform the same transformation from Stratonovich to Itô integral: where, using Equations (46) and (61), we have In the midpoint spatial discretisation, µ A,i depends only on φ i±1 , while B ij = 0 only when j = i, i ± 2. In this case, we thus obtain from Equation (67) that I S→I = 0. This is, however, not generic and due to the specific form of the non-equilibrium chemical potential µ A of AMB. For example, suppose we had written the IEPR in the following equivalent form, which includes the reversible part of the chemical potential µ E (whose contribution toṠ is a total time derivative that gives zero in the large t limit): Then the same line of reasoning shows that the Stratonovich to Itô conversion factor does not vanish even within the spatial midpoint discretisation scheme. However, with either choice of definition forṠ, the discrete dynamics, as formulated above, is unambiguous and necessarily leads to the same final result; this has indeed been checked numerically for AMB [9].
Let us now revisit the computation that we attempted at the continuum level with Equations (64) and (65). If we employ the following definition of the operator ∇ r δ δφ(r) acting on arbitrary functions g of φ and its derivatives: we obtain drI S→I (r, r) = lim ∆x→0 I S→I , where I S→I obeys Equation (67) and depends on the discretisation scale ∆x. It should be observed, however, that Equation (69) remains only a formal relation because the righthand side can be a divergent quantity. This underlines the fact that to avoid all conceptual ambiguities, we should work with a finite discretisation length. The above analysis shows that, when computing numerically the entropy production rate for field theories, care must be taken with not only the temporal but also the spatial discretisation employed. Using the methodology reviewed here,Ṡ was computed numerically within AMB in [9]. Since this quantity is written as a spatial integral, it is natural in the steady state to associate the first integrand in Equation (64), σ(r) = − lim t→∞ T −1 t −1 t 0 µ Aφ (r, s) ds, with a local IEPR density. When the steady-state system is phase separated, it was shown that for a small T, this density is concentrated at the interfaces between the liquid and vapour phases; see Figure 2, where it scales as T 0 . Away from interfaces in the bulk of each fluid, it instead scales as T 1 . Notably, in active field theories that show deterministic currents in the steady state (such as the uniformly aligned state of a polar active liquid crystal [11]), the IEPR density diverges as T −1 . Observing such scalings numerically can give insight into how and where in the system the active dynamics breaks time reversal symmetry [13].

Thermodynamics of Active Field Theories
In this section, we review what happens when active field theories are minimally coupled to chemical degrees of freedom [12]. The latter can describe the energy flows underlying activity so long as the active motion itself results from locally weak departures from the thermal equilibrium. This allows the recreation of a first law. We will show in these extended theories, analogous ambiguities to those encountered in the previous section arise not just when computing the IEPR but even in defining the stochastic dynamics itself. (This is because multiplicative noise arises in the off-diagonal couplings between the two sectors.) As we shall see, these ambiguities are likewise resolved by careful discretisation.
Our interest is in fluctuating hydrodynamic models of complex fluids in which the activity of a conserved scalar field stems from local consumption of chemical fuel. Prototype examples for such active systems are bacterial suspensions [41,42], acto-myosin networks [43], and self-propelling Janus colloids [44][45][46]. At the continuum level, we therefore address below Active Model B+, as presented in Section 3, which is the leading-order theory of this type. Activity is assumed to be sustained by connecting the active system to reservoirs of fuel and its products; see Figure 3. Our approach relies on systematically constructing the dynamics of the underlying chemical driving field from that of the active field dynamics based on the force-current relations of Linear Irreversible Thermodynamics (LIT), which obey Onsager reciprocal relations [47]. This physically requires that the activity stems from relatively small departures from the local chemical equilibrium. The more microscopic the scale of activity or self-propulsion, the more likely this is to be true: our focus is thus on subcellular systems, or perhaps Janus colloids, rather than collections of animals [48].  Figure 3. Schematic representation of an active system (blue) put in contact with reservoirs of chemical fuel (red) and product (green), which set a constant, homogeneous chemical potential difference ∆µ in the active system. Within our framework, ∆µ embodies the driving parameter which controls the nonequilibrium terms in the dynamics Equations (85) and (87) for the active density field φ and the rate of fuel consumptionṅ. The active system and the chemical reservoirs are surrounded by the thermostat (yellow), which maintains a fixed temperature T. The fluctuations of φ and n lead to the dissipation of heat Q into the thermostat, which quantifies the energetic cost to maintain the whole system away from equilibrium. Note that the physical separation of the reservoirs from the active system, as illustrated, is conceptually helpful but not necessary: in practice, the fuel, active particles and products can all share the same physical domain. Adapted from [12]. Importantly, in some cases, we can construct the extended model (and its discretisation) so that the evolution of the original active fields remains independent of the additional chemical dynamics. This is what we mean by 'embedding' the active field theory into a larger model for which thermodynamic consistency and the first law can reappear; we are not changing the active field theory, just placing it into a more general setting. By accounting for the driving mechanism, we find that the rate of heat production for the active system follows from the full entropy production rate (EPR) measuring the irreversibility of both the active and driving fields, which can however be evaluated from the fluctuations of active fields only. Importantly, the heat rate is distinct from the IEPR,Ṡ [φ], which quantifies the irreversibility of the active field dynamics alone, as previously described.
As stated above, we will find that the coupling of an active field to its driving mechanism generally results in multiplicative noise [29]. It is well known that when dealing with multiplicative (state-dependent) noise, one has to define the specific way in which the noise is evaluated, which affects the time discretisation scheme [18] and generally results in spurious drift terms as we discussed for finite dimensional systems, rather than fields, in Section 2. Moreover, for the reasons already described in Section 3, we also need to pay careful attention to spatial discretisation.

Onsager Coupling in Two-Dimensional System
Before constructing our thermodynamic active field theory, it is instructive to consider a simple example of a two-particle system in which the single particle dynamics seem to be additive, but Onsager reciprocal relations result in multiplicative noise due to crosscoupling in the noise terms.
As a minimal particle-based model for this, let us consider the following dynamicṡ where {Γ x , Γ y } are mobilities, C an arbitrary function of {x, y}, T the temperature, and U the potential. Here, {ν x , ν y } are spurious drift terms that will be defined precisely below. The terms {ξ x , ξ y } are Gaussian white noises with zero mean and correlations given by The dynamics in Equations (71) and (72) can be written in a compact form as where T denotes transpose, and we have introduced the Onsager matrix L given by Such a form for linear coupling between the velocities {ẋ,ẏ} and the forces {−∂ x U, −∂ y U} is inspired by the seminal work of Onsager [47], which demonstrated that L must be positive semi-definite (i.e., det L ≥ 0) for stability.
Due to the fact that the correlations between ξ x and ξ y depend explicitly on {x, y} through C, one has to specify the time discretisation of Equation (71). Changing time discretisation affects the explicit expression of the spurious drift terms {ν x , ν y }, which depend on {Γ x , Γ y , C} and derivatives of C. In practice, we choose the spurious drift terms at a given time discretisation to ensure that the corresponding Fokker-Planck Equation (FPE) for the probability density P(x, y) readṡ Then, the steady-state solution is given by the Boltzmann distribution, P s ∼ e −U/T , as expected for any equilibrium dynamics.
To compute the EPR,Ṡ = lim t→∞ (A R − A)/t, it is convenient to express the dynamic action A (and its time-reversed counterpart A R ) associated with dynamics (71) using the Stratonovich convention (SC). Using Equation (38), the spurious drift terms in SC can then be written as where M is defined by M 2 = L. In practice, decomposing L in terms of the diagonal matrix D (with eigenvalues of L as entries) and of the projector P (constructed from eigenvectors of L), one obtains M = P −1 D 1/2 P. The action follows as The term A conv is a result of the stochastic time integral in the dynamic action and depends on its interpretation (see [18] and Section 2.6.2). It is, however, invariant under time reversal and thus does not contribute to the EPR. Notably, because we use SC, the spurious drift terms {ν x , ν y } do not appear in the first term of A [18]. We deduce the EPR aṡ Note that the product in the integrand is written within SC. Then, we can use the standard chain ruleU =ẋ∂ x U +ẏ∂ y U, leading toṠ = lim t→∞ (U(0) − U(t))/(Tt), which vanishes provided that U does not change in time. Therefore, we have shown that the dynamics (71) are associated with vanishing EPR, as expected at equilibrium.

Spatial Discretisation in Stochastic Field-Theories
The example above makes it clear that our construction of the underlying driving field using LIT is likely to result in multiplicative noise due to cross-coupling noise terms. Therefore, prior to actually constructing our theory, it is useful to discuss the space-discretisation issue that arises in stochastic field theories with multiplicative noise. This issue is very similar to the one encountered in Section 3 for additive noise in the dynamic action of a stochastic field theory, but here, the problem appears already at the Langevin dynamics. To present the discretisation issue in the simplest case, we consider the 1D functional diffusion equation for the density φ of a (thermodynamically) ideal gas with density-dependent diffusivity. In Appendix A, we provide a more general form of the spurious drift terms within LIT.
The 1D functional FPE for the density of an ideal gas is [29,49] which have a steady-state solution P s ∼ exp(−F /T), with F = T dx[φ ln φ − φ] being the ideal-gas free energy. Here, D(x, [φ]) is a functional of the density field φ, which we take to be purely local so that D(x, . This locality will lead below to strong dependence on the discretisation scale along lines seen already in Section 3. The corresponding Itô-Langevin equation is [49] φ( where ξ is a zero mean Gaussian noise with variance ξ(x, t)ξ(x , t ) = δ(x − x )δ(t − t ) and M 2 = Dφ/T, such that the fluctuation dissipation theory is obeyed. The second term in the right hand side of Equation (81) is the spurious drift in the Itô convention, which depends on the noise convention and therefore on the time-discretisation scheme [18,21]. For example, in the Stratonovich convention, this term is changed to We already see that evaluating the spurious drift above in the continuum description is problematic [12,49]. The same issue also arose in Section 3, but only at the level of the dynamic action. Discretising the dynamics in space solved the issue and revealed the actual meaning of ∂ x [δD(x)/δφ(x)]; see Equation (69). Now that we understand the meaning of Equation (81), and specifically the spurious drift term, it is straightforward to show how different choices of spatial discretisation result in different spurious drifts. As a purely mathematical example, consider a system that obeys Equation (80) with D(x, [φ]) =D + ∂ x (∂ x φ) 2 , with some constantD. The nonconstant part D −D can be written as either ∂ x (∂ x φ) 2 or 2(∂ x φ)∂ 2 x φ, which, after discretisation, become, respectively: These of course coincide in the continuum limit, ∆x → 0. A priori, one might expect the spurious drift terms to be independent of this choice of implementation, yet we now show that this is not the case. For D (1) , we obtain where we have used A ij = −A ji . Taking A ij = (δ i+1,j − δ i−1,j )/(2∆x), we deduce [A 2 ] ik A ik = 0, so that Equation (83) is zero. Substituting Equation (83) into Equation (81), we conclude that there is no spurious drift associated with D (1) . (However, this no longer holds when considering higher-order schemes for the gradient matrix A). Choosing instead D (2) , we obtain where we used again A ij = −A ji . Given that A is anti-symmetric, any odd (even) power of A is anti-symmetric (symmetric), so that [A 3 ] ii = 0 and [A 2 ] ii = 0. Then, Equation (84) is always non-zero for any form of the gradient matrix A. This simple example of a 1D ideal gas with density-gradient-dependent diffusivity illustrates that the choice of spatial discretisation can drastically affect the form of the spurious drift terms. Although the chosen form for D is somewhat contrived in this context, we will see that precisely the same discretisation choice will enter our discussion below of spurious drift terms for Active Model B+.

Thermodynamics of a Conserved Active Scalar Field
We now consider the fluctuating hydrodynamics of a conserved active scalar field. Suitable models can be either obtained from explicit coarse-graining of microscopic dynamics [23,49,50] or written from symmetry arguments [8,51]-the prototypical example of the latter route being Active Model B+, Equation (55). The key to embedding such models within a thermodynamic framework is to realize that they omit degrees of freedom (chemical or other), which provide the drive needed to sustain nonequilibrium activity, as described in Figure 3 [12]. Therefore, our approach consists of introducing an additional field, associated in this case with chemical reactions that drive the dynamics away from equilibrium. We then identify the nonequilibrium terms in the original dynamics as a coupling to chemical reservoirs following the framework of LIT [52].
The dynamics of a conserved scalar field φ representing the density of active components for an isotropic material can generally be written as: where F is the free energy, Γ is the mobility, the activity term C is a vector-valued function of φ and its gradients, T is the temperature of the surrounding heat bath, and ν a spurious drift discussed below. The driving force for activity is ∆µ, the chemical potential difference between fuel and products [53][54][55]; see Figure 3. (This is not connected with the chemical potential of the φ field, as defined in Section 3, and here denoted δF /δφ.) An example of such a reaction is the decomposition of hydrogen peroxide involved in the self-propulsion of Janus colloids [44][45][46]. In what follows, n is described as a field fluctuating in space and time, while ∆µ is kept constant and homogeneous. This would be an appropriate approximation for large fuel/product reservoirs and when the chemical fuel and products diffuse much faster than the active particles within the active system [12]. Note that for Active Model B, we have ∆µC = −Γλ∇|∇φ| 2 .
To account for the chemical reactions, we introduce the chemical coordinate n, which is (half) the difference between the local number density of product molecules and that of the fuel molecules. Because the active system is a part of a large nonequilibrium system that relaxes (slowly) towards equilibrium, the explicit dynamics of n can be deduced from LIT [21,52,[56][57][58], in which the thermodynamic fluxes are written as a linear combination of the thermodynamic forces. Identifying J and −∇(δF /δφ) as the current and the thermodynamic force associated with φ, respectively, LIT states that (in the absence of noise) where L is the Onsager matrix. It is clear from Equation (85) that the factor coupling the current J and the force ∆µ is directly given by C (similarly to what we have seen in Section 4.1). Note that, though LIT states linear relations between forces and currents, the coupling factor C need not be linear in φ or its gradients. Accordingly, and because φ is even under time-reversal, Onsager reciprocity relations require that the coupling factor between the currentṅ and the force −∇(δF /δφ) is also C [47]. The dynamics of n follows asṅ where γ is the chemical mobility, which we take constant in what follows. As a result of this assumption, the equation for φ is autonomous and does not rely on knowing the fluctuations of the chemical field n.
In the above, the noises Λ and ξ are Gaussian with zero mean and their correlations are given by The terms Tν in Equation (85) and Tχ in Equation (87) are direct generalizations of the spurious drifts that appears in ordinary stochastic differential equations with multiplicative noise (see Section 2). Their expression is determined by that of C; they depend on both time and space discretisations, as explained in Section 4.2. Both obviously vanish when fluctuations are ignored (T = 0).
The dynamics (85) have been used extensively to reproduce the phase separation of active particles [9,10,31,[59][60][61], with Active Model B+ as a leading example of such theories. In these works, the dynamics of the driving chemicals were not considered so that the noise Λ seems to be purely additive. For this reason, and because previous studies were not concerned with thermodynamic consistency, the term Tν was missing. Where possible, the simplest approach to embedding Equation (85) unchanged within a larger, thermodynamically consistent model is therefore to seek a discretisation scheme (that is, an interpretation of the original stochastic field theory) in which this spurious drift becomes zero.
To date, we did not specify the explicit form of ν and χ. As explained above, to do so requires the discretised version of the dynamics, (85) and (87), where we focus on 1D for simplicity:φ Here, ψ i = (∂F /∂φ i )/∆x, and the coupling term C i = C(φ i , ∑ j A ij φ j , . . . ) depends on φ and its gradients. The discrete noise terms {Λ i , ξ i } are Gaussian with zero mean and correlations given by Given that the correlations between Λ i and ξ i depend on the variable φ i through the coupling term C i , one has to specify the temporal discretisation scheme of Equation (91).
In what follows, we choose the Stratonovich convention, which allows one to use the standard rules of differential calculus [29]. As found in Section 2.5 above, there are compelling reasons to prefer this choice when deriving the expression of the heat rate or EPR.
The associated FPE for the probability density P({φ i , n i }, t) can then be derived following standard methods [29] aṡ where we have introduced the matrix M i defined by M i M T i = L i . In the continuum limit of small ∆x, it follows using Equation (69) that Equation (93) converges to the standard functional FPE for the probability density P([φ(x), n(x)], t) [49,62]. Importantly, by taking allows us to simplify Equation (94) as The matrix M i can be written as Substituting the expression of M i in Equation (96), we find that ν i vanishes for any C i in d = 1 (it can still potentially be non-zero in higher dimensions), while the expression of χ i is For the specific coupling term x φ corresponding to Active Model B [63] (and in d = 1 AMB+ also), it is possible to write C using at least two different discretisation schemes, for example, those used in Equation (82). The results in Equations (83) and (84) are then also appropriate in our case and illustrate that the choice of spatial discretisation drastically affects the form of the spurious drift terms appearing in the Langevin equations at the field level. Specifically, for the choice in Equation (83), the spurious drift in Equation (98) vanishes, while for the choice Equation (84), it does not vanish but instead diverges as 1/∆x. Clearly, therefore, any attempt to numerically code the coupled Langevin Equations (85) and (87) that either ignores the spurious drift terms or claims to calculate them without reference to the discretisation scheme used risks very large errors in the simulated dynamics.

Calculation of the Heat Production Rate
We next calculate the heat production rate [12,14] where the average is taken with respect to noise realizations (or P {J,ṅ} t 0 ). Note thaṫ Q/T is the full EPR of our enlarged, thermodynamic model. The conserved field φ and its drivingṅ dynamics can be written aṡ where the noise and spurious drift terms obey Equations (92) and (94), respectively. Generalizing beyond the dynamics in Equations (85) and (87), we now consider an arbitrary Onsager matrix L, with the only constraint that it should be positive semi-definite (det L ≥ 0). Following [18,19] and similarly to the finite-dimensional case considered in Section 2.6.2, the path probability P ∼ e −A associated with Equation (100) is defined by where, as a consequence of the Stratonovich discretisation, no spurious drift terms appear in the expression (101) [18]. Note that, as before (see Sections 3 and 4.1), A conv is even under time-reversal and is not written explicitly in Equation (101) where lim t→∞ 1 t t 0 · ≡ · t is the steady-state time average. In steady state, even in spatially inhomogeneous systems, such as phase separation, the two averages are the same and the temporal one may be omitted. Note that the product above is interpreted here and in what follows with the Stratonovich convention.
Integrating by parts the second term in Equation (103) and usingφ = −∇ · J, we obtain V J · ∇(δF /δφ) dr = d F /dt, which vanishes in steady state, yieldinġ As a result, the steady-state heat rateQ equals the rate of work injected by the nonequilibrium drive ∆µ to sustain the dynamics away from equilibrium. This is equivalent to the first law of thermodynamics, as expected, when the path probabilities include all thermodynamically relevant fields. For equilibrium dynamics where ∆µ derives from the chemical free energy F ch , (∆µ = −δF ch /δn), the heat rate rate vanishes in steady state (Q = −d F ch /dt = 0), as expected. Activity is instead introduced by the fact that ∆µ is held away from equilibrium. (Note that the expression (104) would actually be the same if insteadṅ was held constant and ∆µ allowed to fluctuate.) Substituting the chemical dynamics (87) into Equation (104), we deducė Hence, the heat rate can be separated into (i) a homogeneous contribution γV∆µ 2 corresponding to a background term independent of the fluctuations of the active and chemical fields {φ, n} and (ii) a contribution determined only by the fluctuations of the active field φ, with no contribution from the fluctuations of the chemical coordinate n. The presence of n is nonetheless crucial in determining the form of the heat production rate. Interestingly, the homogeneous contribution is eliminated when considering differences in the heat rates at constant ∆µ, for example, comparing a state of uniform φ with a phaseseparated one and/or finding the effect on heat rate of changing parameters in the free energy F . We continue by comparing the heat rate from Equation (105) with the IEPR as introduced in Section 3 and used in previous works [9,54,64]. Substituting into Equation (105) the expression of ∇(δF /δφ) taken from the dynamics (85) yieldṡ where the IEPRṠ of the φ dynamics reads [9,54,64] S = ∆µ (Note that for AMB, this equates by partial integration to Equations (59) and/or (68) given above.) Clearly, the second line in Equation (106) depends on the spurious drift terms, but it also depends directly on the evaluation of the stochastic integral V C · Λ dr and thereby on the discretisation scheme used to evaluate the heat rate (including spatial discretisation). We show below that for AMB(+) in d = 1, a discretisation scheme can be found for which {ν, χ} = {0, 0} and V C · Λ dr = 0. In this and other cases for which all these terms vanish, we arrive at a simple relation involving the heat rateQ and the IEPR,Ṡ: From the semi-positivity of the Onsager matrix L, which ensures det L = λγ − C 2 ≥ 0, it then follows that TṠ is a lower bound toQ. The bound is saturated when J andṅ are proportional (det L = 0): In such a case, the fluctuations ofṅ are determined by that of J, so the irreversibility of the whole dynamics can be found from the trajectories of J alone. As noted in Section 3,Ṡ can be written as the spatial integral of a local quantity σ(r), and we see that so can be the chemical contribution in Equation (108). Thus,Q = q(r) dr withq a local heat production rate density. As we found with the IEPR, it is interesting to examine where, in phase-separated system, this density is large or small (see Figure 4).
A specific choice of discretisation for which {ν, χ} = {0, 0} and V C · Λ dr = 0, such that Equation (108) holds for AMB+ in d = 1, is that of Equation (83). To establish this, we evaluate ∑ i C i Λ i transforming it into an Itô product (see Section 2.6.1) where we have used again that C i is independent of n i . From Equations (98), (106), and (109), it follows that the relation between the heat rate and the IEPR depends on ∑ j A ij (∂C i /∂φ j ), which vanishes for the discretisation of Equation (83). In this case, a direct comparison of the heat-rate with previous results [9,10,31,[59][60][61], and specifically with the results of Section 3, which did not have spurious drift terms is valuable. In Figure 4, we provide such a comparison. For a phase-separated profile, as shown in Figure 4a,b, the leading order ofQ − γV∆µ 2 scales like T 0 , and it reaches a finite value at T = 0. Hence, the heat rateQ is not only determined by the background term γV∆µ 2 at zero temperature; it now also depends on the mean-field density profile. In contrast, TṠ scales like T and thus vanishes at T = 0, see Figure 4c, as already reported in [9] and in Figure 2. Notably, while the IEPR is maximal on the interface between phases, showing maximal irreversibility of the fluctuating φ dynamics, the heat rate density is suppressed there. This suggests that the chemical reactions are, in the interfacial zone, producing less heat because they are instead doing local work against F to sustain the nonequilibrium coexistence. Thus, bothṠ andQ can differently reveal useful insights into the dynamics of the system.
These results also fully confirm that the IEPR, which considers the irreversibility of the φ dynamics alone, does not capture the full energetic cost of creating phase separation away from equilibrium, as the heat-rateQ does. In fact, if TṠ was indeed a measure of the full energetic cost, a nonequilibrium active phase separation could be sustained at zero energy cost as T → 0, contradicting the basic thermodynamic notion that activity is powered by constant input energy that is ultimately dissipated as heat.

Concluding Remarks
In this paper, we have addressed several conceptual issues arising from the stochastic PDEs (SPDEs) used by physicists to describe the fluctuating hydrodynamics of complex fluids. These conceptual issues arise because the continuum limit, while implicit in the notation used to write down these SPDEs, is generally either nonexistent or at least problematic [6]. The usual physicist's defence is to protest that there is always a short-scale cutoff (set by molecular physics), so the SPDEs are really only a short-hand for a discretised version of the same equations. Rarely are such versions closely examined, and often, they are not even specified unless numerical work is actually undertaken (sometimes, not even then). We hope to have convinced the reader that a more careful study of the meaning of these equations based on careful and consistent discretisation strategies is warranted.
In particular, attention must be paid to achieving detailed balance at the discrete level in the case of equilibrium systems. This is not a new remark (see, e.g., [65]) but is brought into sharper focus by the desire to numerically evaluate the entropy production rate (EPR). This desire, driven by recent work on active rather than equilibrium complex fluids, requires careful study of the discretisation scheme used to establish the path weights (or dynamical action), from which, via the laws of stochastic thermodynamics, the EPR can be calculated. Although the scheme to embed SPDEs within thermodynamically consistent description is not unique a priori, our framework provides a minimal approach to do so without LIT. Interestingly, LIT is also the starting point for a large class of active field theories, known as active gels, which have been extremely successful in capturing the dynamics of complex biological systems, such as acto-myosin networks and living tissues [55,56,58,66].
For SPDEs with purely additive noise (such as Model B and its active variants), these problems are first encountered in computing the informatic EPR (IEPR) for fluctuating active fields, which quantifies the irreversibility of the coarse-grained order parameter dynamics without concern for the underlying heat flows. However, the same problems are accentuated further when one addresses these heat flows by minimally coupling the active order parameter fields to an underlying chemical process governed by linear irreversible thermodynamics. In this case, the active terms in the stochastic hydrodynamic equations for the order parameters become off-diagonal Onsager couplings in the enlarged model. The result is that the off-diagonal noise is multiplicative, even when the original noise in the order parameter sector was not. This necessitates the treatment of spurious drift terms directly in the Langevin dynamics; like similar terms in the dynamical action, these are dependent on both temporal and spatial discretisation schemes. Moreover, unless they can be eliminated altogether by careful design of such schemes, these terms diverge in the spatial continuum limit, ∆x → 0. In this setting, and presumably also in other models of fluctuating complex fluids that involve multiplicative noise (for example, Model B with a composition-dependent mobility), relatively minor oversights in numerical implementation could therefore lead to errors in the generation of Langevin trajectories that are not merely O(1), but unbounded, as the continuum limit is approached.
Note that M can be defined as symmetric, in which case, it is the square root of M s [18]. Note also that the choice of different time-discretisation schemes only affects the dissipative part of the generalized mobility matrix (its symmetric part M s ), while its reactive part M a (the antisymmetric part) [21,58] contributes a term that is unaffected by time-discretisation.
Importantly, and as explained in detail in Sections 3.3 and 4.2, there are problems with the continuous description of the spurious drift, specifically in cases where M involves spatial gradients. Therefore, to make sense of the expressions in this Appendix, they must be discretised. Following principles already laid down in the main text, we discretise space and write ψ a (r) → ψ a;(i,j,k) with {a, b, c} letters denoting the various fields and {i, j, k} referring to the spatial discretisation r → (i∆x, j∆y, k∆z). Then, the spatially discretised spurious drift in the Stratonovich convention is written as where we suppress the time dependence.