Time-reversal symmetry violations and entropy production in field theories of polar active matter

For the computations we present here, we will assume that the steady-state dynamics relaxes onto a 'ground-state' trajectory ${\left({\rho }_{0}\left(\boldsymbol{x},t\right),{\boldsymbol{P}}_{0}\left(\boldsymbol{x},t\right)\right)}_{t\in \left(-\infty ,\infty \right)}$ ³ in the limit D → 0. That is, we assume that the probability distribution over trajectories concentrates on a single path as D → 0, and express this by

$$\begin{equation}\left(\rho ,\boldsymbol{P}\right)D\to 0{\to }\enspace \left({\rho }_{0},{\boldsymbol{P}}_{0}\right),\end{equation} \tag{ 20 }$$

where the limit is understood in the almost sure sense. We also assume that (ρ₀, P ₀) solves (1)–(3) at D = 0 and that this limit is unique up to possible degeneracies arising from rotational invariance. Firstly our aim will be to classify the ground-states that satisfy the pathwise equilibrium condition

$$\begin{equation}\dot {\mathcal{S}}\left[{\rho }_{0},{\boldsymbol{P}}_{0}\right]=0.\end{equation} \tag{ 21 }$$

In particular, if both (20) and (21) hold, the dynamics must have χ > −1. Clearly, the pathwise equilibrium ground-states include those that are invariant under $\mathcal{T}$ in (14), meaning that they satisfy

$$\begin{equation}\left({\rho }_{0}\left(\boldsymbol{x},t\right),{\boldsymbol{P}}_{0}\left(\boldsymbol{x},t\right)\right)=\left({\rho }_{0}\left(\boldsymbol{x},-t\right),-{\boldsymbol{P}}_{0}\left(\boldsymbol{x},-t\right)\right),\end{equation} \tag{ 22 }$$

which follows from the fact that $\dot {\mathcal{S}}{\circ}\mathcal{T}=-\dot {\mathcal{S}}$. The constant homogeneous isotropic state with ρ₀ = const. and P ₀ = 0 provides an example of such a state. On the other hand, the PL state with ρ₀ > 1 and $\vert {\boldsymbol{P}}_{0}\vert =\sqrt{{\rho }_{0}-1}$ clearly violates $\mathcal{T}$ alone. This is where rotational invariance arises as an important symmetry principle, because it implies that $\mathcal{P}$ (and thus $\dot {\mathcal{S}}$) must be invariant under the parity transformation

$$\begin{equation}P=\begin{cases}\rho \left(\boldsymbol{x},t\right){\mapsto}\rho \left(-\boldsymbol{x},t\right),\quad \hfill \\ \boldsymbol{P}\left(\boldsymbol{x},t\right){\mapsto}-\boldsymbol{P}\left(-\boldsymbol{x},t\right),\quad \hfill \end{cases}\end{equation} \tag{ 23 }$$

which translates to the statement that a flock is equally likely to travel to the left as to the right. Now, if the ground-state trajectory (ρ₀, P ₀) is P $\mathcal{T}$-symmetric, i.e. it satisfies

$$\begin{equation}\left({\rho }_{0}\left(\boldsymbol{x},t\right),{\boldsymbol{P}}_{0}\left(\boldsymbol{x},t\right)\right)=\left({\rho }_{0}\left(-\boldsymbol{x},-t\right),{\boldsymbol{P}}_{0}\left(-\boldsymbol{x},-t\right)\right),\end{equation} \tag{ 24 }$$

then it follows that it is pathwise equilibrium. Perhaps surprisingly then, one realises that the constant homogeneous PL state in fact is pathwise equilibrium since it satisfies P $\mathcal{T}$. However, this is to be expected: a charged particle gyrating at constant frequency in the plane perpendicular to an imposed magnetic field is certainly in equilibrium (although here $\mathcal{T}$ should be replaced by C $\mathcal{T}$ to include charge conjugation). Interestingly, these observations also imply that if rotational symmetry is broken a priori, for example by driving the system with an external electric field, then $\mathcal{P}$ would no longer be P-invariant and the PL state would have χ = −1. We also note that the fact that χ > −1 in both the homogeneous isotropic and PL states also follows directly from (18) by evaluating the integral at constant (ρ₀, P ₀).

In order to go beyond the D = 0 dynamics, we must take account of fluctuations. We do so by assuming that the fluctuating dynamics admit an expansion in small $\sqrt{D}$, following [33, 35], so that

$$\begin{equation}\rho ={\rho }_{0}+{\rho }_{1}{D}^{1/2}+O\left(D\right),\end{equation} \tag{ 25 }$$

$$\begin{equation}\boldsymbol{P}={\boldsymbol{P}}_{0}+{\boldsymbol{P}}_{1}{D}^{1/2}+O\left(D\right).\end{equation} \tag{ 26 }$$

Furthermore, we restrict here to the case where (ρ₀, P ₀) is constant and homogeneous. By substituting (25) and (26) into the equations of motion in (1)–(3) and collecting terms, we obtain at order D^1/2

$$\begin{equation}{\partial }_{t}{\rho }_{1}=-w\nabla \cdot {\boldsymbol{P}}_{1},\end{equation} \tag{ 27 }$$

$$\begin{equation}{\partial }_{t}{\boldsymbol{P}}_{1}+\lambda {\boldsymbol{P}}_{0}\cdot \nabla {\boldsymbol{P}}_{1}=-\frac{\delta {F}_{\mathrm{L}}}{\delta \boldsymbol{P}}\left[{\rho }_{1},{\boldsymbol{P}}_{1}\right]+{\boldsymbol{\eta }}_{1}.\end{equation} \tag{ 28 }$$

Here, F_L is the quadratic functional

$$\begin{equation}{F}_{\mathrm{L}}\left[\rho ,\boldsymbol{P}\right]={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left({f}_{\mathrm{L}}\left(\rho ,\boldsymbol{P}\right)+\frac{1}{2}{\left({\nabla }_{\alpha }{P}_{\beta }\right)}^{2}+\boldsymbol{P}\cdot \nabla {{\Phi}}_{\mathrm{L}}\left(\rho ,\boldsymbol{P}\right)\right),\end{equation} \tag{ 29 }$$

where we have defined the local free energy f_L by

$$\begin{equation}{f}_{\mathrm{L}}\left(\rho ,\boldsymbol{P}\right)=\frac{{a}_{0}}{2}\vert \boldsymbol{P}{\vert }^{2}-\rho {\boldsymbol{P}}_{0}\cdot \boldsymbol{P}+{\left({\boldsymbol{P}}_{0}\cdot \boldsymbol{P}\right)}^{2},\end{equation} \tag{ 30 }$$

in addition to the linearised effective pressure

$$\begin{equation}{{\Phi}}_{\mathrm{L}}={w}_{1}\rho -\kappa {\boldsymbol{P}}_{0}\cdot \boldsymbol{P}.\end{equation} \tag{ 31 }$$

Moreover, a₀ = 1 − ρ₀ + | P ₀|² and P ₀ satisfies a₀ P ₀ = 0, while η ₁ is a mean-zero Gaussian white noise with

$$\begin{equation}\left\langle {\eta }_{1\alpha }\left(\boldsymbol{x},t\right){\eta }_{1\beta }\left({\boldsymbol{x}}^{\prime },{t}^{\prime }\right)\right\rangle =2{\delta }_{\alpha \beta }\delta \left(\boldsymbol{x}-{\boldsymbol{x}}^{\prime }\right)\delta \left(t-{t}^{\prime }\right).\end{equation} \tag{ 32 }$$

We may perform a similar procedure in order to obtain an expansion of $\dot {\mathcal{S}}$ from (18) in small D of the form

$$\begin{equation}\dot {\mathcal{S}}\left(D\right)={\dot {\mathcal{S}}}_{-1}{D}^{-1}+{\dot {\mathcal{S}}}_{0}+{\dot {\mathcal{S}}}_{1}D+O\left({D}^{2}\right).\end{equation} \tag{ 33 }$$

In appendix B we show these are the only possible terms that could enter in the expansion of $\dot {\mathcal{S}}$, i.e. that there are no terms of half-integer order in D. Also, from the asymptotic scaling relation (19) and the positivity of the EPR, we know that ${\dot {\mathcal{S}}}_{k}=0$ for k < χ and that the leading order coefficient ${\dot {\mathcal{S}}}_{\chi }{ >}0$. By explicitly computing this expansion of $\dot {\mathcal{S}}$ to order D⁰, it is straightforward to show that

$$\begin{equation}\chi \quad \begin{cases}1,\quad \text{isotropic},\quad \hfill \\ =0,\quad \text{polar}\;\text{liquid},\quad \hfill \end{cases}\end{equation} \tag{ 34 }$$

and so χ > −1 in the homogeneous phases as argued for above. For the explicit calculation of (34), we refer to appendix C for the details. From simulations we further find that in fact χ = 1 in the isotropic phase, as shown in figure 2, although we do not explicitly compute ${\dot {\mathcal{S}}}_{1}$. In the PL case, we obtain an explicit expression for ${\dot {\mathcal{S}}}_{0}$ given by

$$\begin{align}\hfill {\dot {\mathcal{S}}}_{0}& ={P}_{0}^{2}\left(2{w}_{1}-\kappa +\lambda \right){\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {\rho }_{1}{\partial }_{{\Vert}}{P}_{{\Vert}}\right\rangle +{P}_{0}\left(w-2{P}_{0}^{2}\kappa \right){\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {P}_{{\Vert}}{\partial }_{\perp }{P}_{\perp }\right\rangle \hfill \\ \hfill & \quad +2{P}_{0}\kappa {\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \left({\partial }_{\perp }{P}_{\perp }\right)\left({\nabla }^{2}{P}_{{\Vert}}\right)\right\rangle .\hfill \end{align} \tag{ 35 }$$

In (35) we use subscripts ∥ and ⊥ to denote components of P ₁ and ∇ that are parallel and perpendicular to P ₀ respectively.

Zoom In Zoom Out Reset image size

Consistently with (34), equation (35) implies that ${\dot {\mathcal{S}}}_{0}\sim {P}_{0}=\vert {\boldsymbol{P}}_{0}\vert$ for P₀ ≪ 1. In fact, this could have been predicted without explicitly performing the systematic expansion of $\dot {\mathcal{S}}$ in small D. Indeed, if we were to imagine expanding the integral expression for $\dot {\mathcal{S}}$ in (18) using the series representations (25) and (26) at P ₀ = 0, we see from simple power counting that the only combinations of fields that could possibly appear within the integrand at order D⁰ are of the form

$$\begin{equation}\langle {\rho }_{2}\rangle ,\hspace{25.0pt}\langle \nabla \cdot {\boldsymbol{P}}_{2}\rangle ,\hspace{25.0pt}\langle {\rho }_{1}^{2}\rangle ,\hspace{25.0pt}\langle \vert {\boldsymbol{P}}_{1}{\vert }^{2}\rangle ,\hspace{25.0pt}\langle {\rho }_{1}\nabla \cdot {\boldsymbol{P}}_{1}\rangle ,\dots \end{equation} \tag{ 36 }$$

Now, the symmetry $\dot {\mathcal{S}}{\circ}\mathcal{T}=-\dot {\mathcal{S}}$ excludes all of ρ₂, ${\rho }_{1}^{2}$ and | P ₁|² from entering, while ∇ ⋅ P ₂ would just integrate to zero over $\mathcal{V}$. For the final average in (36), observe that (27) implies

$$\begin{equation}\langle {\rho }_{1}\nabla \cdot {\boldsymbol{P}}_{1}\rangle =-{w}^{-1}\langle {\rho }_{1}{\partial }_{t}{\rho }_{1}\rangle =\frac{-1}{2w}{\partial }_{t}\langle {\rho }_{1}^{2}\rangle =0.\end{equation} \tag{ 37 }$$

Hence, there are in fact no nontrivial contributions that could enter in the expansion of $\dot {\mathcal{S}}$ at order D⁰ when P ₀ = 0, so we must have χ ⩾ 1 in the isotropic phase.

Since $\dot {\mathcal{S}}$ is O(D) in the isotropic phase, we in fact recover effective equilibrium in the limit D → 0. To see this, we transform the linearised equations of motion (27) and (28) to Fourier space. Throughout we use the convention that the Fourier coefficients $\hat{h}\left(\boldsymbol{q}\right)$ of a function h( x ) are given by

$$\begin{equation}\hat{h}\left(\boldsymbol{q}\right)={\mathcal{V}}^{-1}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace h\left(\boldsymbol{x}\right)\mathrm{exp}\left(-\mathrm{i}\boldsymbol{x}\cdot \boldsymbol{q}\right),\end{equation} \tag{ 38 }$$

where we with slight abuse of notation denote by $\mathcal{V}=\mathrm{v}\mathrm{o}\mathrm{l}\left(\mathcal{V}\right)$. We thus obtain the set of equations

$$\begin{equation}{\partial }_{t}{\hat{\rho }}_{1}=-iw\boldsymbol{q}\cdot {\hat{\boldsymbol{P}}}_{1},\end{equation} \tag{ 39 }$$

$$\begin{equation}{\partial }_{t}{\hat{\boldsymbol{P}}}_{1}=-{\Gamma}\left(q\right){\hat{\boldsymbol{P}}}_{1}-i{w}_{1}\boldsymbol{q}{\hat{\rho }}_{1}+{\hat{\boldsymbol{\eta }}}_{1},\end{equation} \tag{ 40 }$$

where we have defined the damping coefficient Γ(q) = a₀ + q² and the noise term ${\hat{\boldsymbol{\eta }}}_{1}\left(\boldsymbol{q},t\right)$ is mean zero, Gaussian and white with autocovariance

$$\begin{equation}\langle {\hat{\eta }}_{1\alpha }\left(\boldsymbol{q},t\right){\hat{\eta }}_{1\beta }^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =2{\mathcal{V}}^{-1}{\delta }_{\alpha \beta }{\delta }_{\boldsymbol{q},{\boldsymbol{q}}^{\prime }}\delta \left(t-{t}^{\prime }\right).\end{equation} \tag{ 41 }$$

Using the mapping

$$\begin{equation}{\hat{\rho }}_{1}\left(\boldsymbol{q},t\right)=x\left(\boldsymbol{q},t\right)+iy\left(\boldsymbol{q},t\right),\end{equation} \tag{ 42 }$$

$$\begin{equation}{\partial }_{t}{\hat{\rho }}_{1}\left(\boldsymbol{q},t\right)={v}_{x}\left(\boldsymbol{q},t\right)+i{v}_{y}\left(\boldsymbol{q},t\right),\end{equation} \tag{ 43 }$$

and setting X = (x, y), V = (v_x , v_y ) we immediately see that these follow standard equilibrium Langevin equations for an underdamped particle in a harmonic potential [2],

$$\begin{equation}\dot {\boldsymbol{X}}=\boldsymbol{V},\end{equation} \tag{ 44 }$$

$$\begin{equation}\dot {\boldsymbol{V}}=-{\Gamma}\boldsymbol{V}-{\nabla }_{X}U+\sqrt{2{\Gamma}T}\boldsymbol{\zeta },\end{equation} \tag{ 45 }$$

where the potential U = w w₁ q²| X |²/2. The final degrees of freedom in (39) and (40) are captured by the transverse component V _T = (v_Tx , v_Ty ) of ${\hat{\boldsymbol{P}}}_{1}$ with respect to q , i.e.

$$\begin{equation}{v}_{\mathrm{T}x}\left(\boldsymbol{q},t\right)+i{v}_{\mathrm{T}y}\left(\boldsymbol{q},t\right)=-iw{\boldsymbol{q}}_{\perp }\cdot {\hat{\boldsymbol{P}}}_{1}\left(\boldsymbol{q},t\right),\end{equation} \tag{ 46 }$$

and q _⊥ is perpendicular to q with | q _⊥| = q. Again this follows an equilibrium Langevin equation,

$$\begin{equation}{\dot {\boldsymbol{V}}}_{\mathrm{T}}=-{\Gamma}{\boldsymbol{V}}_{\mathrm{T}}+\sqrt{2{\Gamma}T}{\boldsymbol{\zeta }}_{\mathrm{T}}.\end{equation} \tag{ 47 }$$

In (45) and (47) the noise terms ζ , ζ _T are mean zero unit variance Gaussian white noises, and interestingly the effective temperature T is defined by

$$\begin{equation}T\left(q\right)=\frac{1}{2\mathcal{V}}\frac{{w}^{2}{q}^{2}}{{a}_{0}+{q}^{2}}.\end{equation} \tag{ 48 }$$

Since the modes V , V _T are independent for all q , the dependence of the effective temperature T on q does not lead to any current in phase space. At higher order in D, however, modes ${\hat{\rho }}_{1}\left(\boldsymbol{q},t\right)$ and ${\hat{\boldsymbol{P}}}_{1}\left(\boldsymbol{q},t\right)$ are coupled at different wavevectors q via the nonlinear terms in the equation for the polar density in (3). In particular, these terms couple heat baths at different temperatures T(q), driving the dynamics at the next order away from equilibrium.

Linear theory also allows us to make quantitative predictions about ${\dot {\mathcal{S}}}_{0}$ from (35) in the PL phase. Indeed, transforming this to Fourier space we obtain

$$\begin{equation}{\dot {\mathcal{S}}}_{0}\left({\Lambda}\right)/\mathcal{V}=\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\left\langle {\left({\hat{\boldsymbol{u}}}^{\text{pl}}\right)}^{{\dagger}}{\dot {\sigma }}^{\text{pl}}{\hat{\boldsymbol{u}}}^{\text{pl}}\right\rangle =\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\mathrm{Tr}\left({\dot {\sigma }}^{\text{pl}}{\mathcal{C}}^{\text{pl}}\right),\end{equation} \tag{ 49 }$$

where the Hermitian matrix ${\dot {\sigma }}^{\text{pl}}$ is given by

$$\begin{equation}{\dot {\sigma }}^{\text{pl}}=\frac{i{P}_{0}}{2}\left(\begin{matrix}\hfill 0\hfill & \hfill {P}_{0}\left(2{w}_{1}-\kappa +\lambda \right){q}_{{\Vert}}\hfill & \hfill 0\hfill \\ \hfill {P}_{0}\left(\kappa -2{w}_{1}-\lambda \right){q}_{{\Vert}}\hfill & \hfill 0\hfill & \hfill \left(w-2\left({P}_{0}^{2}+{q}^{2}\right)\kappa \right){q}_{\perp }\hfill \\ \hfill 0\hfill & \hfill \left(2\left({P}_{0}^{2}+{q}^{2}\right)\kappa -w\right){q}_{\perp }\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ 50 }$$

In addition, in (49) we have defined the vector

$$\begin{equation}{\hat{\boldsymbol{u}}}^{\text{pl}}={\left({\hat{\rho }}_{1},{\hat{P}}_{{\Vert}},{\hat{P}}_{\perp }\right)}^{\mathrm{T}}\end{equation} \tag{ 51 }$$

of Fourier modes, as well as the matrix ${\mathcal{C}}^{\text{pl}}\equiv \left({\mathcal{C}}_{ij}^{\text{pl}}\right)$ of equal-time correlators by

$$\begin{equation}{\mathcal{C}}_{ij}^{\text{pl}}\left(q\right){\delta }_{\boldsymbol{q},{\boldsymbol{q}}^{\prime }}=\left\langle {\hat{u}}_{i}^{\text{pl}}\left(\boldsymbol{q},t\right){\left({\hat{u}}_{j}^{\text{pl}}\left({\boldsymbol{q}}^{\prime },t\right)\right)}^{{\ast}}\right\rangle .\end{equation} \tag{ 52 }$$

The sum in (49) runs over modes with wavenumbers smaller than the ultraviolet cutoff Λ, which is introduced since the sum is divergent with Λ → ∞. This is sometimes seen in field theories of this kind, since they are often derived based on the assumption that they are only valid down to a certain length scale. In appendix C we show from this that ${\dot {\mathcal{S}}}_{0}\left({\Lambda}\right)$ as predicted by the linear theory diverges in the ultraviolet as Λ². Although the closed form expressions for the correlators entering in (49) are too algebraically involved to report explicitly, they may be calculated straightforwardly by numerical methods. By doing this, we may calculate the corresponding sum in (49) and quantitatively compare the results with measurements of $\dot {\mathcal{S}}$ from simulations. In figure 2 we plot the results obtained from simulations, which show good agreement with the analytical results. Plot a) in figure 2 demonstrates that the EPR is O(D) in the isotropic phase, as well as the predicted O(D⁰) scaling in the PL phase. In particular, both remain finite at fixed Λ as D → 0, with $\dot {\mathcal{S}}\to 0$ in the isotropic phase.

Previous studies have investigated the nonlinear solutions to (1)–(3) at D = 0, and particularly interesting to our present context are the seminal contributions by Solon et al [43, 49] on the structure of the banded profiles. These are effectively one-dimensional travelling wave solutions that are invariant along the direction perpendicular to the motion. We thus write P ₀ = (P₀, 0) without loss of generality, and look for solutions of the form

$$\begin{equation}{\rho }_{0}\left(\boldsymbol{x},t\right)\equiv {\tilde {\rho }}_{0}\left(x-ct\right),\end{equation} \tag{ 53 }$$

$$\begin{equation}{P}_{0}\left(\boldsymbol{x},t\right)\equiv {\tilde {P}}_{0}\left(x-ct\right).\end{equation} \tag{ 54 }$$

Direct substitution then allows us to deduce a set of equations for ${\tilde {\rho }}_{0}$ and ${\tilde {P}}_{0}$ in terms of the variable z = x − ct given explicitly by

$$\begin{equation}{\tilde {\rho }}_{0}={\rho }_{g}+\frac{w}{c}{\tilde {P}}_{0},\end{equation} \tag{ 55 }$$

$$\begin{equation}{\tilde {P}}_{0}^{{\prime\prime}}=-\left(c-\frac{w{w}_{1}}{c}-\lambda {\tilde {P}}_{0}\right){\tilde {P}}_{0}^{\prime }-\left({\rho }_{g}-1+\frac{w}{c}{\tilde {P}}_{0}-{\tilde {P}}_{0}^{2}\right){\tilde {P}}_{0},\end{equation} \tag{ 56 }$$

where primes denote differentiation with respect to z. Equation (56) can be mapped onto a Newton problem for a particle in a potential under the influence of a nonlinear drag, and all stable orbits in the $\left({\tilde {P}}_{0},{\tilde {P}}_{0}^{\prime }\right)$ plane with ${\tilde {P}}_{0}{\geqslant}0$ can be uniquely identified with a pair (c, ρ_g ). In terms of the stochastic equations in (1)–(3), it is assumed that the noise selects the stable steady-state profile (of which there are infinitely many [43, 49]).

Importantly, these solutions to (54) break both $\mathcal{T}$ and P $\mathcal{T}$-symmetry. Thus we expect the microphase-separated steady-state to have χ = −1. Indeed, using the travelling wave ansatz in (53) and (54) we deduce two expressions for ${\dot {\mathcal{S}}}_{-1}={\mathrm{lim}}_{D\to 0}\enspace D\dot {\mathcal{S}}\left(D\right)$, that are

$$\begin{equation}{\dot {\mathcal{S}}}_{-1}/\mathcal{V}=\frac{1}{2L}{\int }_{-L}^{L}\mathrm{d}z{\left(\left(1-{\tilde {\rho }}_{0}\right){\tilde {P}}_{0}+{\tilde {P}}_{0}^{3}-{\tilde {P}}_{0}^{{\prime\prime}}\right)}^{2}\end{equation} \tag{ 57 }$$

$$\begin{equation}=-\frac{\lambda }{2L}{\int }_{-L}^{L}\mathrm{d}z\enspace {\left({\tilde {P}}_{0}^{\prime }\right)}^{3}.\end{equation} \tag{ 58 }$$

The first of these is most straightforwardly derived from (C.3) in appendix C by using the ODE (56) for ${\tilde {P}}_{0}$, and verifies that ${\dot {\mathcal{S}}}_{-1}{\geqslant}0$ as must be the case. The latter, i.e. (58), follows immediately from (18) after integrating out total derivatives. Now, from the final expression in (58) one observes that ${\dot {\mathcal{S}}}_{-1}$ vanishes identically for even distributions, i.e. those that satisfy ${\tilde {P}}_{0}\left({z}_{0}+z\right)={\tilde {P}}_{0}\left({z}_{0}-z\right)$ for some z₀, which is exactly the P $\mathcal{T}$-symmetry in (24). However, the travelling wave profiles are clearly asymmetric with a steeper front than tail, leading in general to the observed ${\dot {\mathcal{S}}}_{-1}\ne 0$. In figure 3 we illustrate this, and in particular how the banded profile transforms under P $\mathcal{T}$. Finally, a sanity check also verifies that both expressions (57) and (58) are invariant under P alone (${\tilde {P}}_{0}\left(z\right)\to -{\tilde {P}}_{0}\left(-z\right)$) as they should be since $\dot {\mathcal{S}}$ is, as remarked previously, oblivious to whether the wave is moving left or right (in (56) this must be complemented by c → −c).

Zoom In Zoom Out Reset image size

Interestingly, this mode of TRS violation at D = 0 is a collective emergent phenomenon and does not have any counterpart for a single active particle. On the microscopic scale, it is associated with different rates of promotion and demotion of spins at the head and tail of the nonlinear profile respectively, leading to a difference in the rate of dissipation from the active alignment at these edges—a point which we will explore further in a separate paper. In the following section, this point will become more explicit when we see how TRS can also be violated at D = 0 by explicitly tracking mass currents in addition to the local polar density. In particular, in this case there is an analogous mode of TRS violation on the level of a single active particle.

We include in figure 2 the scaling of the EPR $\dot {\mathcal{S}}$ in both the microphase-separated and polar cluster regimes. As shown, both are truly nonequilibrium within our classification scheme, with $\dot {\mathcal{S}}\sim {D}^{-1}$. Although we do not possess explicit polar cluster solutions to the dynamics at D = 0, it is straightforward to argue heuristically that this is what one should expect due to the highly inhomogeneous nature of the phase. For future work, we aim to investigate this in more detail.

So far we have seen that the EPR of the HVM at small noise satisfies the asymptotic scaling relation $\dot {\mathcal{S}}\sim {D}^{\chi }$ for D ≪ 1, where the exponent χ ∈ {−1, 0, 1, ...}. By performing a small noise expansion, we may systematically investigate the coefficients that appear at each order to determine the lowest order nontrivial contribution, and thus χ. We also find that symmetries effectively bound χ from below; when the ground-state dynamics is pathwise equilibrium, the contribution ${\dot {\mathcal{S}}}_{-1}/D=\dot {\mathcal{S}}\left[{\rho }_{0},{\boldsymbol{P}}_{0}\right]$ is locked out and χ > −1. In the isotropic case, the fact that ⟨ρ₁∇ ⋅ P ₁⟩ = 0 further constrains χ > 0, and the small noise limit becomes an effective equilibrium regime. In section 4 we look at the entropy production in a generalised model, where the constitutive equation for the density current in (2) is modified to include a diffusive contribution.

In the following, we add a diffusive contribution and noise to the constitutive equation for the density current J . Specifically, we consider J = J _d + ξ as in [19], where the noise ξ is mean zero, Gaussian and white with covariance

$$\begin{equation}\langle {\xi }_{\alpha }\left(\boldsymbol{x},t\right){\xi }_{\beta }\left({\boldsymbol{x}}^{\prime },{t}^{\prime }\right)\rangle =2{D}_{\rho }{\delta }_{\alpha \beta }\delta \left(\boldsymbol{x}-{\boldsymbol{x}}^{\prime }\right)\delta \left(t-{t}^{\prime }\right).\end{equation} \tag{ 59 }$$

Furthermore, we take the deterministic part J _d of the current J to be of the form

$$\begin{equation}{\boldsymbol{J}}_{\mathrm{d}}=w\boldsymbol{P}-{\gamma }^{-1}\nabla \mu ,\end{equation} \tag{ 60 }$$

where γ is a constant friction coefficient. Here, μ serves an analogous purpose with the chemical potential known from equilibrium diffusive systems. However, since TRS is broken in this model, there is a priori no reason that it should be the functional variation of a free-energy. Notwithstanding, we will for simplicity ignore this issue and assume that we may write

$$\begin{equation}\mu =\frac{\delta F}{\delta \rho },\end{equation} \tag{ 61 }$$

with the same functional F as that which appears in the equation for P . We only make minor changes to F for stability purposes by modifying the local free-energy f(ρ, P ) to include a quadratic term in ρ, so that now

$$\begin{equation}f\left(\rho ,\boldsymbol{P}\right)=\frac{{a}_{\rho }}{2}{\rho }^{2}+\frac{1}{2}\left(1-\rho \right)\vert \boldsymbol{P}{\vert }^{2}+\frac{1}{4}\vert \boldsymbol{P}{\vert }^{4}.\end{equation} \tag{ 62 }$$

In addition we add a square gradient contribution, giving us

$$\begin{equation}F\left[\rho ,\boldsymbol{P}\right]={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left(f\left(\rho ,\boldsymbol{P}\right)+\frac{{\nu }_{\rho }}{2}\vert \nabla \rho {\vert }^{2}+\frac{1}{2}{\left({\nabla }_{\alpha }{P}_{\beta }\right)}^{2}+\boldsymbol{P}\cdot \nabla {\Phi}\left(\rho ,\boldsymbol{P}\right)\right).\end{equation} \tag{ 63 }$$

Observe, however, that these new terms do not change the equation for P in (8), since they both drop out when considering the functional variation of F with respect to P . With this choice, we have that

$$\begin{equation}\mu ={a}_{\rho }\rho -{\nu }_{\rho }{\nabla }^{2}\rho -\frac{1}{2}\vert \boldsymbol{P}{\vert }^{2}-{w}_{1}\nabla \cdot \boldsymbol{P}.\end{equation} \tag{ 64 }$$

Again, (60) is motivated by coarse graining, and the diffusive contribution arises for example in cases where interactions such as steric repulsion are included in the microscopic model [19]. Notably, there is a kind of paradigm shift when breaking the local linear relation J ∝ P , which implies that W = P /ρ should no longer be considered the local average direction of the velocity of particles. Physically, this reflects a situation on the microscopic scale where the bare self-propulsion may be thwarted by e.g. repulsive forces, so that mass currents may move against the local polar order. More importantly in the context of entropy production, this means that a trajectory of the system in which J and P do not point in the same direction is realisable in the forward time ensemble, since fluctuations alone can now reverse J at fixed P even if highly unlikely.

Numerical integration of the dynamics with J = J _d + ξ , hereafter referred to as the DFM, allows us to investigate the resulting phase diagram as in section 2. On the other hand, achieving analytical progress to a comparable extent as with the HVM is more difficult. Notably, however, from a linear analysis we do in fact find a finite wavelength instability in the region where ρ₀ < 1, in which the coarsening dynamics develop a polar crystalline structure as illustrated in figure 4. Our analysis also provides us with the isotropic-to-crystal phase-boundary, and we find that in the case w₁ = w/2 the system is unstable to perturbations when

$$\begin{equation}{w}_{1}^{2}{ >}{w}_{c}^{2}\equiv {a}_{\rho }+{\nu }_{\rho }\left(1-{\rho }_{0}\right)+4\gamma {\nu }_{\rho }+2\sqrt{{\nu }_{\rho }\left(1-{\rho }_{0}+2\gamma \right)\left({a}_{\rho }+2\gamma {\nu }_{\rho }\right)}.\end{equation} \tag{ 65 }$$

More specifically, when the inequality (65) holds there is a finite range q ∈ [q₋, q₊] of modes that are unstable, where the exact expressions for q_± are provided in appendix A. In contrast, we do not observe significant changes to the phase diagram in the region where ρ₀ > 1. We explain this behaviour by observing that when ρ₀ < 1, and the local polar order is weak, the diffusive dynamics is significant while it is overpowered by advective transport when the polar order is strong [19]. At D = 0 the state is stationary and has both ∂_t ρ = 0 and |∂_t P | = 0. In fact, from simulations we observe that the stronger condition |⟨ J _d⟩| = 0 is met, meaning that at D = 0 we expect

$$\begin{equation}w\boldsymbol{P}={\gamma }^{-1}\nabla \mu .\end{equation} \tag{ 66 }$$

From simulations we observe that the local polar order is directed such that it points in towards low density, as illustrated in figure 4. Equation (66) then tells us that the advective transport induced by P is compensated by a reversed 'diffusion' running up gradients in ρ.

Zoom In Zoom Out Reset image size

In the following we also restrict to the case where D_ρ = D/γ for simplicity [19], in which case the Freidlin–Wentzell action for the DFM takes the form

$$\begin{equation}{\mathcal{A}}_{\mathrm{D}\mathrm{F}}\left[\rho ,\boldsymbol{P}\right]=\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left[\gamma {\left\vert {\nabla }^{-1}\left({\partial }_{t}\rho +\nabla \cdot {\boldsymbol{J}}_{\mathrm{d}}\right)\right\vert }^{2}+{\left\vert {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta F}{\delta \boldsymbol{P}}\right\vert }^{2}\right],\end{equation} \tag{ 67 }$$

and the path transition density ${\mathcal{P}}_{\mathrm{D}\mathrm{F}}\left[\rho ,\boldsymbol{P}\right]$ is constructed as before by setting ${\mathcal{P}}_{\mathrm{D}\mathrm{F}}\propto \mathrm{exp}\left(-{\mathcal{A}}_{\mathrm{D}\mathrm{F}}/D\right)$. Note that in (67), we have defined the inverse gradient operator ∇⁻¹ = ∇⁻²∇, i.e. with gauge choice |∇ × ∇⁻¹ h( x )| = 0 [33]. Crucially, with the added density current fluctuations, we are now free to define time-reversal without the polarity flip used in (14). Specifically, we let

$$\begin{equation}{\mathcal{T}}^{{\pm}}=\begin{cases}\rho \left(\boldsymbol{x},t\right){\mapsto}\rho \left(\boldsymbol{x},-t\right)\quad \hfill \\ \boldsymbol{P}\left(\boldsymbol{x},t\right){\mapsto}{\pm}\boldsymbol{P}\left(\boldsymbol{x},-t\right)\quad \hfill \end{cases}\end{equation} \tag{ 68 }$$

and observe that both compositions ${\mathcal{A}}_{\mathrm{D}\mathrm{F}}{\circ}{\mathcal{T}}^{{\pm}}$ are now well defined on the full space of trajectories. As before, this means that when we define the two time-reversed ensembles to ${\mathcal{P}}_{\mathrm{D}\mathrm{F}}$ by setting ${{\leftarrow}{\mathcal{P}}}_{\mathrm{D}\mathrm{F}}^{{\pm}}={\mathcal{P}}_{\mathrm{D}\mathrm{F}}{\circ}{\mathcal{T}}^{{\pm}}$, all trajectories that are observable under ${\mathcal{P}}_{\mathrm{D}\mathrm{F}}$ are also observable under ${{\leftarrow}{\mathcal{P}}}_{\mathrm{D}\mathrm{F}}^{{\pm}}$. Now, comparing the time-forward ensemble with each of these gives rise to two different definitions of the EPR, given by

$$\begin{equation}{\dot {\mathcal{S}}}^{{\pm}}\equiv \underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{1}{2\tau }\enspace \mathrm{log}\enspace \frac{{\mathcal{P}}_{\mathrm{D}\mathrm{F}}\left[\rho ,\boldsymbol{P}\right]}{{{\leftarrow}{\mathcal{P}}}_{\mathrm{D}\mathrm{F}}^{{\pm}}\left[\rho ,\boldsymbol{P}\right]}.\end{equation} \tag{ 69 }$$

In section 4.2 we will attempt to understand how this choice of polar signature may alter the scaling of the EPR at low noise [35, 55].

Analogously with our treatment in the previous section, we observe that when P is odd under time-reversal, the functional F splits into even and odd pieces F^S and F^A respectively. However, in our current setting this has further consequences as well, since it also implies that we should not consider μ a chemical potential like quantity either. Indeed, we see that μ splits into contributions

$$\begin{equation}\mu ={\mu }^{S}+\frac{\delta {F}^{A}}{\delta \rho }.\end{equation} \tag{ 70 }$$

Furthermore, the deterministic part of the current, J _d, also splits into a P -like odd piece under time-reversal and a ∇μ^S -like even piece. That is, we write ${\boldsymbol{J}}_{\mathrm{d}}={\boldsymbol{J}}_{\mathrm{d}}^{S}+{\boldsymbol{J}}_{\mathrm{d}}^{A}$, where we define

$$\begin{equation}{\boldsymbol{J}}_{\mathrm{d}}^{S}=-{\gamma }^{-1}\nabla {\mu }^{S},\end{equation} \tag{ 71 }$$

$$\begin{equation}{\boldsymbol{J}}_{\mathrm{d}}^{A}=w\boldsymbol{P}-{\gamma }^{-1}\nabla \frac{\delta {F}^{A}}{\delta \rho }.\end{equation} \tag{ 72 }$$

Consequently, since F does not remain invariant under time-reversal, and therefore neither μ nor J _d either, it could not feature in an equilibrium theory and violates TRS.

With these definitions, we may again deduce explicit integral expressions for the EPRs ${\dot {\mathcal{S}}}^{{\pm}}$ by using (67)–(69) and the definitions of ${\mathcal{P}}_{\mathrm{D}\mathrm{F}}$, ${{\leftarrow}{\mathcal{P}}}_{\mathrm{D}\mathrm{F}}^{{\pm}}$, and we refer to appendix C for the details. There we show that

$$\begin{equation}{\dot {\mathcal{S}}}^{+}={D}^{-1}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \gamma {w}^{2}\vert K\boldsymbol{P}{\vert }^{2}-w\boldsymbol{P}\cdot \nabla \mu +\lambda \left[\left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}\right]\cdot \left(\lambda \left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}+\frac{\delta F}{\delta \boldsymbol{P}}\right)\right\rangle \end{equation} \tag{ 73 }$$

and

$$\begin{equation}{\dot {\mathcal{S}}}^{-}={D}^{-1}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {\boldsymbol{J}}_{\mathrm{d}}^{A}\cdot \nabla {\mu }^{S}-\left(\lambda \left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right)\cdot \frac{\delta {F}^{S}}{\delta \boldsymbol{P}}\right\rangle .\end{equation} \tag{ 74 }$$

where K is a matrix operator with entries ${K}_{\alpha \beta }={\nabla }_{\alpha }^{-1}{\nabla }_{\beta }={\nabla }^{-2}{\nabla }_{\alpha }{\nabla }_{\beta }$. Note that in (73) we have performed an average over noise histories in order to obtain the given expression, which explains why the symmetry ${\dot {\mathcal{S}}}^{+}{\circ}\mathcal{T}=-{\dot {\mathcal{S}}}^{+}$ does not seem to hold pathwise any longer. However, it is recovered when writing the expression out as in e.g. (C.14) in appendix C. As in section 3, we proceed to analyse (73) and (74) in the low D limit both analytically and numerically. We also carry over the definitions we employed there, in particular defining the exponents χ^± via the asymptotic scaling relation ${\dot {\mathcal{S}}}^{{\pm}}\left(D\right)\sim {D}^{{\chi }^{{\pm}}}$ for D ≪ 1. As we will see, we find that similar considerations to those made before carry over in a straightforward manner, allowing us to predict the correct scaling in all cases.

Continuing as in section 3, we look for ground-state trajectories (ρ₀, P ₀) that satisfy the pathwise equilibrium condition

$$\begin{equation}{\dot {\mathcal{S}}}^{{\pm}}\left[{\rho }_{0},{\boldsymbol{P}}_{0}\right]=0,\end{equation} \tag{ 75 }$$

in order to determine when we should expect χ^± = −1. Again, it is clear that these include all states that are either ${\mathcal{T}}^{{\pm}}$ or P ${\mathcal{T}}^{{\pm}}$-symmetric, and so the situation remains unchanged when choosing ${\mathcal{T}}^{-}$. Indeed, in this case the isotropic P ₀ = 0 state satisfies both symmetries, while the PL $\vert {\boldsymbol{P}}_{0}\vert =\sqrt{{\rho }_{0}-1}$ state is P ${\mathcal{T}}^{-}$-symmetric only. However, we should also expect a similar situation when choosing the ${\mathcal{T}}^{+}$-protocol for time-reversal; since P ₀ does not flip sign upon time-reversal, any constant trajectory satisfies ${\mathcal{T}}^{+}$ alone. Thus we expect χ^± > −1 for both the isotropic gas and PL, meaning that there is no clear distinction between the two protocols for the homogeneous phases at ground-state level.

On the other hand, the situation changes quite drastically once the density current dynamics are tracked explicitly. In particular, in doing so, we expect that TRS violation at the D = 0 level should become visible from a misalignment of the density current and polar density. To see this, we promote J to a dynamical variable and consider the Freidlin–Wentzell action at this level. That is, we define

$$\begin{equation}{\mathcal{A}}_{\mathrm{D}\mathrm{F}}^{J}\left[\rho ,\boldsymbol{J},\boldsymbol{P}\right]=\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left[\gamma \vert \boldsymbol{J}-{\boldsymbol{J}}_{\mathrm{d}}{\vert }^{2}+{\left\vert {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta F}{\delta \boldsymbol{P}}\right\vert }^{2}\right],\end{equation} \tag{ 76 }$$

if ∂_t ρ + ∇ ⋅ J = 0 and ${\mathcal{A}}_{\mathrm{D}\mathrm{F}}^{J}=\infty$ otherwise. Importantly, J now takes the role that w P had previously in section 3, in the sense that it must be odd on time-reversal. Again, this is due to the constraint imposed by the continuity equation which limits the space of observable trajectories under the action (76). Time-reversal is thus generalised accordingly by setting

$$\begin{equation}{\mathcal{T}}_{J}^{{\pm}}=\begin{cases}\rho \left(\boldsymbol{x},t\right){\mapsto}\rho \left(\boldsymbol{x},-t\right),\quad \hfill \\ \boldsymbol{J}\left(\boldsymbol{x},t\right){\mapsto}-\boldsymbol{J}\left(\boldsymbol{x},-t\right),\quad \hfill \\ \boldsymbol{P}\left(\boldsymbol{x},t\right){\mapsto}{\pm}\boldsymbol{P}\left(\boldsymbol{x},-t\right).\quad \hfill \end{cases}\end{equation} \tag{ 77 }$$

Pathwise there is now a clear distinction between the two protocols ${\mathcal{T}}_{J}^{{\pm}}$. Indeed, when (ρ₀, J ₀, P ₀) is a constant trajectory with both | J ₀| > 0 and | P ₀| > 0, the protocol ${\mathcal{T}}_{J}^{+}$ introduces a discrepancy between J ₀ and P ₀ that cannot be transformed away by parity. On the other hand, since both J ₀ and P ₀ transform the same way under ${\mathcal{T}}_{J}^{-}$, the trajectory (ρ₀, J ₀, P ₀) remains invariant under P ${\mathcal{T}}_{J}^{-}$ as before. This fact is reflected by the EPR computed at the level of the action in (76). To see this, we set ${\mathcal{P}}_{\mathrm{D}\mathrm{F}}^{J}\propto \mathrm{exp}\left(-{\mathcal{A}}_{\mathrm{D}\mathrm{F}}^{J}/D\right)$ and ${{\leftarrow}{\mathcal{P}}}_{\mathrm{D}\mathrm{F}}^{J,{\pm}}={\mathcal{P}}_{\mathrm{D}\mathrm{F}}^{J}{\circ}{\mathcal{T}}_{J}^{{\pm}}$ as before, and define ${\dot {\mathcal{S}}}_{J}^{{\pm}}$ via the usual definition as in (69). In appendix C we show that this computation results in the expression

$$\begin{equation}{\dot {\mathcal{S}}}_{J}^{+}={D}^{-1}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left\langle \gamma {w}^{2}\vert \boldsymbol{P}{\vert }^{2}-w\boldsymbol{P}\cdot \nabla \mu -\lambda \left[\left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}\right]\cdot \left(\lambda \left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}+\frac{\delta F}{\delta \boldsymbol{P}}\right)\right\rangle \end{equation} \tag{ 78 }$$

for ${\dot {\mathcal{S}}}_{J}^{+}$, which differs from (73) by omission of the factor K inside the first term, while in fact ${\dot {\mathcal{S}}}_{J}^{-}={\dot {\mathcal{S}}}^{-}$ remains unchanged from (74).

It is instructive to investigate the difference between the two expressions (73) for ${\dot {\mathcal{S}}}^{+}$ and (78) for ${\dot {\mathcal{S}}}_{J}^{+}$, and in particular to show that indeed the operator K ≠ Id. To this end, we introduce the potential φ as the solution to the Poisson's equation

$$\begin{equation}{\nabla }^{2}\varphi =\nabla \cdot \boldsymbol{P},\end{equation} \tag{ 79 }$$

which is unique up to an additive constant (assuming periodic boundaries on $\mathcal{V}$). Defining P _T ≡ P − ∇φ it follows by construction that

$$\begin{equation}\boldsymbol{P}=\nabla \varphi +{\boldsymbol{P}}_{\mathrm{T}},\end{equation} \tag{ 80 }$$

where P _T is solenoidal, i.e. ∇ ⋅ P _T = 0. Importantly, this construction is in general not the same as the standard Helmholtz decomposition, since P _T is not necessarily the curl of a vector potential. A specific example demonstrating that these are indeed different is provided via simplest case for which P = P ₀ is constant and nonzero. Indeed, in this case, periodic boundaries forces φ = const. and P _T = P ₀, and the latter cannot be written as the curl of a vector field that respects the periodic boundaries. The decomposition in (80) is, however, orthogonal in ${L}^{2}\left(\mathcal{V}\right)$, meaning that

$$\begin{equation}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \vert \boldsymbol{P}{\vert }^{2}={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \vert \nabla \varphi {\vert }^{2}+{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \vert {\boldsymbol{P}}_{\mathrm{T}}{\vert }^{2}.\end{equation} \tag{ 81 }$$

Thus, observing that by definition we have K P = ∇φ, it follows that the difference between ${\dot {\mathcal{S}}}^{+}$ and ${\dot {\mathcal{S}}}_{J}^{+}$ is simply

$$\begin{equation}{\dot {\mathcal{S}}}_{J}^{+}-{\dot {\mathcal{S}}}^{+}=\gamma {w}^{2}{D}^{-1}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \vert {\boldsymbol{P}}_{\mathrm{T}}{\vert }^{2}\right\rangle .\end{equation} \tag{ 82 }$$

This is consistent with our observation that the PL breaks ${\mathcal{T}}_{J}^{+}$ at ground-state level as remarked above, and in particular it follows that we have

$$\begin{equation}\underset{D\to 0}{\mathrm{lim}}\enspace D{\dot {\mathcal{S}}}_{J}^{+}/\mathcal{V}=\gamma {w}^{2}\vert {\boldsymbol{P}}_{0}{\vert }^{2}\end{equation} \tag{ 83 }$$

for the constant homogeneous ground-states. This mechanism by which TRS is broken at D = 0 due to J and P having different polar signatures under time-reversal has an analogue for a single active particle on the microscopic level. To see this we make the identifications ${\dot {\boldsymbol{x}}}_{t}\enspace \delta \left({\boldsymbol{x}}_{t}-\boldsymbol{x}\right)\to \boldsymbol{J}$ and ${\hat{\boldsymbol{n}}}_{t}\enspace \delta \left({\boldsymbol{x}}_{t}-\boldsymbol{x}\right)\to \boldsymbol{P}$, where x _t is the position of the active particle and ${\hat{\boldsymbol{n}}}_{t}$ denotes its polar orientation and direction of self-propulsion. Similarly to the current J , the particle velocity ${\dot {\boldsymbol{x}}}_{t}$ must change sign under time-reversal. The polar orientation ${\hat{\boldsymbol{n}}}_{t}$ need not, on the other hand, and thus generally leads to a bare entropy production associated with the motility.

Continuing as in section 3, we include fluctuations by performing a systematic expansion of the equations of motion and the EPRs ${\dot {\mathcal{S}}}^{{\pm}}$ via the integral expressions in (73) and (74) in small D. Again, we start by assuming that the ground-state (ρ₀, P ₀) is constant and homogeneous. Thus, substituting the expansions in (25) and (26) into the continuity equation (1) with J = J _d + ξ we obtain to lowest nontrivial order

$$\begin{equation}{\partial }_{t}{\rho }_{1}=-\nabla \cdot \left(w{\boldsymbol{P}}_{1}-{\gamma }^{-1}\nabla \frac{\delta {F}_{\mathrm{L}}}{\delta \rho }\left[{\rho }_{1},{\boldsymbol{P}}_{1}\right]+{\boldsymbol{\xi }}_{1}\right),\end{equation} \tag{ 84 }$$

where now with slight abuse of notation

$$\begin{equation}{F}_{\mathrm{L}}\left[\rho ,\boldsymbol{P}\right]={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left({f}_{\mathrm{L}}\left(\rho ,\boldsymbol{P}\right)+\frac{{\nu }_{\rho }}{2}\vert \nabla \rho {\vert }^{2}+\frac{1}{2}{\left({\nabla }_{\alpha }{P}_{\beta }\right)}^{2}+\boldsymbol{P}\cdot \nabla {{\Phi}}_{\mathrm{L}}\left(\rho ,\boldsymbol{P}\right)\right).\end{equation} \tag{ 85 }$$

In addition, the local free energy f_L is given by

$$\begin{equation}{f}_{\mathrm{L}}\left(\rho ,\boldsymbol{P}\right)=\frac{{a}_{\rho }}{2}{\rho }^{2}+\frac{{a}_{0}}{2}\vert \boldsymbol{P}{\vert }^{2}-\rho {\boldsymbol{P}}_{0}\cdot \boldsymbol{P}+{\left({\boldsymbol{P}}_{0}\cdot \boldsymbol{P}\right)}^{2},\end{equation} \tag{ 86 }$$

while Φ_L remains the same as in (31). The noise term ξ ₁ is a mean zero Gaussian white noise process with covariance

$$\begin{equation}\langle {\xi }_{1\alpha }\left(\boldsymbol{x},t\right){\xi }_{1\beta }\left({\boldsymbol{x}}^{\prime },{t}^{\prime }\right)\rangle =2{\gamma }^{-1}{\delta }_{\alpha \beta }\delta \left(\boldsymbol{x}-{\boldsymbol{x}}^{\prime }\right)\delta \left(t-{t}^{\prime }\right),\end{equation} \tag{ 87 }$$

so that the linearised equation (84) is indeed independent of D. Note also that there is no change to the linearised equation for the polar density since this is the same for both models under consideration.

Similarly, an expansion of the EPRs ${\dot {\mathcal{S}}}^{{\pm}}$ in small D allows us to write

$$\begin{equation}{\dot {\mathcal{S}}}^{{\pm}}\left(D\right)={\dot {\mathcal{S}}}_{-1}^{{\pm}}{D}^{-1}+{\dot {\mathcal{S}}}_{0}^{{\pm}}+{\dot {\mathcal{S}}}_{1}^{{\pm}}D+O\left({D}^{2}\right),\end{equation} \tag{ 88 }$$

where ${\dot {\mathcal{S}}}_{k}^{{\pm}}=0$ for k < χ^±, and ${\dot {\mathcal{S}}}_{{\chi }^{{\pm}}}^{{\pm}}{ >}0$ as before. By explicitly computing this expansion one finds that

$$\begin{equation}{\chi }^{{\pm}}=\begin{cases}0,\quad \text{isotropic},\quad \hfill \\ 0,\quad \text{polar}\;\text{liquid}.\quad \hfill \end{cases}\end{equation} \tag{ 89 }$$

In particular, since the linearised continuity equation (84) of the DFM implies that the steady-state expectation ⟨ρ₁∇ ⋅ P ₁⟩ no longer vanishes identically as in (37), we cannot any longer expect that χ⁻ > 0 in the isotropic phase. Since χ^± = 0 in the isotropic and PL phases, we classify both phases as being marginally nonequilibrium for the DFM. This means that the linearised dynamics of the DFM cannot be mapped onto an equilibrium dynamics for any choice of P ₀.

As in section 3, we now present a more in-depth analysis of the leading order term in the expansion (88) of ${\dot {\mathcal{S}}}^{{\pm}}$. Beginning with the isotropic phase where | P ₀| = 0, we find after a Fourier transform that ${\dot {\mathcal{S}}}_{0}^{{\pm}}$ may be written in bilinear form as in (49):

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{{\pm}}\left({\Lambda}\right)/\mathcal{V}=\sum _{0{< }\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\left\langle {\left({\hat{\boldsymbol{u}}}^{\text{iso}}\right)}^{{\dagger}}{\dot {\sigma }}^{{\pm},\mathrm{i}\mathrm{s}\mathrm{o}}{\hat{\boldsymbol{u}}}^{\text{iso}}\right\rangle =\sum _{0{< }\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\mathrm{Tr}\left({\dot {\sigma }}^{{\pm},\mathrm{i}\mathrm{s}\mathrm{o}}{\mathcal{C}}^{\text{iso}}\right).\end{equation} \tag{ 90 }$$

Here, the sum runs over wavevectors q , and an explicit dependence in ${\dot {\mathcal{S}}}_{0}^{{\pm}}\left({\Lambda}\right)$ on the ultraviolet cutoff Λ is introduced in order to study the limit in which it is taken to infinity. We denote the Fourier modes of ρ₁ and P ₁ by ${\hat{\rho }}_{1}$ and ${\hat{\boldsymbol{P}}}_{1}$ respectively, and have defined

$$\begin{equation}{\hat{\boldsymbol{u}}}^{\text{iso}}={\left({\hat{\rho }}_{1},{\hat{P}}_{\mathrm{L}},{\hat{P}}_{\mathrm{T}}\right)}^{\mathrm{T}},\end{equation} \tag{ 91 }$$

where ${\hat{P}}_{\mathrm{L}}=\hat{\boldsymbol{q}}\cdot {\boldsymbol{P}}_{1}$ and ${\hat{P}}_{\mathrm{T}}={\hat{\boldsymbol{q}}}_{\perp }\cdot {\boldsymbol{P}}_{1}$ are respectively the longitudinal and transverse components of P ₁ with respect to $\hat{\boldsymbol{q}}$, and ${\hat{\boldsymbol{q}}}_{\perp }$ is perpendicular to q and of unit length. The equal-time steady-state correlation matrix ${\mathcal{C}}^{\text{iso}}\equiv \left({\mathcal{C}}_{ij}^{\text{iso}}\right)$ is then defined by

$$\begin{equation}{\mathcal{C}}_{ij}^{\text{iso}}\left(q\right){\delta }_{\boldsymbol{q},{\boldsymbol{q}}^{\prime }}=\left\langle {\hat{u}}_{i}^{\text{iso}}\left(\boldsymbol{q},t\right){\left({\hat{u}}_{j}^{\text{iso}}\left({\boldsymbol{q}}^{\prime },t\right)\right)}^{{\ast}}\right\rangle .\end{equation} \tag{ 92 }$$

In addition, the Hermitian matrices ${\dot {\sigma }}^{{\pm},\mathrm{i}\mathrm{s}\mathrm{o}}$ in (90) are given by

$$\begin{equation}{\sigma }^{+,\mathrm{i}\mathrm{s}\mathrm{o}}=\frac{1}{2}\left(\begin{matrix}\hfill 0\hfill & \hfill iwq{{\Gamma}}_{\rho }\hfill & \hfill 0\hfill \\ \hfill -iwq{{\Gamma}}_{\rho }\hfill & \hfill 2\gamma w\tilde {w}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 93 }$$

and

$$\begin{equation}{\sigma }^{-,\mathrm{i}\mathrm{s}\mathrm{o}}=\frac{1}{2}\left(\begin{matrix}\hfill 0\hfill & \hfill iq\left({w}_{1}{\Gamma}-\tilde {w}{{\Gamma}}_{\rho }\right)\hfill & \hfill 0\hfill \\ \hfill -iq\left({w}_{1}{\Gamma}-\tilde {w}{{\Gamma}}_{\rho }\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 94 }$$

where we have defined $\tilde {w}=w\left(1-{w}_{1}{q}^{2}/\gamma w\right)$ as well as damping coefficients Γ = 1 − ρ₀ + q² and Γ_ρ = a_ρ + ν_ρ q². In fact, we may explicitly compute ${\mathcal{C}}^{\text{iso}}$ from the linearised dynamics and we refer to appendix A for the details. By substituting the result of this calculation back into (90), we obtain

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{+}\left({\Lambda}\right)=\sum _{0{< }\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\frac{\gamma {w}^{2}}{{\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }},\end{equation} \tag{ 95 }$$

and

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{-}\left({\Lambda}\right)=\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\frac{{\gamma }^{-1}{q}^{2}{\left({w}_{1}{\Gamma}-\tilde {w}{{\Gamma}}_{\rho }\right)}^{2}}{\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)\left(\tilde {w}{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }{\Gamma}\right)}.\end{equation} \tag{ 96 }$$

Interestingly, we see from direct power counting that ${\dot {\mathcal{S}}}_{0}^{+}\left({\Lambda}\right)$ converges as Λ → ∞, while ${\dot {\mathcal{S}}}_{0}^{-}\left({\Lambda}\right)\sim \mathrm{log}\enspace {\Lambda}$. In fact, expressions (95) and (96) generalise trivially to dimensions d ≠ 2, meaning that

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{+}\left({\Lambda}\right)\sim \begin{cases}1,\quad \hfill & d{< }4\hfill \\ {{\Lambda}}^{d-4},\quad \hfill & d{\geqslant}4\hfill \end{cases}\end{equation} \tag{ 97 }$$

and

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{-}\left({\Lambda}\right)\sim \begin{cases}1,\quad \hfill & d{< }2\hfill \\ {{\Lambda}}^{d-2},\quad \hfill & d{\geqslant}2\hfill \end{cases},\end{equation} \tag{ 98 }$$

where we denote by Λ⁰ a logarithmic divergence.

We may perform an identical procedure in the PL case, although we leave the details of this calculation in appendix C to simplify the presentation. We note, however, that in the PL phase the scaling of ${\dot {\mathcal{S}}}_{0}^{{\pm}}\left({\Lambda}\right)$ with the ultraviolet cutoff Λ changes. For our present case where d = 2, we find that ${\dot {\mathcal{S}}}_{0}^{+}\sim {{\Lambda}}^{2}$ while ${\dot {\mathcal{S}}}_{0}^{-}\sim {{\Lambda}}^{4}$.

Similarly to our treatment in section 3, we find good agreement between predictions and the results from simulations. In figure 5 we demonstrate this comparison for both of the homogeneous phases. Furthermore, all considerations extend straightforwardly to the level where we explicitly track the current J ; the scaling exponent ${\chi }_{J}^{+}$ of the EPR ${\dot {\mathcal{S}}}_{J}^{+}$ is at this level given by

$$\begin{equation}{\chi }_{J}^{+}=\begin{cases}0,\quad \hfill & \text{isotropic}\hfill \\ -1,\quad \hfill & \text{polar}\;\text{liquid}\hfill \end{cases}\end{equation} \tag{ 99 }$$

while ${\chi }_{J}^{-}={\chi }^{-}$. Note, however, that in the isotropic phase the coefficient of the leading order term of ${\dot {\mathcal{S}}}_{J}^{+}$ changes by virtue of the equation (82).

Zoom In Zoom Out Reset image size

Finally, we consider inhomogeneous ground-states (ρ₀, P ₀) (or (ρ₀, J ₀, P ₀) at the level of J ), specifically the nonlinear polar cluster, microphase-separated and polar crystal states. Following the same reasoning as in section 3.2, we conclude that ${\chi }_{J}^{{\pm}}={\chi }^{{\pm}}=-1$ for both the banded profiles and PC. Indeed, these profiles still break both ${\mathcal{T}}^{{\pm}}$ and P ${\mathcal{T}}^{{\pm}}$ as before and are therefore truly nonequilibrium at small noise. On the other hand, for the crystal state we find that χ⁺ = 0, while χ⁻ = −1, as shown in figure 5. We explain this by the observation that at D = 0 the ground-state solution is stationary, i.e. it has both ∂_t ρ₀ = 0 and ∂_t P ₀ = 0. In particular, since (ρ₀, P ₀) is independent of time, it is in fact also invariant under ${\mathcal{T}}^{+}$. This may also be confirmed by inspection of e.g. (C.14), where it is apparent that the stationary condition implies we must have χ⁺ > −1. Similarly, at the level of the current J we also find that ${\chi }_{J}^{+}=0$ for the crystal state. To see why this should be true, note that also | J ₀| = 0 for the polar crystal at ground-state level, meaning that there is no difference between (ρ₀, J ₀, P ₀) and its time-reversal under ${\mathcal{T}}_{J}^{+}$ for this phase. Coincidentally, this also shows that inhomogeneity is not necessarily sufficient alone to make the system truly non-equilibrium. On the other hand, the polar crystal state is clearly not invariant under ${\mathcal{T}}^{-}$ or P ${\mathcal{T}}^{-}$ (nor ${\mathcal{T}}_{J}^{-}$, P ${\mathcal{T}}_{J}^{-}$), and so we conclude that ${\chi }_{J}^{-}={\chi }^{-}=-1$.

Beginning with the isotropic state, we transform the linearised equations (A.3) and (A.6) for the DFM to Fourier space when n = 1. The resulting equations are most conveniently expressed in matrix form as

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill {\hat{\rho }}_{1}\hfill \\ \hfill {\hat{P}}_{\mathrm{L}}\hfill \\ \hfill {\hat{P}}_{\mathrm{T}}\hfill \end{matrix}\right)=-\underset{{\mathcal{L}}^{\text{iso}}\left(q\right)}{\underbrace{\left(\begin{matrix}\hfill {q}^{2}{{\Gamma}}_{\rho }/\gamma \hfill & \hfill i\tilde {w}q\hfill & \hfill 0\hfill \\ \hfill i{w}_{1}q\hfill & \hfill {\Gamma}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\Gamma}\hfill \end{matrix}\right)}}\left(\begin{matrix}\hfill {\hat{\rho }}_{1}\hfill \\ \hfill {\hat{P}}_{\mathrm{L}}\hfill \\ \hfill {\hat{P}}_{\mathrm{T}}\hfill \end{matrix}\right)+\left(\begin{matrix}\hfill -iq{\hat{\xi }}_{\mathrm{L}}\hfill \\ \hfill {\hat{\eta }}_{\mathrm{L}}\hfill \\ \hfill {\hat{\eta }}_{\mathrm{T}}\hfill \end{matrix}\right).\end{equation} \tag{ A.8 }$$

Here, ${\hat{\xi }}_{\mathrm{L}}=\hat{\boldsymbol{q}}\cdot {\hat{\boldsymbol{\xi }}}_{1}$ is the longitudinal component of the Fourier coefficient ${\hat{\boldsymbol{\xi }}}_{1}$, while ${\hat{\eta }}_{\mathrm{L}}=\hat{\boldsymbol{q}}\cdot {\hat{\boldsymbol{\eta }}}_{1}$, ${\hat{\eta }}_{\mathrm{T}}={\hat{\boldsymbol{q}}}_{\perp }\cdot {\hat{\boldsymbol{\eta }}}_{1}$ are the longitudinal and transverse components of ${\hat{\boldsymbol{\eta }}}_{1}$ respectively. Because of the convention (38) we have chosen for Fourier transforms, noise correlations in Fourier space contain an extra factor of ${\mathcal{V}}^{-1}$, i.e.

$$\begin{equation}\langle {\hat{\xi }}_{\mathrm{L}}\left(\boldsymbol{q},t\right){\hat{\xi }}_{\mathrm{L}}^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =2{\left(\gamma \mathcal{V}\right)}^{-1}{\delta }_{\boldsymbol{q},{\boldsymbol{q}}^{\prime }}\delta \left(t-{t}^{\prime }\right),\end{equation} \tag{ A.9 }$$

$$\begin{equation}\langle {\hat{\eta }}_{\mathrm{L}}\left(\boldsymbol{q},t\right){\hat{\eta }}_{\mathrm{L}}^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =\langle {\hat{\eta }}_{\mathrm{T}}\left(\boldsymbol{q},t\right){\hat{\eta }}_{\mathrm{T}}^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =2{\mathcal{V}}^{-1}{\delta }_{\boldsymbol{q},{\boldsymbol{q}}^{\prime }}\delta \left(t-{t}^{\prime }\right),\end{equation} \tag{ A.10 }$$

while $\langle {\hat{\eta }}_{\mathrm{L}}\left(\boldsymbol{q},t\right){\hat{\eta }}_{\mathrm{T}}^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =\langle {\hat{\eta }}_{\mathrm{L}}\left(\boldsymbol{q},t\right){\hat{\xi }}_{\mathrm{L}}^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =\langle {\hat{\eta }}_{\mathrm{T}}\left(\boldsymbol{q},t\right){\hat{\xi }}_{\mathrm{L}}^{{\ast}}\left({\boldsymbol{q}}^{\prime },{t}^{\prime }\right)\rangle =0$. Note also that the damping coefficient Γ(q) = 1 − ρ₀ + q² > 0 ensures that transverse fluctuations ${\hat{\boldsymbol{P}}}_{\mathrm{T}}$ decay on non-hydrodynamic timescales when ρ₀ < 1.

Observe that solutions to (A.8) are stable and decay at an exponential rate when the eigenvalues σ_±, σ_T of the linear operator ${\mathcal{L}}^{\text{iso}}$ have positive real parts. These are straightforwardly found from the characteristic equation of ${\mathcal{L}}^{\text{iso}}$, and are given by

$$\begin{equation}{\sigma }_{{\pm}}=\frac{1}{2}\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }{\pm}\sqrt{{\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)}^{2}-4{q}^{2}\left(\tilde {w}{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }{\Gamma}\right)}\right),\end{equation} \tag{ A.11 }$$

$$\begin{equation}{\sigma }_{\mathrm{T}}={\Gamma}.\end{equation} \tag{ A.12 }$$

Equation (A.11) for σ_± may be unravelled by first understanding its behaviour at γ = ∞. In this limit, we find that

$$\begin{equation}{\sigma }_{{\pm}}\left(\gamma \to \infty \right)=\frac{1}{2}\left({\Gamma}{\pm}\sqrt{{{\Gamma}}^{2}-4{q}^{2}w{w}_{1}}\right).\end{equation} \tag{ A.13 }$$

Since w, w₁ > 0 it follows that Re σ_±(γ → ∞) > 0 if and only if σ_T = Γ > 0. Thus, at γ = ∞, the isotropic state is linearly stable when ρ₀ < 1, while it is unstable otherwise.

From this, it is fairly straightforward to see that a similar condition holds for finite γ. Specifically, Re σ_± > 0 only when ρ₀ < 1. However, we see from (A.11) that in this case stability also requires that

$$\begin{equation}\tilde {w}{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }{\Gamma}{ >}0.\end{equation} \tag{ A.14 }$$

After some algebra, one finds that this latter condition holds for all q if and only if

$$\begin{equation}{w}_{1}^{2}{< }{a}_{\rho }+{\nu }_{\rho }\left(1-{\rho }_{0}\right)+4\gamma {\nu }_{\rho }+2\sqrt{{\nu }_{\rho }\left(1-{\rho }_{0}+2\gamma \right)\left({a}_{\rho }+2\gamma {\nu }_{\rho }\right)}.\end{equation} \tag{ A.15 }$$

Thus, the phase diagram of the DFM in the region where ρ₀ < 1 as predicted by the linear theory is no longer trivial. Indeed, when condition (A.15) is broken, there is a finite range of wave numbers q ∈ [q₋, q₊] for which the corresponding modes ${\hat{\rho }}_{1}$, ${\hat{P}}_{\mathrm{L}}$ grow in time, where

$$\begin{equation}2{\nu }_{\rho }{q}_{{\pm}}^{2}={w}_{1}^{2}-{\alpha }_{\rho }-{\nu }_{\rho }\left(1-{\rho }_{0}\right)\end{equation} \tag{ A.16 }$$

$$\begin{equation}\quad {\pm}\sqrt{{\left({w}_{1}^{2}-{\alpha }_{\rho }-{\nu }_{\rho }\left(1-{\rho }_{0}\right)\right)}^{2}-4{\nu }_{\rho }\left(\gamma w{w}_{1}+{\alpha }_{\rho }\left(1-{\rho }_{0}\right)\right)}.\end{equation} \tag{ A.17 }$$

From simulations, we find that this instability leads to the polar crystal phase reported in section 4.

We are also interested in calculating the equal-time correlation functions of the linear theory in order to make analytical predictions about the EPR. Since the dynamics is linear, the correlation matrix ${\mathcal{C}}^{\text{iso}}$ in (92) solves the algebraic Riccati equation

$$\begin{equation}{\mathcal{L}}^{\text{iso}}{\mathcal{C}}^{\text{iso}}+{\mathcal{C}}^{\text{iso}}{\left({\mathcal{L}}^{\text{iso}}\right)}^{{\dagger}}=2\mathcal{D},\end{equation} \tag{ A.18 }$$

where we have defined the diffusion matrix $\mathcal{D}$ by

$$\begin{equation}\mathcal{D}=\frac{1}{\mathcal{V}}\left(\begin{matrix}\hfill {q}^{2}/\gamma \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right).\end{equation} \tag{ A.19 }$$

To derive this, one simply needs to apply the chain rule to the left-hand side of

$$\begin{equation}\left\langle \frac{\mathrm{d}}{\mathrm{d}t}\left({\hat{u}}_{i}^{\text{iso}}\left(\boldsymbol{q},t\right){\left({\hat{u}}_{j}^{\text{iso}}\left({\boldsymbol{q}}^{\prime },t\right)\right)}^{{\ast}}\right)\right\rangle =0,\end{equation} \tag{ A.20 }$$

where ${\hat{\boldsymbol{u}}}^{\text{iso}}$ is defined in (51). It is well known that the solution to the Riccati equation (A.18) may be expressed in integral form. However, we find that for our present purposes it is less cumbersome to tackle it straight on. First we observe that the linear system in (A.8) reduces to the two-dimensional coupled dynamics of $\left({\hat{\rho }}_{1},{\hat{P}}_{\mathrm{L}}\right)$, in addition to the one-dimensional dynamics of ${\hat{P}}_{\mathrm{T}}$. Thus, clearly, we may treat these separately. Beginning with the former, the Riccati equation (A.18) may be re-expressed as a four-by-four linear system, specifically

$$\begin{equation}\underset{{\mathcal{R}}^{\text{iso}}\left(q\right)}{\underbrace{\left(\begin{matrix}\hfill 2{q}^{2}{{\Gamma}}_{\rho }/\gamma \hfill & \hfill -i\tilde {w}q\hfill & \hfill i\tilde {w}q\hfill & \hfill 0\hfill \\ \hfill -i{w}_{1}q\hfill & \hfill {\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }+{\Gamma}\hfill & \hfill 0\hfill & \hfill i\tilde {w}q\hfill \\ \hfill i{w}_{1}q\hfill & \hfill 0\hfill & \hfill {\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }+{\Gamma}\hfill & \hfill -i\tilde {w}q\hfill \\ \hfill 0\hfill & \hfill i{w}_{1}q\hfill & \hfill -i{w}_{1}q\hfill & \hfill 2{\Gamma}\hfill \end{matrix}\right)}}\left(\begin{matrix}\hfill {\mathcal{C}}_{11}^{\text{iso}}\hfill \\ \hfill {\mathcal{C}}_{12}^{\text{iso}}\hfill \\ \hfill {\mathcal{C}}_{21}^{\text{iso}}\hfill \\ \hfill {\mathcal{C}}_{22}^{\text{iso}}\hfill \end{matrix}\right)=\frac{2}{\mathcal{V}}\left(\begin{matrix}\hfill {q}^{2}/\gamma \hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{matrix}\right).\end{equation} \tag{ A.21 }$$

To find the correlators of the linear theory we therefore simply invert the matrix ${\mathcal{R}}^{\text{iso}}$, and one may check that the solution is given by

$$\begin{equation}{\mathcal{C}}_{11}^{\text{iso}}=\frac{1}{\mathcal{V}}\frac{\tilde {w}w+{\gamma }^{-1}{\Gamma}\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)}{\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)\left(\tilde {w}{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }{\Gamma}\right)},\end{equation} \tag{ A.22 }$$

$$\begin{equation}{\mathcal{C}}_{12}^{\text{iso}}={\left({C}_{21}^{\text{iso}}\right)}^{{\ast}}=\frac{1}{\mathcal{V}}\frac{i{\gamma }^{-1}q\left({w}_{1}{\Gamma}-\tilde {w}{{\Gamma}}_{\rho }\right)}{\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)\left(\tilde {w}{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }{\Gamma}\right)},\end{equation} \tag{ A.23 }$$

$$\begin{equation}{\mathcal{C}}_{22}^{\text{iso}}=\frac{1}{\mathcal{V}}\frac{w{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)}{\left({\Gamma}+{\gamma }^{-1}{q}^{2}{{\Gamma}}_{\rho }\right)\left(\tilde {w}{w}_{1}+{\gamma }^{-1}{{\Gamma}}_{\rho }{\Gamma}\right)}.\end{equation} \tag{ A.24 }$$

The final non-trivial component of ${\mathcal{C}}^{\text{iso}}$ is found straightforwardly from the dynamics of ${\hat{P}}_{\mathrm{T}}$, and is given by

$$\begin{equation}{\mathcal{C}}_{33}^{\text{iso}}=\frac{1}{\mathcal{V}{\Gamma}}.\end{equation} \tag{ A.25 }$$

To recover the correlators in the linearised HVM, we simply take the limit γ → ∞ in (A.22)–(A.25). We obtain that

$$\begin{equation}{\mathcal{C}}^{\text{iso}}\left(q,\gamma \to \infty \right)=\frac{1}{\mathcal{V}{\Gamma}}\left(\begin{matrix}\hfill w/{w}_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right).\end{equation} \tag{ A.26 }$$

In particular, this result verifies that $\langle {\hat{\rho }}_{1}{\hat{P}}_{\mathrm{L}}^{{\ast}}\rangle =0$ as advertised in (37).

The linearised dynamics about the homogeneous PL phase may be treated similarly to the isotropic case considered above. Transforming (A.3) and (A.6) for n = 1 and | P ₀|² = ρ₀ − 1 to Fourier space we find that

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill {\hat{\rho }}_{1}\hfill \\ \hfill {\hat{P}}_{{\Vert}}\hfill \\ \hfill {\hat{P}}_{\perp }\hfill \end{matrix}\right)=-{\mathcal{L}}^{\text{pl}}\left(\boldsymbol{q}\right)\left(\begin{matrix}\hfill {\hat{\rho }}_{1}\hfill \\ \hfill {\hat{P}}_{{\Vert}}\hfill \\ \hfill {\hat{P}}_{\perp }\hfill \end{matrix}\right)+\left(\begin{matrix}\hfill -iq{\hat{\xi }}_{\mathrm{L}}\hfill \\ \hfill {\hat{\eta }}_{{\Vert}}\hfill \\ \hfill {\hat{\eta }}_{\perp }\hfill \end{matrix}\right),\end{equation} \tag{ A.27 }$$

where we have defined

$$\begin{equation}{\mathcal{L}}^{\text{pl}}\left(\boldsymbol{q}\right)=\left(\begin{matrix}\hfill {q}^{2}{{\Gamma}}_{\rho }/\gamma \hfill & \hfill i\tilde {w}{q}_{{\Vert}}-{q}^{2}{P}_{0}/\gamma \hfill & \hfill i\tilde {w}{q}_{\perp }\hfill \\ \hfill i{w}_{1}{q}_{{\Vert}}-{P}_{0}\hfill & \hfill {{\Gamma}}_{{\Vert}}+i\lambda {P}_{0}{q}_{{\Vert}}\hfill & \hfill i\kappa {P}_{0}{q}_{\perp }\hfill \\ \hfill i{w}_{1}{q}_{\perp }\hfill & \hfill -i\kappa {P}_{0}{q}_{\perp }\hfill & \hfill {{\Gamma}}_{\perp }+i\lambda {P}_{0}{q}_{{\Vert}}\hfill \end{matrix}\right),\end{equation} \tag{ A.28 }$$

in addition to new damping coefficients

$$\begin{equation}{{\Gamma}}_{{\Vert}}=2\left({\rho }_{0}-1\right)+{q}^{2},\end{equation} \tag{ A.29 }$$

$$\begin{equation}{{\Gamma}}_{\perp }={q}^{2}.\end{equation} \tag{ A.30 }$$

As before, ${\hat{\xi }}_{\mathrm{L}}=\hat{\boldsymbol{q}}\cdot \hat{\boldsymbol{\xi }}$ is the longitudinal component of the noise $\hat{\boldsymbol{\xi }}$, while η_∥ and η_⊥ are respectively the parallel and perpendicular components of η ₁ with respect to P ₀. Furthermore, all noise terms ${\hat{\xi }}_{\mathrm{L}}$, ${\hat{\eta }}_{{\Vert}}$ and ${\hat{\eta }}_{\perp }$ are independent and of the same covariance as in the isotropic case.

Rather than solving the full cubic polynomial characteristic equation of ${\mathcal{L}}^{\text{pl}}$, we will study its roots only for wave vectors that are parallel or perpendicular to P ₀, i.e. those for which either q_⊥ = 0 or q_∥ = 0 respectively. For both of these cases, we further assume that the limit γ → ∞ has been taken. In other words, we will study the roots of the two polynomial equations

$$\begin{equation}\mathrm{det}\left({\sigma }_{{\Vert}}-{\mathcal{L}}^{\text{pl}}\left({q}_{{\Vert}},0,\gamma \to \infty \right)\right)=0,\end{equation} \tag{ A.31 }$$

$$\begin{equation}\mathrm{det}\left({\sigma }_{\perp }-{\mathcal{L}}^{\text{pl}}\left(0,{q}_{\perp },\gamma \to \infty \right)\right)=0,\end{equation} \tag{ A.32 }$$

and aim to deduce the conditions under which σ_∥ ≡ σ_∥(q_∥) and σ_⊥ ≡ σ_⊥(q_⊥) have positive real parts.

Starting with (A.31), it is straightforward to show that σ_∥ solves

$$\begin{equation}\left({\sigma }_{{\Vert}}-{{\Gamma}}_{\perp }-i\lambda {P}_{0}{q}_{{\Vert}}\right){g}_{{\Vert}}\left({\sigma }_{{\Vert}}\right)=0,\end{equation} \tag{ A.33 }$$

where the polynomial g_∥ is given by

$$\begin{equation}{g}_{{\Vert}}\left(\sigma \right)={\sigma }^{2}-\left({{\Gamma}}_{{\Vert}}+i\lambda {P}_{0}{q}_{{\Vert}}\right)\sigma +w{w}_{1}{q}_{{\Vert}}^{2}+iw{P}_{0}{q}_{{\Vert}}.\end{equation} \tag{ A.34 }$$

In order to determine the number of roots of g_∥ in the halfplane {Re σ > 0} we apply the argument principle (for an introduction to this method, see e.g. [56] chapter 5, theorem 5.1.4, or [57] chapter 5). It states that the number Z_R of zeros of g_∥ inside the semi-circle contour C_R = I_R ∪ A_R , where

$$\begin{equation}{I}_{R}=\left[-iR,iR\right],\end{equation} \tag{ A.35 }$$

$$\begin{equation}{A}_{R}=\left\{R\enspace {\text{e}}^{\text{i}}\theta :\theta \in \left[-\pi /2,\pi /2\right]\right\},\end{equation} \tag{ A.36 }$$

is given by the change in the argument of g_∥ as we trace C_R counterclockwise, i.e.

$$\begin{equation}{Z}_{R}={\frac{1}{2\pi i}\oint }_{{C}_{R}}\mathrm{d}\sigma \enspace \frac{{g}_{{\Vert}}^{\prime }\left(\sigma \right)}{{g}_{{\Vert}}\left(\sigma \right)}=\frac{{\Delta}{\mathrm{arg}}_{{C}_{R}}\left({g}_{{\Vert}}\right)}{2\pi i},\end{equation} \tag{ A.37 }$$

which is illustrated in figure A1. Thus, by determining Z_R in the limit R → ∞, we may in particular deduce the number of roots of g_∥ with positive real part. Moreover, the choice of contour is quite convenient when dealing with polynomials such as (A.34), since we know that on the arc A_R , g_∥(R e^iθ ) ∼ R² e^2iθ + O(R). In other words,

$$\begin{align}\hfill {\Delta}{\mathrm{arg}}_{{A}_{R}}\left({g}_{{\Vert}}\right)& =\mathrm{log}\enspace {g}_{{\Vert}}\left(R\enspace {\text{e}}^{\text{i}\pi /2}\right)-\mathrm{log}\enspace {g}_{{\Vert}}\left(R\enspace {\text{e}}^{-\text{i}\pi /2}\right)\hfill \\ \hfill & =2\pi i+O\left({R}^{-1}\right).\hfill \end{align} \tag{ A.38 }$$

To compute Z_R in the limit R → ∞, we are therefore left with having to find the change in the argument of g_∥ along I_R .

Zoom In Zoom Out Reset image size

Stability clearly requires that Z_R → 2 as R → ∞, or combining (A.37) and (A.38),

$$\begin{equation}\underset{R\to \infty }{\mathrm{lim}}\enspace {\Delta}{\mathrm{arg}}_{{I}_{R}}\left({g}_{{\Vert}}\right)=2\pi i.\end{equation} \tag{ A.39 }$$

This occurs if and only if the image g_∥(I_R ) wraps around the origin once, as illustrated in figure A1, or equivalently that the winding number of g_∥(I_R ) about the origin is 1. To investigate when this occurs, we decompose g_∥(iy) where y ∈ [−R, R] into its real and imaginary parts, thus

$$\begin{equation}\mathrm{R}\mathrm{e}\enspace {g}_{{\Vert}}\left(iy\right)=-{y}^{2}+\lambda {P}_{0}{q}_{{\Vert}}y+w{w}_{1}{q}_{{\Vert}}^{2},\end{equation} \tag{ A.40 }$$

$$\begin{equation}\mathrm{I}\mathrm{m}\enspace {g}_{{\Vert}}\left(iy\right)=-{{\Gamma}}_{{\Vert}}y+w{P}_{0}{q}_{{\Vert}}.\end{equation} \tag{ A.41 }$$

Firstly, from (A.41) we see that we must require Γ_∥ > 0, or the winding number of g_∥(I_R ) could only be 0 or −1 (recall that we are tracing the line segment I_R = [−iR, iR] from iR to −iR since C_R is traced counterclockwise). This is ensured so long as ρ₀ > 1, which we assume in the following. Furthermore, from (A.40), it follows that the quadratic Re g_∥(iy) has two distinct real roots for all q_∥ ≠ 0, given by

$$\begin{equation}{y}_{{\pm}}=\frac{{q}_{{\Vert}}}{2}\left(\lambda {P}_{0}{\pm}\sqrt{{\left(\lambda {P}_{0}\right)}^{2}+4w{w}_{1}}\right).\end{equation} \tag{ A.42 }$$

Thus, we only need to require that Im g_∥(iy₊) < 0 and Im g_∥(iy₋) > 0. One may show by standard means that this occurs if and only if

$$\begin{equation}{\rho }_{0}{ >}1+\frac{1}{2}\frac{w}{\lambda +2{w}_{1}}.\end{equation} \tag{ A.43 }$$

In conclusion therefore, it follows that all the roots of the characteristic equation of ${\mathcal{L}}^{\text{pl}}$ are positive only when (A.43) is satisfied.

Proceeding with (A.32), i.e. the characteristic equation of ${\mathcal{L}}^{\text{pl}}$ for wave vectors that are perpendicular to P ₀, we apply the argument principle once more. In this case, the roots σ_⊥ solve the cubic polynomial equation

$$\begin{equation}{g}_{\perp }\left({\sigma }_{\perp }\right)\equiv {\sigma }_{\perp }^{3}-{c}_{2}{\sigma }_{\perp }^{2}+{c}_{1}{\sigma }_{\perp }-{c}_{0}=0,\end{equation} \tag{ A.44 }$$

where the coefficients c₂, c₁ and c₀ are given by

$$\begin{equation}{c}_{2}={{\Gamma}}_{{\Vert}}+{{\Gamma}}_{\perp },\end{equation} \tag{ A.45 }$$

$$\begin{equation}{c}_{1}={{\Gamma}}_{{\Vert}}{{\Gamma}}_{\perp }+{q}_{\perp }^{2}\left(w{w}_{1}-{\kappa }^{2}{P}_{0}^{2}\right),\end{equation} \tag{ A.46 }$$

$$\begin{equation}{c}_{0}=w{q}_{\perp }^{2}\left({w}_{1}{{\Gamma}}_{{\Vert}}-\kappa {P}_{0}^{2}\right).\end{equation} \tag{ A.47 }$$

By considering the image of C_R under g_⊥, we see that along the semi-circle arc A_R , the change in the argument of g_⊥ is 3πi + O(R⁻¹). Thus, in this case we must require that the winding number of g_⊥(I_R ) is $\frac{3}{2}$ in order to have Z_R → 3. This occurs if and only if g_⊥(I_R ) wraps around the origin in the way illustrated in figure A1. To uncover the conditions under which this occurs, we again decompose g_⊥(iy) into its real and imaginary parts:

$$\begin{equation}\mathrm{R}\mathrm{e}\enspace {g}_{\perp }\left(iy\right)={c}_{2}{y}^{2}-{c}_{0},\end{equation} \tag{ A.48 }$$

$$\begin{equation}\mathrm{I}\mathrm{m}\enspace {g}_{\perp }\left(iy\right)=-{y}^{3}+{c}_{1}y.\end{equation} \tag{ A.49 }$$

Assuming that ρ₀ > 1, we have c₂ > 0 which is necessary in order for the winding number of g_⊥(I_R ) to be positive. Furthermore, from (A.48) we see that we must require c₀ > 0 so that the quadratic Re g_⊥(iy) has two distinct real roots. This holds if and only if

$$\begin{equation}\kappa {< }2{w}_{1}.\end{equation} \tag{ A.50 }$$

Similarly, from (A.49) we deduce that we must have c₁ > 0, or equivalently

$$\begin{equation}{\kappa }^{2}{< }2+\frac{w{w}_{1}}{{\rho }_{0}-1},\end{equation} \tag{ A.51 }$$

so that the cubic Im g_⊥(iy) has three distinct real roots. Finally, requiring that $\mathrm{R}\mathrm{e}\enspace {g}_{\perp }\left({\pm}i\sqrt{{c}_{1}}\right){ >}0$ one straightforwardly verifies that we are indeed guaranteed that Z_R → 3. This last condition is equivalent to having

$$\begin{equation}{c}_{2}{c}_{1}-{c}_{0}{ >}0.\end{equation} \tag{ A.52 }$$

One may show that (A.52) holds if and only if either

$$\begin{equation}{\kappa }^{2}{< }3+\frac{w{w}_{1}}{2\left({\rho }_{0}-1\right)}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad 2{P}_{0}^{2}\left(2-{\kappa }^{2}\right)+w\kappa { >}0\end{equation} \tag{ A.53 }$$

or

$$\begin{equation}{\left(2{P}_{0}^{2}\left(3-{\kappa }^{2}\right)+w{w}_{1}\right)}^{2}{< }8{P}_{0}^{2}\left(2{P}_{0}^{2}\left(2-{\kappa }^{2}\right)+w\kappa \right).\end{equation} \tag{ A.54 }$$

For all simulations we perform, conditions (A.51) and (A.52) are satisfied. Because of this, we will only be concerned the stability requirement in (A.50). It is also worth highlighting once more that, despite the rather involved analysis undertaken here, we have not solved the full cubic characteristic equation of ${\mathcal{L}}^{\text{pl}}$ and thus have not fully identified all necessary and sufficient conditions for linear stability.

Lastly, we also investigate the structure of the correlators of the linearised theory about the PL state. In analogy with the isotropic calculation, the matrix ${\mathcal{C}}^{\text{pl}}$ in (52) solves an algebraic Riccati equation as in (A.18), although in this case all three modes $\left({\hat{\rho }}_{1},{\hat{P}}_{{\Vert}},{\hat{P}}_{\perp }\right)$ remain coupled for general q . Finding its solution can be achieved by solving the linear system

$$\begin{equation}{\mathcal{R}}^{\text{pl}}{\boldsymbol{C}}^{\text{pl}}=2\boldsymbol{D},\end{equation} \tag{ A.55 }$$

where ${\boldsymbol{C}}^{\text{pl}}={\left({\mathcal{C}}_{11}^{\text{pl}},{\mathcal{C}}_{12}^{\text{pl}},\dots ,{\mathcal{C}}_{33}^{\text{pl}}\right)}^{\mathrm{T}}$, $\boldsymbol{D}={\left({\mathcal{D}}_{11},{\mathcal{D}}_{12},\dots ,{\mathcal{D}}_{33}\right)}^{\mathrm{T}}$. We avoid explicitly writing out the nine-by-nine matrix ${\mathcal{R}}^{\text{pl}}$ here for sake of clarity of presentation, although it may be found fairly straightforwardly from the Riccati equation. Analytical inversion of (A.55) may be done using standard computer algebra systems that perform symbolic computations. Due to the algebraic complexity of the resulting expressions, we choose to only state the result in certain limiting cases. More specifically, we look at the asymptotic behaviour of the components ${\mathcal{C}}_{ij}^{\text{pl}}$ in the limit q → ∞ at fixed q_∥/q_⊥, as well as when γ → ∞ is taken first. From this, we deduce in appendix C the dependence of the EPR on the ultraviolet cutoff Λ quoted in the main text. For the six independent components of ${\mathcal{C}}^{\text{pl}}$ we summarize the results in table A1.

Table A1. Leading order asymptotic expressions for the components ${\mathcal{C}}_{ij}^{\text{pl}}$ of the equal-time correlation matrix ${\mathcal{C}}^{\text{pl}}$ in the limit where q → ∞, as well as when γ → ∞ is taken first.

${\mathcal{C}}_{ij}^{\text{pl}}$	γ → ∞, q → ∞	q → ∞
${\mathcal{C}}_{11}^{\text{pl}}$	$w/\left(\mathcal{V}{w}_{1}{q}^{2}\right)$	$1/\left(\mathcal{V}{\nu }_{\rho }{q}^{2}\right)$
${\mathcal{C}}_{12}^{\text{pl}}$	${P}_{0}w\left({\mathrm{cos}}^{2}\enspace \theta +\left(1+\kappa {w}_{1}\right){\mathrm{sin}}^{2}\enspace \theta \right)/\left(\mathcal{V}{w}_{1}{q}^{4}\right)$	$i{w}_{1}\enspace \mathrm{cos}\enspace \theta /\left(\mathcal{V}{\nu }_{\rho }{q}^{3}\right)$
${\mathcal{C}}_{13}^{\text{pl}}$	$-{P}_{0}\kappa w\enspace \mathrm{cos}\enspace \theta \enspace \mathrm{sin}\enspace \theta /\left(\mathcal{V}{q}^{4}\right)$	$i{w}_{1}\enspace \mathrm{sin}\enspace \theta /\left(\mathcal{V}{\nu }_{\rho }{q}^{3}\right)$
${\mathcal{C}}_{22}^{\text{pl}}$	$1/\left(\mathcal{V}{q}^{2}\right)$	$1/\left(\mathcal{V}{q}^{2}\right)$
${\mathcal{C}}_{23}^{\text{pl}}$	$-i{P}_{0}\kappa \enspace \mathrm{sin}\enspace \theta /\left(\mathcal{V}{q}^{3}\right)$	$-i{P}_{0}\kappa \enspace \mathrm{sin}\enspace \theta /\left(\mathcal{V}{q}^{3}\right)$
${\mathcal{C}}_{33}^{\text{pl}}$	$1/\left(\mathcal{V}{q}^{2}\right)$	$1/\left(\mathcal{V}{q}^{2}\right)$

Starting from the definitions (17) for the EPR $\dot {\mathcal{S}}$ and (13) for the path transition probability density $\mathcal{P}$ of the HVM, one finds that we may equivalently express $\dot {\mathcal{S}}$ in terms of the $\mathcal{T}$-antisymmetric part of the Freidlin–Wentzell action $\mathcal{A}$, i.e.

$$\begin{equation}\dot {\mathcal{S}}=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{{\leftarrow}{\mathcal{A}}-\mathcal{A}}{2D\tau },\end{equation} \tag{ C.1 }$$

where ${\leftarrow}{\mathcal{A}}=\mathcal{A}{\circ}\mathcal{T}$ as before. By applying $\mathcal{T}$ to $\mathcal{A}$ as given in (12), we find that the action ${\leftarrow}{\mathcal{A}}$ for the time-reversed ensemble may be expressed more explicitly by

$$\begin{equation}{\leftarrow}{\mathcal{A}}=\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace {\left\vert {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}-\frac{\delta {F}^{S}}{\delta \boldsymbol{P}}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right\vert }^{2}\quad \text{if}\enspace {\partial }_{t}\rho +w\nabla \cdot \boldsymbol{P}=0,\end{equation} \tag{ C.2 }$$

and ${\leftarrow}{\mathcal{A}}=\infty$ otherwise. In particular, from (C.1) one straightforwardly deduces that

$$\begin{equation}\dot {\mathcal{S}}=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{-1}{2D\tau }{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left(\odot {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right)\cdot \frac{\delta {F}^{S}}{\delta \boldsymbol{P}}.\end{equation} \tag{ C.3 }$$

In (C.3)—and throughout this appendix—we indicate by ⊙ that the product with ∂_t P should be interpreted in the Stratonovich sense, i.e. employ mid-point discretisation in time.

To transform this into the expression given in (18), we simply have to observe that some terms in (C.3) are integrable. Specifically, we find that we may write

$$\begin{equation}\dot {\mathcal{S}}=\underset{\tau \to \infty }{\mathrm{lim}}\left[\frac{-{\Delta}\mathcal{I}}{2D\tau }+\frac{1}{2D\tau }{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left(\rho \boldsymbol{P}\odot {\partial }_{t}\boldsymbol{P}-\left(\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right)\cdot \frac{\delta {F}^{S}}{\delta \boldsymbol{P}}\right)\right],\end{equation} \tag{ C.4 }$$

where we have defined $\mathcal{I}$ to be the functional

$$\begin{equation}\mathcal{I}\left[\boldsymbol{P}\right]={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left(\frac{1}{2}\vert \boldsymbol{P}{\vert }^{2}+\frac{1}{2}{\left({\nabla }_{\alpha }{P}_{\beta }\right)}^{2}+\frac{1}{4}\vert \boldsymbol{P}{\vert }^{4}\right),\end{equation} \tag{ C.5 }$$

i.e. as the part of F^S which does not explicitly depend on the density ρ. Since we assume that the moments of P and its higher order spatial derivatives are finite in steady-state, it follows that ${\Delta}\mathcal{I}/\tau \to 0$ as τ → ∞ so that the first term in (C.4) may safely be ignored. In a similar vein, an integration by parts allows us to write the integral over the first term appearing in the integrand in (C.4) as

$$\begin{equation}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \rho \odot {\partial }_{t}\vert \boldsymbol{P}{\vert }^{2}={\left.{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \rho \vert \boldsymbol{P}{\vert }^{2}\right\vert }_{-\tau }^{\tau }+w{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left(\nabla \cdot \boldsymbol{P}\right)\vert \boldsymbol{P}{\vert }^{2},\end{equation} \tag{ C.6 }$$

where we have used the continuity equation ∂_t ρ = −w∇ ⋅ P . Again, after dividing by 4Dτ, the first term on the right-hand side of (C.6) goes away in the limit τ → ∞ since it grows sublinearly in τ. Thus, assuming that we may replace temporal averages by averages over noise realizations, we finally arrive at the expression

$$\begin{equation}\dot {\mathcal{S}}={D}^{-1}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \frac{w}{2}\vert \boldsymbol{P}{\vert }^{2}\left(\nabla \cdot \boldsymbol{P}\right)-\left(\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right)\cdot \frac{\delta {F}^{S}}{\delta \boldsymbol{P}}\right\rangle \end{equation} \tag{ C.7 }$$

reported in the main text.

Next, we investigate the expansion of (C.7) about the constant homogeneous ground-states. To do this, we substitute the expansions (25) and (26) into (C.7) and collect terms at equal order in D. It is fairly straightforward to check that there are no contributions to $\dot {\mathcal{S}}$ at orders D⁻¹ or D^−1/2. At order D⁰ we find after some fairly tedious algebra that

$$\begin{align}\hfill {\dot {\mathcal{S}}}_{0}& ={\boldsymbol{P}}_{0}\cdot {\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle w\left(\nabla \cdot {\boldsymbol{P}}_{1}\right){\boldsymbol{P}}_{1}+{\rho }_{1}\left(\left(2{w}_{1}-\kappa \right)\nabla \left({\boldsymbol{P}}_{0}\cdot {\boldsymbol{P}}_{1}\right)+\lambda \left({\boldsymbol{P}}_{0}\cdot \nabla \right){\boldsymbol{P}}_{1}\right)\right.\hfill \\ \hfill & \quad +\left.2\kappa \left(\nabla \cdot {\boldsymbol{P}}_{1}\right)\left({\nabla }^{2}{\boldsymbol{P}}_{1}-\vert {\boldsymbol{P}}_{0}{\vert }^{2}{\boldsymbol{P}}_{1}\right)\right\rangle .\hfill \end{align} \tag{ C.8 }$$

In order to arrive at this expression we have only assumed that P ₀ satisfies the zeroth order equation (A.2), and used (A.1) repeatedly. If we further assume that P ₀ is a polarised solution with | P ₀| > 0, we may write P ₀ = (P₀, 0) without loss of generality, from which (35) follows immediately after integrating out total derivatives.

After transforming (C.8) to Fourier space, we obtain (49) as stated in the main text. Using our results from appendix A, we may now investigate the scaling of this expression with Λ → ∞. To determine this, we note that the scaling of the sum in this limit is determined by the corresponding integral in q-space, i.e.

$$\begin{equation}{\dot {\mathcal{S}}}_{0}/\mathcal{V}=\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\mathrm{Tr}\left({\dot {\sigma }}^{\text{pl}}\enspace {C}^{\text{pl}}\right)\sim \frac{\mathcal{V}}{{\left(2\pi \right)}^{2}}{\int }_{0}^{2\pi }\mathrm{d}\theta {\int }_{0}^{{\Lambda}}\mathrm{d}q\enspace q\mathrm{Tr}\left({\dot {\sigma }}^{\text{pl}}\enspace {C}^{\text{pl}}\right),\end{equation} \tag{ C.9 }$$

where ${\dot {\sigma }}^{\text{pl}}$ is given by (50). In particular, direct substitution from table A1 and subsequently performing the integrals over q and θ yields

$$\begin{equation}{\dot {\mathcal{S}}}_{0}/\mathcal{V}\sim \frac{{P}_{0}^{2}{\kappa }^{2}}{4\pi }{{\Lambda}}^{2}.\end{equation} \tag{ C.10 }$$

In analogy with the above calculation, we may compute the EPRs ${\dot {\mathcal{S}}}^{{\pm}}$ of the DFM from the ${\mathcal{T}}^{{\pm}}$-antisymmetric part of the Freidlin–Wentzell action ${\mathcal{A}}_{\mathrm{D}\mathrm{F}}$, specifically

$$\begin{equation}{\dot {\mathcal{S}}}^{{\pm}}=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{{{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{{\pm}}-{\mathcal{A}}_{\mathrm{D}\mathrm{F}}}{2D\tau },\end{equation} \tag{ C.11 }$$

where we have defined ${{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{{\pm}}={\mathcal{A}}_{\mathrm{D}\mathrm{F}}{\circ}{\mathcal{T}}^{{\pm}}$. Writing out the actions for the time-reversed ensembles explicitly, we find that

$$\begin{equation}{{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{+}=\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left[\gamma {\left\vert {\nabla }^{-1}\left({\partial }_{t}\rho -\nabla \cdot {\boldsymbol{J}}_{\mathrm{d}}\right)\right\vert }^{2}+{\left\vert {\partial }_{t}\boldsymbol{P}-\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}-\frac{\delta F}{\delta \boldsymbol{P}}\right\vert }^{2}\right],\end{equation} \tag{ C.12 }$$

and

$$\begin{align}\hfill {{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{-}& =\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left[\gamma {\left\vert {\nabla }^{-1}\left({\partial }_{t}\rho -\nabla \cdot {\boldsymbol{J}}_{\mathrm{d}}^{S}+\nabla \cdot {\boldsymbol{J}}_{\mathrm{d}}^{A}\right)\right\vert }^{2}\right.\left.+{\left\vert {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}-\frac{\delta {F}^{S}}{\delta \boldsymbol{P}}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right\vert }^{2}\right].\hfill \end{align} \tag{ C.13 }$$

Thus, proceeding as above we straightforwardly find that

$$\begin{align}\hfill {\dot {\mathcal{S}}}^{+}& =\underset{\tau \to \infty }{\mathrm{lim}}\left[\frac{-{\Delta}F}{2D\tau }\right.\left.+\frac{1}{2D\tau }{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left(\gamma w\left(\nabla \cdot \boldsymbol{P}\right){\nabla }^{-2}\odot {\partial }_{t}\rho -\lambda \left[\left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}\right]\odot {\partial }_{t}\boldsymbol{P}\right)\right]\hfill \end{align} \tag{ C.14 }$$

and

$$\begin{align}\hfill {\dot {\mathcal{S}}}^{-}& =\underset{\tau \to \infty }{\mathrm{lim}}\left[\frac{-{\Delta}{F}^{S}}{2D\tau }\right.\left.+\frac{1}{2D\tau }{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left({\boldsymbol{J}}_{\mathrm{d}}^{A}\cdot \nabla {\mu }^{S}-\left(\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right)\cdot \frac{\delta {F}^{S}}{\delta \boldsymbol{P}}\right)\right].\hfill \end{align} \tag{ C.15 }$$

The expression for ${\dot {\mathcal{S}}}^{-}$ immediately reduces to (74) provided in the main text upon replacing the temporal average by an average over noise realizations and noting again that ΔF^S /τ → 0 as τ → ∞. On the other hand, (C.14) requires a bit more massaging. Firstly, we observe that the equal-time expectation

$$\begin{equation}\left\langle \left(\nabla \cdot \boldsymbol{P}\right){\nabla }^{-2}\nabla \odot \boldsymbol{\xi }\right\rangle =0,\end{equation} \tag{ C.16 }$$

since here the Stratonovich product coincides with the corresponding Ito product (i.e. there is no spurious drift term). It follows that

$$\begin{equation}\left\langle \left(\nabla \cdot \boldsymbol{P}\right){\nabla }^{-2}\odot {\partial }_{t}\rho \right\rangle =-w\left\langle \left(\nabla \cdot \boldsymbol{P}\right){\nabla }^{-2}\nabla \cdot \boldsymbol{P}\right\rangle +{\gamma }^{-1}\left\langle \left(\nabla \cdot \boldsymbol{P}\right)\mu \right\rangle .\end{equation} \tag{ C.17 }$$

To treat the product with ⊙∂_t P in (C.14) we must explicitly compute the spurious drift. We will show that, in fact, the spurious drift is a total derivative and therefore does not contribute to the EPR ${\dot {\mathcal{S}}}^{+}$. To do this, we consider a finite discretisation of the process in Fourier space with | q | ⩽ Λ, which is consistent with our numerical scheme. By applying standard stochastic calculus, we then obtain

$$\begin{align}\hfill {\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \left[\left(\boldsymbol{P}\cdot \nabla \right)\boldsymbol{P}\right]\odot \boldsymbol{\eta }\right\rangle & =\mathcal{V}\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\sum _{\vert \boldsymbol{k}\vert {\leqslant}{\Lambda}}i{k}_{\beta }\left\langle {P}_{\beta }\left(\boldsymbol{q}-\boldsymbol{k}\right){P}_{\alpha }\left(\boldsymbol{k}\right)\odot {\eta }_{\alpha }\left(-\boldsymbol{q}\right)\right\rangle \hfill \\ \hfill & =D\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\sum _{\vert \boldsymbol{k}\vert {\leqslant}{\Lambda}}i{k}_{\beta }\left\langle \left({\delta }_{\boldsymbol{k},\boldsymbol{0}}{P}_{\beta }\left(\boldsymbol{k}\right)+2{\delta }_{\boldsymbol{q},\boldsymbol{k}}{P}_{\beta }\left(\boldsymbol{q}-\boldsymbol{k}\right)\right)\right\rangle \hfill \\ \hfill & =0.\hfill \end{align} \tag{ C.18 }$$

Here, the second equality follows from the transformation rule between Ito and Stratonovich processes [58], i.e.

$$\begin{equation}\left\langle h\left(\boldsymbol{P}\left({\boldsymbol{q}}_{1}\right),\dots ,\boldsymbol{P}\left({\boldsymbol{q}}_{n}\right)\right)\odot {\eta }_{\alpha }\left(\boldsymbol{k}\right)\right\rangle =D\sum _{m=1}^{n}\left\langle \frac{\partial h}{\partial {P}_{\beta }\left({\boldsymbol{q}}_{m}\right)}\right\rangle {\delta }_{{\boldsymbol{q}}_{m},-\boldsymbol{k}}{\delta }_{\alpha \beta },\end{equation} \tag{ C.19 }$$

and the fact that δ_αα = 2. To see why the final equality holds, note also that

$$\begin{equation}\sum _{\vert \boldsymbol{k}\vert {\leqslant}{\Lambda}}{k}_{\alpha }=0\end{equation} \tag{ C.20 }$$

since the sum is finite. From this, and taking τ → ∞ in (C.14), the result (73) reported in the main text follows immediately.

For the DFM, we may additionally calculate the EPRs ${\dot {\mathcal{S}}}_{J}^{{\pm}}$ at the level of the fluctuating density current J . At this level, the actions ${{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{J,{\pm}}$ for the time-reversed ensembles may be expressed as

$$\begin{equation}{{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{J,+}=\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left[\gamma {\left\vert \boldsymbol{J}+{\boldsymbol{J}}_{\mathrm{d}}\right\vert }^{2}+{\left\vert {\partial }_{t}\boldsymbol{P}-\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}-\frac{\delta F}{\delta \boldsymbol{P}}\right\vert }^{2}\right]\end{equation} \tag{ C.21 }$$

and

$$\begin{equation}{{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{J,-}=\frac{1}{4}{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left[\gamma {\left\vert \boldsymbol{J}+{\boldsymbol{J}}_{\mathrm{d}}^{S}-{\boldsymbol{J}}_{\mathrm{d}}^{A}\right\vert }^{2}+{\left\vert {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}-\frac{\delta {F}^{S}}{\delta \boldsymbol{P}}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right\vert }^{2}\right]\end{equation} \tag{ C.22 }$$

if ∂_t ρ + ∇ ⋅ J = 0 and ${{\leftarrow}{\mathcal{A}}}_{\mathrm{D}\mathrm{F}}^{J,{\pm}}=\infty$ otherwise. By direct substitution we then find that

$$\begin{equation}{\dot {\mathcal{S}}}_{J}^{+}=\frac{1}{2D\tau }{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left(\gamma {\boldsymbol{J}}_{\mathrm{d}}\odot \boldsymbol{J}-\left(\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta F}{\delta \boldsymbol{P}}\right)\odot {\partial }_{t}\boldsymbol{P}\right),\end{equation} \tag{ C.23 }$$

in addition to

$$\begin{align}\hfill {\dot {\mathcal{S}}}_{J}^{-}& =\frac{1}{2D\tau }{\int }_{-\tau }^{\tau }\mathrm{d}t{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\left(\gamma \left(\odot \boldsymbol{J}-{\boldsymbol{J}}_{\mathrm{d}}^{A}\right)\cdot {\boldsymbol{J}}_{\mathrm{d}}^{S}\right.\left.-\left(\odot {\partial }_{t}\boldsymbol{P}+\lambda \boldsymbol{P}\cdot \nabla \boldsymbol{P}+\frac{\delta {F}^{A}}{\delta \boldsymbol{P}}\right)\cdot \frac{\delta {F}^{S}}{\delta \boldsymbol{P}}\right).\hfill \end{align} \tag{ C.24 }$$

Now, it is fairly easy to see that

$$\begin{align}\hfill {\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {\boldsymbol{J}}_{\mathrm{d}}\odot \boldsymbol{J}\right\rangle & ={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle w\boldsymbol{P}\odot \boldsymbol{J}-{\gamma }^{-1}\odot \boldsymbol{J}\cdot \nabla \mu \right\rangle \hfill \\ \hfill & ={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {w}^{2}\vert \boldsymbol{P}{\vert }^{2}-{\gamma }^{-1}w\boldsymbol{P}\cdot \nabla \mu -{\gamma }^{-1}\frac{\delta F}{\delta \rho }\odot {\partial }_{t}\rho \right\rangle ,\hfill \end{align} \tag{ C.25 }$$

where the second equality follows from an integration by parts and the fact that ⟨ P ⊙ J ⟩ = ⟨ P ⋅ J _d⟩. Substituting this back into (C.23) gives the desired result for ${\dot {\mathcal{S}}}_{J}^{+}$, stated in (78). Similarly, we have that

$$\begin{equation}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \left(\odot \boldsymbol{J}-{\boldsymbol{J}}_{\mathrm{d}}^{A}\right)\cdot {\boldsymbol{J}}_{\mathrm{d}}^{S}\right\rangle =-{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {\gamma }^{-1}\frac{\delta {F}^{S}}{\delta \rho }\odot {\partial }_{t}\rho +{\boldsymbol{J}}_{\mathrm{d}}^{A}\cdot {\boldsymbol{J}}_{\mathrm{d}}^{S}\right\rangle ,\end{equation} \tag{ C.26 }$$

from which the fact that ${\dot {\mathcal{S}}}_{J}^{-}={\dot {\mathcal{S}}}^{-}$ follows upon substitution back into (C.24).

Again we may linearize the expressions (73) and (74) about the homogeneous isotropic and PL states by substituting in an expansion of the form (25) and (26). Treating this as above for the HVM, we find that

$$\begin{align}\hfill {\dot {\mathcal{S}}}_{0}^{+}& ={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle \gamma {w}^{2}\vert K{\boldsymbol{P}}_{1}{\vert }^{2}+w\left(\nabla \cdot {\boldsymbol{P}}_{1}\right)\left({a}_{\rho }{\rho }_{1}-{\nu }_{\rho }{\nabla }^{2}{\rho }_{1}\right)-{P}_{0}w{P}_{{\Vert}}\left({\partial }_{\perp }{P}_{\perp }\right)\right.\hfill \\ \hfill & \quad \left.-w{w}_{1}{\left(\nabla \cdot {\boldsymbol{P}}_{1}\right)}^{2}+{P}_{0}^{2}{\lambda }^{2}\vert {\partial }_{{\Vert}}{\boldsymbol{P}}_{1}{\vert }^{2}-{P}_{0}^{2}\lambda {\rho }_{1}\left({\partial }_{{\Vert}}{P}_{{\Vert}}\right)\right.\left.+{P}_{0}\lambda {w}_{1}\left({\partial }_{{\Vert}}{\rho }_{1}\right)\left(\nabla \cdot {\boldsymbol{P}}_{1}\right)\right\rangle ,\hfill \end{align} \tag{ C.27 }$$

$$\begin{align}\hfill {\dot {\mathcal{S}}}_{0}^{-}& ={\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{x}\enspace \left\langle {P}_{0}{P}_{{\Vert}}{\partial }_{\perp }\left(w{P}_{\perp }+{\gamma }^{-1}{w}_{1}{\nabla }^{2}{P}_{\perp }\right)\right.\left.-\left({a}_{\rho }{\rho }_{1}-{\nu }_{\rho }{\nabla }^{2}{\rho }_{1}\right)\nabla \cdot \left(w{\boldsymbol{P}}_{1}+{\gamma }^{-1}{w}_{1}{\nabla }^{2}{\boldsymbol{P}}_{1}\right)\right.\hfill \\ \hfill & \quad \left.+{P}_{0}^{2}\left(2{w}_{1}-\kappa +\lambda \right){\rho }_{1}{\partial }_{{\Vert}}{P}_{{\Vert}}-2{P}_{0}\kappa \left({\partial }_{\perp }{P}_{\perp }\right)\left({P}_{0}^{2}{P}_{{\Vert}}-{\nabla }^{2}{P}_{{\Vert}}\right)\right\rangle ,\hfill \end{align} \tag{ C.28 }$$

where both expressions hold for general P₀ ⩾ 0. We may equivalently express (C.27) and (C.28) in Fourier space, and for the constant homogeneous ground-states we obtain expressions analogous to those which we encountered for the HVM (49). Specifically, we find that for the isotropic and PL states,

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{{\pm}}/\mathcal{V}=\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\mathrm{Tr}\left({\dot {\sigma }}^{{\pm},\mathrm{i}\mathrm{s}\mathrm{o}}{C}^{\text{iso}}\right),\end{equation} \tag{ C.29 }$$

and

$$\begin{equation}{\dot {\mathcal{S}}}_{0}^{{\pm}}/\mathcal{V}=\sum _{\vert \boldsymbol{q}\vert {\leqslant}{\Lambda}}\mathrm{Tr}\left({\dot {\sigma }}^{{\pm},\mathrm{p}\mathrm{l}}{C}^{\text{pl}}\right),\end{equation} \tag{ C.30 }$$

respectively. Furthermore, we may choose to write the sum such that ${\dot {\sigma }}^{{\pm},\mathrm{i}\mathrm{s}\mathrm{o}}$ and ${\dot {\sigma }}^{{\pm},\mathrm{p}\mathrm{l}}$ are Hermitian. Taking P₀ = 0 in (C.27) and (C.28) and transforming to Fourier space we find that ${\dot {\sigma }}^{{\pm},\mathrm{i}\mathrm{s}\mathrm{o}}$ are given by (93) and (94) as advertised. For σ^±,pl we list the six independent components of each matrix in table C1. Finally, from (C.30) in addition to tables A1 and C1, we straightforwardly deduce that in the PL phase, ${\dot {\mathcal{S}}}_{0}^{+}/\mathcal{V}\sim {P}_{0}^{2}{\lambda }^{2}{{\Lambda}}^{2}/\left(4\pi \right)$ and ${\dot {\mathcal{S}}}_{0}^{-}/\mathcal{V}\sim {w}_{1}^{2}{{\Lambda}}^{4}/\left(8\pi \gamma \right)$, while the exact results in the isotropic phase are presented in the main text.

Table C1. Independent components of the Hermitian bilinear EPR coupling matrices ${\dot {\sigma }}^{{\pm},\mathrm{p}\mathrm{l}}$.

(i, j)	${\dot {\sigma }}_{ij}^{+,\mathrm{p}\mathrm{l}}$	${\dot {\sigma }}_{ij}^{-,\mathrm{p}\mathrm{l}}$
(1, 1)	0	0
(1, 2)	$\frac{i}{2}\left(w{{\Gamma}}_{\rho }-{P}_{0}^{2}\lambda -i{P}_{0}\lambda {w}_{1}{q}_{{\Vert}}\right){q}_{{\Vert}}$	$\frac{i}{2}\left({w}_{1}{{\Gamma}}_{{\Vert}}-\tilde {w}{{\Gamma}}_{\rho }+{P}_{0}^{2}\lambda \right){q}_{{\Vert}}$
(1, 3)	$\frac{i}{2}\left(w{{\Gamma}}_{\rho }-i{P}_{0}\lambda {w}_{1}{q}_{{\Vert}}\right){q}_{\perp }$	$\frac{i}{2}\left({w}_{1}{{\Gamma}}_{\perp }-\tilde {w}{{\Gamma}}_{\rho }+{P}_{0}^{2}\kappa \right){q}_{\perp }$
(2, 2)	$\left(\gamma w\tilde {w}+{P}_{0}^{2}{\lambda }^{2}{q}^{2}\right){q}_{{\Vert}}^{2}/{q}^{2}$	0
(2, 3)	$-\frac{i}{2}w\left({P}_{0}{q}^{2}+2i\gamma \tilde {w}{q}_{{\Vert}}\right){q}_{\perp }/{q}^{2}$	$\frac{i}{2}{P}_{0}\left(\tilde {w}-\kappa \left({{\Gamma}}_{{\Vert}}+{{\Gamma}}_{\perp }\right)\right){q}_{\perp }$
(3, 3)	$\left(\gamma w\tilde {w}+{P}_{0}^{2}{\lambda }^{2}{q}^{2}{\left(\frac{{q}_{{\Vert}}}{{q}_{\perp }}\right)}^{2}\right){q}_{\perp }^{2}/{q}^{2}$	0