# Active XY model on a substrate: Density fluctuations and phase ordering

Astik Haldar<EMAIL_ADDRESS>Theory Division, Saha Institute of
Nuclear Physics, HBNI, 1/AF Bidhannagar, Calcutta 700064, West Bengal, India
Apurba Sarkar<EMAIL_ADDRESS>School of Mathematical &
Computational Sciences, Indian Association for the Cultivation of Science,
Kolkata-700032, West Bengal, India  Swarnajit Chatterjee
<EMAIL_ADDRESS>Center for Biophysics & Department for
Theoretical Physics, Saarland University, 66123 Saarbrücken, Germany  Abhik
Basu<EMAIL_ADDRESS><EMAIL_ADDRESS>Theory Division, Saha
Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Calcutta 700064, West
Bengal, India

###### Abstract

We explore the generic long wavelength properties of an active XY model on a
substrate, consisting of collection of nearly phase-ordered active XY spins in
contact with a diffusing, conserved species, as a representative system of
active spinners with a conservation law. The spins rotate actively in response
to the local density fluctuations and local phase differences, on a solid
substrate. We investigate this system by Monte-Carlo simulations of an agent-
based model, which we set up, complemented by the hydrodynamic theory for the
system. We demonstrate that this system can phase-synchronize without any
hydrodynamic interactions. Our combined numerical and analytical studies show
that this model, when stable, displays hitherto unstudied scaling behavior: As
a consequence of the interplay between the mobility, active rotation and
number conservation, such a system can be stable over a wide range of the
model parameters characterized by a novel correspondence between the phase and
density fluctuations. In different regions of the phase space where the phase-
ordered system is stable, it shows phase ordering which is generically either
logarithmically stronger than the conventional quasi long range order (QLRO)
found in its equilibrium limit, together with “miniscule number fluctuations”,
or logarithmically weaker than QLRO along with “giant number fluctuations”,
showing a novel one-to-one correspondence between phase ordering and density
fluctuations in the ordered states. Intriguingly, these scaling exponents are
found to depend explicitly on the model parameters. We further show that in
other parameter regimes there are no stable, ordered phases. Instead, two
distinct types of disordered states with short range phase-order are found,
characterized by the presence or absence of stable clusters of finite sizes.
In a surprising connection, the hydrodynamic theory for this model also
describes the fluctuations in a Kardar-Parisi-Zhang (KPZ) surface with a
conserved species on it, or an active fluid membrane with a finite tension,
without momentum conservation and a conserved species living on it. This
implies the existence of stable fluctuating surfaces that are only
logarithmically smoother or rougher than the Edward-Wilkinson surface at two
dimensions ($2d$) can exist, in contrast to the $2d$ pure KPZ-like “rough”
surfaces.

##  I Introduction

Recent studies on nonequilibrium systems reveal a fundamental phenomenological
difference with equilibrium systems that nonequilibrium steady states (NESS)
can display orders, which are prohibited in equilibrium systems due to the
basic laws of equilibrium statistical mechanics. A well-known example is a
two-dimensional ($2d$) flock of self-propelled particles (e.g., a flock of
birds, a school of fish and their artificial imitations). These systems can
show true orientational long range order Toner and Tu (1995), in contravention
with $2d$ equilibrium systems with short range interactions and continuous
symmetries (e.g., the equilibrium XY model) that can only show quasi-long
range order (QLRO) but never a long range order Mermin and Wagner (1966).
Nonequlibrium systems can also display instabilities and phase separations not
found in equilibrium systems. For instance, active nematic liquid crystals are
found to be intrinsically phase separated Mishra and Ramaswamy (2006). Another
example of such contrasting behavior between an equilibrium system and its
nonequilibrium counterpart is a $2d$ superfluid: an equilibrium $2d$
superfluid displays algebraic order or QLRO Chaikin and Lubensky (1995),
whereas its active analogue can show such an order only if the system is
anisotropic Altman  _et al._ (2015). Yet another example is that of a
fluctuating surface: An “Edward-Wilkinson” (EW) Barabási  _et al._ (1995)
surface, that follows an equilibrium dynamics, is “logarithmically rough” at
$2d$, i.e., the variance of the height fluctuations grows logarithmically with
the system size, in an exact correspondence with the QLRO of the $2d$ XY model
at low temperature. In contrast, a surface that satisfies the well-known
Kardar-Parisi-Zhang (KPZ) equation Kardar  _et al._ (1986); Barabási  _et al._
(1995), a paradigmatic nonequilibrium model, is “algebraically rough”, i.e.,
the height fluctuation variance grows as a power law in the system size at
$2d$. Exploration of the natures of order, universality and steady states are
a primary goal in the general studies of nonequilibrium systems.

Universal scaling in driven nonequilibrium systems has been a topic of many
theoretical studies. Studying the active or nonequilibrium versions of models
of equilibrium statistical mechanics has been a useful and rewarding approach
that has yielded unexpectedly rich nonequilibrium physics. A particularly
fascinating example of this kind is the driven or active $2d$ XY model. Wide
ranging natural systems of physical, chemical or biological origin are
described in terms of driven or active XY spins. This includes, e.g., active
superfluids, which though generalize equilibrium superfluidity in ultra low-
temperature condensed matter systems Chaikin and Lubensky (1995), can also
occur in very different systems like live bacterial suspensions Takatori and
Brady (2017); Guo  _et al._ (2018); López  _et al._ (2015), often marked by
resistanceless flows. Yet another class of systems where (active) XY-like
phenomenologies is expected to be important includes many natural systems of
diverse origin are described in terms of collections of interacting
oscillators Strogatz (2018) that have same internal symmetry as the XY spin.
There are many such examples, e.g., chemical oscillators Miyakawa and Isikawa
(2002); Toiya  _et al._ (2010, 2008); Fukuda  _et al._ (2003); Kiss  _et al._
(2002); Migliorini (2008); Strogatz  _et al._ (1989); Paoletti  _et al._
(2006), synthetic genetic oscillators Danino  _et al._ (2010), and
biologically relevant systems Jülicher and Prost (1997); Fürthauer  _et al._
(2013); Fürthauer and Ramaswamy (2013); Leoni and Liverpool (2014); Naganathan
_et al._ (2014); Jiang  _et al._ (2000); Oates  _et al._ (2012); Gregor  _et
al._ (2010); Laskar  _et al._ (2013). Another physically different but related
example is the KPZ equation for surface growth, which is intimately connected
to the $2d$ driven XY model and does not admit a “smooth” phase at $2d$
Barabási  _et al._ (1995).

The paramount issue in $2d$ driven XY model is the existence of a stable
orientationally ordered phase: What are the effects of conservation laws and
finite mobility on phase ordering, i.e., do mobility and conservation laws
together facilitate or hinder phase-ordering? Apart from its theoretical
motivation, it assumes importance as a prototype active model with diffusive
motion having rotational invariance that displays QLRO phase order and normal
number fluctuations in the equilibrium limit. Thus it offers an intuitive
understanding of the impact of drive and mobility on the equilibrium states of
$2d$ models with rotational invariance. For instance, isotropic active
superfluids without number conservation on substrates, e.g., driven fluids of
exciton polaritons, excitations of two-dimensional quantum wells in optical
cavities Carusotto and Ciuti (2013), have no algebraic order Altman  _et al._
(2015), in contrast to its equilibrium analog; see also Refs. Winfree (1980);
Arenas  _et al._ (2008); Sakaguchi  _et al._ (1987); Daido (1988); Hong  _et
al._ (2005). Can number conservation with mobility restore order, and if so,
at what conditions and what is the nature of this “order”? Further,
biologically relevant rather ubiquitous experimentally realizable examples
with possible ordered phases, which motivate our study, include “active
carpets” of cilia or bacterial flagella, in which active carpets of cilia F.
Meng and Golestanian (2021) or bacteria Uchida and Golestanian (2010a, b),
move chemical substances or other materials, such as sperm in fallopian tubes
Li and Winuthayanon (2017) and mucous in the respiratory tract D. R. Brumley
and Goldstein (2012), trapped near the carpet. In vitro magnetic cilia carpets
Gu  _et al._ (2020) is another related example that inspires the present
study. We are also motivated by some recent studies Uriu  _et al._ (2010);
Uriu and Morelli (2014), which suggest that cell movements can indeed promote
synchronization (i.e., phase-ordering) of coupled genetic oscillators.
Similarly, a collection of interacting mobile oscillators reveal the
significant effects of the oscillator mobility on the steady states, with the
possibility of synchronization at $2d$ Frasca  _et al._ (2008); Peruani  _et
al._ (2010); Levis  _et al._ (2017, 2019). Although it is generally believed
that elasto-hydrodynamic interactions are crucial in synchronization of
biological flagella Elfring and Lauga (2009); Golestanian  _et al._ (2011);
Chakrabarti and Saintillan (2019), recently, other mechanisms like “cell-
rocking” Geyer  _et al._ (2013) and phase-dependent forces Bennett and
Golestanian (2013) are argued to be crucial in the synchronization of swimming
unicellular green alga Chlamydomonas, whereas direct hydrodynamic interactions
are shown to contribute insignificantly. In a recent study on a linearized
hydrodynamic model of diffusive oscillators Banerjee and Basu (2017), it has
been found by using a linearized treatment that a very large oscillator
diffusivity can either promote or destroy synchronization. Lastly, it remains
a paramount theoretically unresolved issue as to how to suppress the
instability of the $2d$ KPZ equation, that is formally related to the $2d$
driven XY model. All these results call for systematic, generic study on the
order and steady states of a collection of active XY model with number
conservation in $2d$.

Finding generic conditions either favorable or disadvantageous to the phase-
ordering in a $2d$ active XY model with a conservation law remain
theoretically challenging, but poorly understood till now. The fundamental
issue here is the existence (or lack of it) of ordered states of the active
$2d$ XY model: Given the applications of the active $2d$ XY model as a
possible basis for theoretical understanding of wide-ranging phenomena of
diverse origin, it would be clearly desirable to explore the concept of
ordering and the extent of universality in the active $2d$ XY model, and the
role of drive on them, and find the scaling laws for such a system. We do so
below by focusing on the specific case of a $2d$ collection of active phase-
ordered XY spins in contact with a conserved density. We show such a system
can synchronize without any long-range hydrodynamic interactions, and both the
local rotation of the spins, and phase difference-dependent current of the
conserved species play a crucial role in synchronization.

In this article, we consider a nearly phase-ordered active XY model
interacting with a conserved, diffusively mobile species on a $2D$ substrate.
We study the stability of small fluctuations around uniform states in this
model, and the scaling properties of these fluctuations in the stable states.
Inspired by the phenomenologies of bacterial quorum sensing and
synchronization in response to complex and dynamic environments Monte  _et
al._ (2007); Mukherjee and Bassler (2019); Schwab  _et al._ (2012), we assume
the microscopic phase and density dynamics to depend on the local phase and
density inhomogeneities in the model. The activity of the XY spins arises from
their tendency to rotate actively in response to local number density of the
diffusive species and magnitude of the nearest-neighbor phase differences,
processes which have no equilibrium analogs, and by local alignments with
neighbors. Further, the conserved current depends on the local phase
differences. From the stand point of active matter systems, our model forms a
new class of “polar” active matter, known as active spinners without directed
walking. These are distinct from the celebrated “moving XY model” for polar
ordered active fluids on substrates Toner and Tu (1995, 1998); Vicsek  _et
al._ (1995). The latter model is characterized by “self-propulsion” of the
particles. This is absent in our model; instead, the particles can rotate
actively, resulting into stable phase-ordered NESS with novel parameter-
dependent scaling behavior, or lack of phase-ordered states with distinct
density morphologies. See Fig. 1 for a schematic model diagram.

Figure 1:  A schematic diagram demonstrating the different local microscopic
dynamical update rules in an active model of 2D XY spins interacting with a
conserved, mobile species on a substrate. The microscopic processes involve
five distinct local processes. The phase changes by local alignments due to
the spin stiffness (process marked $1$), and active rotations proportional to
functions of the local phase differences ($2$ & $2^{\prime}$) and local
concentrations ($3$ & $3^{\prime}$). Unlike the local alignments in $1$, which
always try to minimize the phase differences between neighbouring spins, the
phase updates caused by local phase inhomogeneities in $2$ & $2^{\prime}$
introduce a net rotation of the spins, which either decrease ($2$) or increase
($2^{\prime}$) the phase differences between neighboring spins. The dashed-
outlined arrows represent the updated phase of the spins in response to the
aforementioned factors. The mobile species moves by diffusion ($4$), and also
in response to local phase heterogeneities ($5$).

We study our active XY model by Monte-Carlo simulations (MCS). To understand
the results from the MCS studies theoretically, we write down the hydrodynamic
theory for this system based on symmetry principles and conservation laws, and
analyze and solve the dynamical equations for the hydrodynamic variables.

In interesting connections, our hydrodynamic theory also formally applies to a
generalized KPZ-like surface growth model given by a single-valued non-compact
height field, coupled to the conserved species living on the surface, or an
active fluid membrane with a finite surface tension but without momentum
conservation having a conserved species living on it.

Our MCS studies of the agent-based model, backed by our hydrodynamic theory,
predict that this system of active XY spins can be in stable phase-ordered
states over wide-ranging model parameters. We attribute the ordering to the
interplay between the mobility and the number conservation of the mobile
species, and therefore not found in its immobile analog, or without number
conservation. For wide-ranging choices of the model parameters, the variance
of the phase fluctuations is found to grow slower or faster logarithmically
with the system size compared to the $2d$ equilibrium XY model. The density
fluctuations concomitantly are either miniscule (i.e., hyperuniform) or giant,
compared to their equilibrium counterparts in a gas with short range
interactions at a non-critical temperature, in the steady states. We also find
breakdown of conventional dynamic scaling by logarithmic modulations of the
diffusive dynamic scaling. We further find that depending upon the model
parameters, the hydrodynamic theory can get unstable by two distinct ways,
viz. linear instability, and nonlinear effects, corresponding to the existence
of two distinct types of disordered states in the agent-based model,
characterized by the presence or absence of stable density clusters.

A short account of this problem and the principal results are available in
Ref. Haldar  _et al._ (2021). Detailed derivations and many additional results
are presented here.

The rest of this article is organized as follows. We start with a short
summary of the technical results, which include a brief overview of the
simulations of the agent-based model, and also the governing hydrodynamic
equations themselves, in Sec. II. We set up the agent-based discrete model in
Sec. III. Next, in Sec. IV, we present detailed numerical results on the
phase-ordered states and the disordered states from our simulations of the
agent-based model. In Sec. V, we derive the hydrodynamic equations. Then in
Sec. VI, we study the linearized equations, after expanding about the mean
density $c_{0}$, and obtain the scaling properties of the correlation
functions. Next in Sec. VII.1, we set up the perturbative renormalization
group (RG) calculations to study the nonlinear effects going beyond the
linearized approximations. In Sec. VII.2, we extract the scaling properties
and discuss the nature of order in the nonlinear theory by using perturbative
RG methods, and calculate the various correlation functions in real space in
the equal-time limit. In Sec. VII.3, we discuss the possible routes to
destabilization of the ordered states in the hydrodynamic theory. We further
discuss the phase diagram of the hydrodynamic model in Sec. VII.4. In Sec.
VIII, we briefly discuss and speculate on the role of the topological defects
in the steady states. Finally in Sec. IX, we conclude this paper with a
summary and discussion of the results. Some of the technical details and
several additional results are given in Appendix for the interested readers.

##  II Summary of the technical results

Before we summarize our main technical results, it is instructive to first
briefly recall the $2d$ XY model in equilibrium, which is used as a reference
to our results on the $2d$ active XY model. At $2d$, the classical XY model in
thermal equilibrium is quite remarkable in not having a low temperature phase
with a long-range order Mermin and Wagner (1966). Instead it has a low
temperature phase with QLRO, that undergoes a phase transition to a disordered
state or a state with short-range order (SRO) via unbinding of topological
defects Kosterlitz and Thouless (1973); Chaikin and Lubensky (1995). In the
QLRO phase, the variance of the fluctuations in the local phase $\theta({\bf
x},t)$ in the XY model at $2d$ is

$\langle\theta^{2}({\bf
x},t)\rangle_{eq}=\frac{T}{2\pi\kappa_{0}}\ln\left(\frac{L}{a_{0}}\right),$
(1)

concomitant with an algebraic decay of the spin correlation functions Chaikin
and Lubensky (1995). Here, $T$ is the temperature, $\kappa_{0}$ the spin
stiffness, $L$ the system size, and $a_{0}$ the small scale microscopic
cutoff, and we have set the Boltzmann constant $k_{B}=1$;
$\langle...\rangle_{eq}$ implies averages over thermal noises or thermodynamic
averages. Equation (1) also holds for a fluctuating Edward-Wilkinson (EW)
surface Barabási  _et al._ (1995), an equilibrium dynamics of a surface.
Furthermore, in an equilibrium gas at $2d$ away from any critical points, the
standard deviation $\sigma(N_{0})$ of the number fluctuations in an open area
having $N_{0}$ particles on average follow the well-known “square root $N$”
law Pathria and Beale (2011)

$\sigma(N_{0})\propto\sqrt{N_{0}}.$ (2)

Equations (1) and (2) serve as the benchmark of the results that we derive
here from our model.

In this work, we introduce an agent-based model with novel microscopic
dynamics to systematically study phase ordering, disorder and clustering. We
consider a $2d$ collections of active rotating XY spins of fixed length (or
rotors) on a square lattice with periodic boundary conditions and a mean
number density $c_{0}\equiv N_{\text{tot}}/L^{2}$ of a diffusively mobile
species with local concentration $c_{\ell}$ at site $\ell$. The microscopic
dynamics of the XY spins and the conserved particles consist of simple rules
on the phase updates of the spins together with random fluctuations, and on
the positions of the particles with random fluctuations. The general forms of
the update rules of density $c_{\ell}$ and phase $\theta_{\ell}$ at site
$\ell$ must be invariant under spatial translation and rotation, and
individual rotations of the XY spins by $2\pi$. Motivated by bacterial quorum
sensing and synchronization in complex and dynamic environments Monte  _et
al._ (2007); Mukherjee and Bassler (2019); Schwab  _et al._ (2012) we assume
the phase updates to depend upon the environmental inhomogeneities, e.g.,
local phase inhomogeneities and local concentration of the conserved species,
and also on the lack of clockwise-anticlockwise symmetry due to individual
spin rotations. Position update rules of the conserved species depend upon the
local density inhomogeneities as well as local phase inhomogeneities.

Our simulations clearly reveal the existence of distinct ordered and
disordered states.

To quantitatively characterize the steady states, we calculate (a) the
variance of the orientation fluctuations
$\langle\theta^{2}\rangle\equiv\langle\sum_{\ell}(\theta_{\ell}-\overline{\theta}(t))^{2}\rangle/L^{2}$,
$\overline{\theta}(t)\equiv\sum_{\ell}\theta_{\ell}(t)/L^{2}$, and (b) the
standard deviation $\sigma(N_{0})\equiv\sqrt{\langle N^{2}\rangle-\langle
N\rangle^{2}}$, where $N$ is the total number of particles in an open square
box of linear length $<L$, and $N_{0}\equiv\langle N\rangle$ is its average.
Remarkably, the stable, phase-ordered states show a variety of scaling:

(i) Our systematic MCS studies show that depending upon the model parameters,
the variance of the phase fluctuations $\langle\theta^{2}\rangle$ can grow
with the system size $L$ faster or slower than $\ln L$, i.e., faster or slower
than QLRO or the equilibrium limit result. Our subsequent analytical
approaches to the problem suggests $\langle\theta^{2}\rangle_{L^{2}}\sim(\ln
L)^{\gamma_{1}}$. Here, the exponent $\gamma_{1}$ is found to be positive, and
vary with the model parameters that determine the various active dynamical
processes. The dependence of the scaling exponents on the model parameters is
truly novel and unique feature of this model. Qualitatively speaking, although
like QLRO [see Eq. (1) above], $\langle\theta^{2}\rangle$ grows with the
system size $L$, it does so either logarithmically slower than QLRO (i.e.,
stronger order than QLRO or “SQLRO”) for $0<\gamma_{1}<1$, or faster than QLRO
(i.e. weaker order than QLRO or “WQLRO”) for $\gamma_{1}>1$, depending upon
the model parameters.

(ii) Similarly, the standard deviation of the density fluctuations
$\sigma(N_{0})$ can be bigger or smaller than $\sqrt{N_{0}}$, which holds in
the equilibrium limit. Again, our analytical theory reveals that
$\sigma(N_{0})\sim\sqrt{N_{0}/(\ln N_{0})^{\gamma_{2}}}$. The exponent
$\gamma_{2}$, that can be positive or negative, varies with the same model
parameters as for $\gamma_{1}$. Thus qualitatively speaking, like the
equilibrium “square root $N$-law” or normal number fluctuations (NNF)
$\sigma(N_{0})$ grows with $N_{0}$, it does so either logarithmically slower
than NNF (i.e., fluctuations suppressed in comparison to the equilibrium
limit, known as miniscule number fluctuations (MNF) or hyperuniformity), or
logarithmically faster than NNF (i.e., enhanced fluctuations in comparison to
the equilibrium limit, known as giant number fluctuations (GNF)).

(iii) In a surprising correspondence between the phase and density
fluctuations, our agent-based model shows simultaneous occurrence of SQLRO and
MNF (WQLRO and GNF).

(iv) Our simulations further detect two distinct types of disordered states,
i.e., states with short-range phase order (SRO), distinguished by the presence
or absence of stable density clusters of typical sizes.

In order to understand and rationalize the results from MCS studies of the
agent-based model, we construct the hydrodynamic theory of the nearly phase-
ordered states of an active XY model coupled with a diffusively mobile
conserved density. One of the hydrodynamic variables reflects the
“orientational order” in the order parameter space, or “phase-order” of the XY
spins. This is described by the non-conserved, broken symmetry field or
Goldstone mode $\theta({\bf r},t)$ Chaikin and Lubensky (1995), a suitable
coarse-grained version of individual spin (particle) based orientation of the
spins with respect to some (arbitrary) reference axis. Our second hydrodynamic
variable is the conserved density fluctuations $\delta c({\bf r},t)\equiv
c({\bf r},t)-c_{0}$ of the surface density $c({\bf r},t)$ of a diffusing
species about its mean $c_{0}$. In our hydrodynamic theory, we study small
fluctuations around a uniform reference state with $\theta({\bf
r},t)=\theta_{0}$, measured with respect to an arbitrary reference axis, and
density $c({\bf r},t)=c_{0}$. The hydrodynamic equations for these variables
will be systematically derived in Sec. V below. These are obtained by gradient
expansions around the uniform reference state, retaining the lowest order
symmetry-permitted nonlinear terms. These equations are

$\displaystyle\frac{\partial\theta}{\partial
t}=\kappa\nabla^{2}\theta+\Omega_{1}\delta
c+\frac{\lambda}{2}({\bm{\nabla}}\theta)^{2}+\Omega_{2}(\delta
c)^{2}+f_{\theta},$ (3a) $\displaystyle\frac{\partial\delta c}{\partial
t}=\lambda_{0}\tilde{\Omega}_{0}\nabla^{2}\theta+D\nabla^{2}\delta
c+\lambda_{0}\tilde{\Omega}_{1}{\bm{\nabla}}\cdot(\delta
c{\bm{\nabla}}\theta)+f_{c}.$ (3b)

Here,
$\kappa>0,D>0,\Omega_{1},\lambda,\Omega_{2},\lambda_{0},\tilde{\Omega}_{0},\tilde{\Omega}_{1}$
are phenomenological parameters of our model. Physical interpretations of the
different terms in (3a) and (3b) are discussed in Sec. V below. Equations (3a)
and (3b) may also be interpreted as the coupled hydrodynamic equations for a
local single-valued height field $\theta$, measured with respect to an
arbitrary base plane, and a conserved density $\delta c$ that exists on it;
see also below.

The Gaussian noises $f_{\theta}$ and $f_{c}$ have zero mean and variances

$\displaystyle\langle f_{\theta}({\bf
x},t)f_{\theta}(0,0)\rangle=2D_{\theta}\delta^{2}({\bf x})\delta(t),$ (4)
$\displaystyle\langle f_{c}({\bf
x},t)f_{c}(0,0)\rangle=2D_{c}(-\nabla^{2})\delta^{2}({\bf x})\delta(t).$ (5)

The hydrodynamic equations (3a) and (3b) together with (4) and (5) form the
active spinner analog of the celebrated Toner-Tu model Toner and Tu (1995) for
dry active flocks.

We now list the principal results obtained by using (3a) and (3b) together
with (4) and (5). It predicts either stable uniform NESS with diverse scaling
properties of small fluctuations around them, or novel instabilities
signalling the loss of any order. We show that the hydrodynamic equations that
we derive here when linearized about $c_{0}$, the mean particle number
density, can be linearly stable or unstable. We further show that in the
linearly stable case, the NESS of such a system admits underdamped propagating
waves having linear dispersion with the wavespeed being proportional to
$c_{0}$. If linearly unstable, then the growth rate of the modes is
proportional to the wavevector. Focusing on the NESS for large system sizes
when the propagating modes clearly dominate over damping in the long
wavelength limit, we show that the NESS can be characterized by the variances
of the phase fluctuations, density fluctuations and the associated correlation
functions. The stable states in the linear theory have scaling laws for the
phase and density fluctuations indistinguishable from its equilibrium
counterparts. Nonetheless, the system remains out of equilibrium even in the
linear theory as can be seen from the non-zero phase-density cross-correlator
in the linear theory. The nonlinear effects either leave scaling of the
linearly stable states essentially unchanged modulating only by introducing
logarithmic corrections, or make the linearly stable states unstable,
destroying any ordered states. Nonetheless, these logarithmic corrections in
the stable case actually make the scaling regime rich and diverse, either by
strengthening or weakening the phase order vis-à-vis the QLRO phase order in
equilibrium or as predicted by the linear theory, and at the same time
suppressing or enhancing standard deviations of the density fluctuations vis-
à-vis again the equilibrium limit result of NNF, or as predicted by the linear
theory. We obtain these results by applying one-loop dynamic renormalization
group methods on the hydrodynamic equations. These results are consistent with
those obtained from the MCS studies on the agent-based model.

The NESS of the nonlinearly stable states are characterized by a variety of
quantities. For instance, the variance of the phase fluctuations

$\displaystyle\Delta_{\theta}^{R}\equiv\langle\theta^{2}({\bf
x},t)\rangle_{R}\approx\frac{\overline{D}(0)}{4\pi\Gamma(0)}\left[\ln\left(\frac{L}{a_{0}}\right)\right]^{1-\eta},$
(6)

where $\eta$ is a scaling exponent that surprisingly depends on the
dimensionless ratios $\mu_{1},\,\mu_{2}$ of the “bare” or unrenormalized
nonlinear coupling constants in the hydrodynamic equations of motion that we
derive later; $R$ refers to renormalized quantities. Further, $\eta$ can be
positive or negative, with $\eta<1/3$ always. Furthermore, $\overline{D}(0)$
and $\Gamma(0)$ are the unrenormalized (i.e., the microscopic or bare) values
of the effective noise strength and damping coefficient. Since $\eta$ can be
positive or negative, $\Delta_{\theta}$ can grow, respectively, slower or
faster than $\ln(L/a_{0})$, its form in QLRO or the equilibrium limit. When
$\theta({\bf x},t)$ interpreted as a single-valued height field, (6) implies a
fluctuating surface that is logarithmically smoother or rougher than an EW
surface, respectively, for $\eta>0$ and $\eta<0$. Similarly, the equal-time
renormalized phase-difference correlation function $C_{\theta\theta}^{R}(r)$
for the XY spins, an exact analogue of the phase-difference correlation
function in the XY model, scales as

$C_{\theta\theta}^{R}(r)\equiv\langle[\theta({\bf x+r},t)-\theta({\bf
x},t)]^{2}\rangle_{R}\approx\frac{\overline{D}(0)}{2\pi\Gamma(0)}[\ln(r/{a_{0}})]^{1-\eta},$
(7)

in the limit of large $r$. This means, the renormalized spin correlation
function $C_{ZZ}^{R}(r)$, defined analogously with the spin-spin correlation
function for the usual (equilibrium) XY model (see later for a formal
definition), decays slower or faster for $\eta>0$ or $\eta<0$, respectively,
than its equilibrium limit (i.e., $\eta=0$), as we show later in details. Once
again, for a fluctuating surface, Eq. (7) means the equal-time height
difference correlation function decays logarithmically slower or faster than
in an EW surface as the separation $r$ grows. The above results on phase
fluctuations clearly indicate order logarithmically stronger than QLRO for
$\eta>0$ and logarithmically weaker than QLRO for $\eta<0$. We call these
SQLRO and WQLRO phase order, respectively. Clearly, for $\eta=0$, QLRO is
retrieved, which is identical to the equilibrium results.

We further obtain the equal-time renormalized correlation function of the
density difference $\delta c({\bf x},t)\equiv c({\bf x},t)-c_{0}$, where
$c_{0}$ is the mean density. We find in the Fourier space,
$C_{cc}^{R}(k)\equiv\langle|\delta c({\bf k},t)|^{2}\rangle_{R}$ picks up a
weak $k$-dependence displaying unusual scaling behavior:

$C_{cc}^{R}(k)=\frac{\overline{D}(0)}{2\Gamma(0)}\left[\ln\left(\frac{\Lambda}{k}\right)\right]^{-\eta}$
(8)

in the hydrodynamic limit $k\rightarrow 0$; where $\delta c({\bf k},t)$ is the
Fourier transform of $\delta c({\bf x},t)$, and $\Lambda=2\pi/a_{0}$ is an
upper wavevector cutoff. Clearly, for $\eta>0$, the density fluctuations are
strongly suppressed in the long wavelength limit $k\rightarrow 0$ vis-à-vis
the equilibrium result for a gas with short range interactions away from any
critical points. This is a manifestation of miniscule number fluctuations
(MNF) (also known as hyperuniformity), an exotic state of matter Torquato and
Stillinger (2003); Torquato (2018) rarely encountered in ordered active matter
Oppenheimer  _et al._ (2022). On the other hand, if $\eta<0$ the density
fluctuations are hugely enhanced in the long wavelength limit, in comparison
with NNF. These are giant number fluctuations (GNF), often encountered in the
context of orientationally ordered active fluids Marchetti  _et al._ (2013),
and also in equilibrium superfluids Chaikin and Lubensky (1995); $\eta=0$
corresponds to the equilibrium result, as before. Furthermore, the equal-time
renormalized density autocorrelation function
$C^{R}_{cc}(r)\equiv\langle\delta c({\bf x+r},t)\delta c({\bf
x},t)\rangle_{R}$, which is just the inverse Fourier transform of
$C_{cc}^{R}(k)$, has the form

$\displaystyle C^{R}_{cc}(r)$ $\displaystyle\equiv\langle\delta c({\bf
x+r},t)\delta c({\bf x},t)\rangle_{R}$
$\displaystyle\approx\frac{\overline{D}(0)}{4\pi\Gamma(0)}\frac{-1}{r^{2}}\eta[\ln(\Lambda
r)]^{(-1-\eta)},$ (9)

for large $r$, $r\Lambda\gg 1$. Thus for $\eta>(<)0$, $C^{R}_{cc}(r)$ falls of
relatively faster (slower). This further implies that the standard deviation
of the number fluctuations $\sigma(N_{0})$ in an area of linear size
$\tilde{L}$ containing $N_{0}$ particles on average (i.e., $N_{0}=\langle
N\rangle$ as measured in the area $\tilde{L}\times\tilde{L}$) is

$\sigma(N_{0})\propto\sqrt{N_{0}/(\ln N_{0})^{\eta}}<(>)\sqrt{N_{0}},$ (10)

for $\eta>(<)0$. Noting that $\sigma(N_{0})\propto\sqrt{N_{0}}$ is the well-
known “square root mean” law of fluctuations of equilibrium statistical
mechanics, which are “normal number fluctuations” (NNF), expected in
equilibrium systems with short range interactions away from any critical
point, (10) again implies MNF, i.e., hyperuniform density and GNF,
respectively, for $\eta>(<)0$. This is consistent with the discussions on MNF
and GNF in this model above. NNF is retrieved when $\eta=0$. We thus show that
SQLRO phase-ordering is associated with MNF, whereas WQLRO phase-ordering
comes with GNF. This establishes a novel one-to-one correspondence between
phase and density fluctuations in the ordered states, and forms a principal
conclusion from this study.

We also find breakdown of conventional dynamic scaling in the stable ordered
states, where time-scale $t$ no longer shows simple scaling with length-scale
$r=|{\bf r}|$ via $t\propto r^{z}$, but is modulated by logarithmic
corrections

$t\propto r^{2}/[\ln(r/a_{0})]^{(1-\eta)/2},$ (11)

for large $r$. Therefore, the dynamics is necessarily faster than ordinary
diffusion, albeit only logarithmically. This remains true even when $\eta=0$,
i.e., when QLRO is observed. These results define the scaling behavior of this
model. The model parameter-dependence of the scaling exponents is a key novel
aspect, that reflects breakdown of the conventional notion of universality
here.

In other regions of the phase space, no uniformly ordered states are found.
The hydrodynamic theory predicts that in the parameter regimes corresponding
to no stable uniform order, uniformly ordered states can get destabilized by
two distinct mechanism - (a) due to linear instability, and (b) nonlinear
instability of the linearly stable states. Results from our numerical
simulations of the agent-based model also reveal the existence of two distinct
kinds of the disordered states: (a) one kind has an intrinsic length-scale
that should control the morphology of the eventual steady states, whereas (b)
disordered states of the second type are visually similar to the rough phase
of the $2d$ KPZ equation, lacking any intrinsic length-scale.

In the immobile limit with uniform density, or in the absence of number
conservation, e.g., due to birth-death events, the number density fluctuations
relax fast, and is not a hydrodynamic variable. Therefore, the hydrodynamic
equation for $\theta({\bf x},t)$ reduces to the well-known Kardar-Parisi-Zhang
(KPZ) equation for surface growth phenomena Altman  _et al._ (2015). The $2d$
KPZ equation does not have a stable “smooth phase” (equivalent to a phase
ordered state in the spin language), but has only an algebraically rough phase
Barabási  _et al._ (1995), which for phase $\theta$ should correspond to short
range order (SRO), in agreement with the results from the agent-based model.

The hydrodynamic theory predictions provide a generic theoretical framework to
understand the numerical results from the agent-based model. In particular,
both the agent-based model and the continuum hydrodynamic theory confirm the
variety of scaling of the numerical analogs of $\Delta_{\theta}$ with
$\ln(L/a_{0})$, and the standard deviation $\sigma(N_{0})$ with $N_{0}$,
demonstrating NESS with order both stronger and weaker than QLRO, together
with MNF and GNF, respectively. This reinforces the hydrodynamic theory
predictions on the occurrence of SQLRO with MNF and WQLRO with GNF, the
principal qualitative feature of scaling in the phase-ordered states of this
model. We further claim that the existence of two types of the disordered
states of the agent-based model can be argued by considering the instabilities
of the phase-ordered states in the hydrodynamic theory.

Our theory should form the basis of understanding of a wide class of driven
systems, e.g., active superfluidity in bacterial suspensions on substrates,
active carpets of cilia and their artificial imitations, and mobile oscillator
synchronization. Our hydrodynamic equations (3a) and (3b) also describe an
active fluid membrane with a finite tension and a conserved species living on
it and without momentum conservation Cagnetta  _et al._ (2020). It is also
intriguing from nonequilibrium statistical mechanics perspectives by showing
how the well-known instability of the $2d$ KPZ equation can be suppressed by
coupling it to a conserved density.

##  III Agent-based model

In this Section, we construct an ab initio agent-based lattice-gas type model
and study it by Monte-Carlo simulations. We use these MCS studies to
characterize the phase-ordered states, and extract the scaling exponents for
both density and phase fluctuations. The agent-based model is further used to
systematically study the disordered phases.

We consider a collection of XY spins, each grafted on a lattice site of a two-
dimensional square lattice of size $L^{2}$, interacting with a conserve,
mobile particles of mean density $c_{0}=N_{tot}/L^{2}$; where $L$ is the
linear system size along the X and Y directions, and $N_{tot}$ is the total
number of mobile particles. On the lattice, a site with $x$\- and
$y$-coordinates is conveniently labeled by a single index $\ell\in[1,L^{2}]$
defined as

$\ell=(y-1)L+x,$ (12)

where $x\in[1,L]$, $y\in[1,L]$. Each site has four nearest neighbors and we
impose periodic boundary conditions along the $x$ and $y$ directions.

The initial condition of the system is then prepared by assigning each site
$\ell$ a fixed number of $c_{0}$ particles and a unit spin vector
${Z}_{\ell}=(\cos\theta_{\ell},\sin\theta_{\ell})$, where
$\theta_{\ell}=\phi\pm r_{\ell}\Delta\phi$, and $\phi\in(-\pi,\pi)$ is same
for all sites, $r_{\ell}\in(0,1)$ is a uniform random number and
$\Delta\phi=\pi/36$. Therefore, the initial configuration of the system is
nearly phase ordered with a given uniform density of the mobile species.

Simulation Details:- Now let us discuss the principal facets of the update
rules that govern the system dynamics; a comprehensive details of the same are
presented in Appendix B. We use a square lattice for simplicity although the
results should be independent of the underlying lattice structure in the long
wavelength limit, as corroborated by our hydrodynamic theory later. We impose
periodic boundary condition. To simulate the model, we make use of the
canonical sampling Monte Carlo (MC) method where in a single MC step, we
sequentially update both the phase and density of the $L^{2}$ sites where the
sites to be updated are chosen in a stochastic manner. Therefore, certain
sites may be chosen more than once, and certain sites may not be updated in a
single MC step.

As mentioned above, inspired by the phenomena of quorum sensing by bacteria
Monte  _et al._ (2007); Mukherjee and Bassler (2019); Schwab  _et al._ (2012),
in which bacteria may try to synchronize and rotate collectively in response
to dynamically changing environmental inhomogeneities (which in the present
model should be just local phase and density inhomogeneities), we define the
update rules of the phases.

(i) Local alignment: This is present in the equilibrium limit as well, and
involves the relaxation of the spins to align themselves locally with the
neighbors. We model this term, by updating the local phase of site along the
average orientation of the nearest neighbors including the site and then
adding a small noise to it that describes the stochastic nature of the
microscopic dynamics (see also Ref. Vicsek  _et al._ (1995) in the context of
active polar flocks on a substrate).

(ii) Phase difference-dependent local active rotation: This is the local phase
difference-dependent rotation of the spins, in response to the phase
inhomogeneities with the nearest neighbor; see also e.g. Ref. Mukherjee and
Bassler (2019). In the absence of any self-propulsion, we take it to be
proportional to a function of the local phase differences, that is generally
arbitrary, constrained only by the symmetry of the model under spatial
translation, rotation and rotation of any spin by $2\pi$, and should break the
clockwise-anticlockwise symmetry for (active) spin rotations. For simplicity,
we choose them to be
$g_{1}\cos(\sum_{\ell^{\prime}}|\theta_{\ell}-\theta_{\ell^{\prime}}|/4)$,
where $\ell^{\prime}$ is a nearest neighbor site of the randomly chosen site
$\ell$; $g_{1}$ measures the relative amplitudes of this active processes. The
parameter $g_{1}$, a simulation parameter, controls the (relative) strength of
this process in the overall update of the phase.

(iii) Number density dependent local active rotation: Lastly, the spins can
undergo local number-dependent rotation; see, e.g., Ref. Monte  _et al._
(2007) for a similar mechanism. This term is modeled by assuming the local
rotation to be proportional to $g_{2}\Omega(c_{\ell})$, taken as a quadratic
function of its argument for simplicity; where $c_{\ell}$ represents the
instantaneous number of interacting mobile particles at site $\ell$.
Simulation parameter $g_{2}$ controls the (relative) amplitude of this
particular process in the overall update of the phase.

Finally, the phase of the random site is updated probabilistically where with
probability $p_{\theta}$, we have a new phase having contributions from the
local alignment (discussed in (i)) and from the phase and density fluctuations
discussed in (ii) and (iii). With probability $1-p_{\theta}$, we only allow
update of $\theta$ from the equilibrium relaxation due to the nearest neighbor
orientations as described in (i). By varying $p_{\theta}$, we can control the
relative strength of the microscopic active processes in the update of the
phases. The probabilistic update is used to encode the inherent stochasticity
of the noisy dynamics.

For both the above mentioned update scenarios, the updating of the phase is
also supplemented by a random noise $\xi$ drawn from a uniform distribution
$[-\xi/2,\xi/2]$, which has an effect analogous to thermal fluctuations
measured by the temperature in equilibrium systems. Here, let us emphasize the
well known analogy between the flocking transition described by the Vicsek
Vicsek  _et al._ (1995) model and the ferromagnetic transition in the Ising
model. The role played by the random noise $\xi$ associated with the alignment
in the Vicsek model is equivalent to the role of temperature in the
ferromagnetic transition of the Ising model; for less $\xi$ the system is
ordered, while the system is disordered for higher $\xi$. The analogy between
the noise $\xi$ and temperature $T$ in equilibrium could also be drawn from
Eq. (1). In the QLRO phase of the equilibrium XY model, the variance of
fluctuations of the local phase $\langle\theta^{2}({\bf x},t)\rangle$ scales
with $T$. It may be noted that our agent-based model reduces to the
equilibrium XY model in the passive limit (i.e., when the “active” or
nonequilibrium ingredients in the phase update rules are dropped); see Fig. 12
in Appendix C and the discussions there.

Likewise, the time-evolution of the local number density occurs via two
different processes:

(i) Random hopping in response to local density inhomogeneities: Particles hop
randomly to the nearest neighbor sites with a lower density, which models
ordinary diffusion and

(ii) By hopping in response to the phase difference between the originating
site $\ell$ and a randomly chosen nearest neighbor target site $v$. This is
reminiscent of the phase-dependent force introduced in a microscopic model of
synchronization in Chlamydomonas Bennett and Golestanian (2013). Again, this
function is arbitrary, being constrained only the above symmetries. Further,
allowing the phase difference-dependent current to depend on the sense of
rotation, we take this function to be $c_{\rm sum}=\Omega^{\prime}(\langle
c\rangle_{v})\sin(\theta_{\ell}-\theta_{v})$, for simplicity, where $\langle
c\rangle_{v}$ is the mean density of sites $\ell$ and $v$. Function
$\Omega^{\prime}(\langle c_{v}\rangle)$ is assumed to be a quadratic function
of its arguments for simplicity that contains two parameters $c_{s}$ and
$c_{s2}$; see Eq. 145 in Appendix B. At each MC step, we calculate $c_{\rm
sum}$: For $c_{\text{sum}}=\pm 1$, particles move from (to) the originating
site to (from) the target site.

In the same spirit of $p_{\theta}$, we define another probability $p_{c}$,
which dictates whether the particle hopping to its nearest neighbor be
determined by pure diffusion with probability $1-p_{c}$ (as described in (i))
or by the phase difference-induced current of the particles with probability
$p_{c}$ (as discussed in (ii)).

It is instructive to compare our agent-based model with the Vicsek model for
flocking Vicsek  _et al._ (1995); Ginelli (2016). In the latter model, the
particles or agents are “self-propelled’, i.e., locally align themselves with
the neighbors, and move along that direction, resulting into flocking at low
noises and high densities. In contrast, in the present model, the spins are
immobile and locally align with the neighbors, but obviously do not move in
that direction. Instead, they actively rotate locally in response to the local
phase difference in a way that breaks the clockwise-anticlockwise symmetry,
and to the local number density as well.

##  IV Monte-Carlo simulations of the agent-based model

###  IV.1 General features

The model which we have defined in Sec. III, is simulated on a square lattice
with linear sizes $L=32,64,128,192$ and $256$ with $c_{0}=5$. We have fixed
the lattice spacing $a_{0}=1$. Starting from a nearly phase-ordered spin
configurations, the Monte Carlo algorithm (see Appendix B) evolves the system
under various control parameters ($e.g.,$ $g_{1}$, $g_{2}$, $c_{s}$, $c_{s2}$,
and $\xi$) until steady states are reached. In our simulations, $\xi=0.1$
(unless otherwise mentioned).

During the simulations, we first let the system to evolve for $t_{ss}=10^{4}$
MC step, so that steady states are reached, and following this, measurements
are carried out in the steady states and noise-averaged data of the phase
fluctuation, $\langle\theta^{2}\rangle$, are recorded up to some maximum time
$t_{m}=10^{6}$ MC step. To achieve better statistics, data obtained are
further averaged over $13$ independent initial spin configurations. The
numerical code is implemented to measure the variance of the average spin
fluctuations $\langle\theta^{2}\rangle$ in the following manner:

$\langle\theta^{2}\rangle=\langle(\theta-\bar{\theta}_{L^{2}})^{2}\rangle_{t,L^{2}},$
(13)

where $\bar{\theta}_{L^{2}}$ is calculated using

$\bar{\theta}_{L^{2}}=\arctan[\langle\sin(\theta)\rangle_{L^{2}}/\langle\cos(\theta)\rangle_{L^{2}}],$
(14)

Here, $\langle...\rangle_{L^{2}}$ denotes the average over all the $L^{2}$
sites at a given time $t$, and $\langle...\rangle_{t,L^{2}}$ denotes averaged
over time and $L^{2}$ sites for a fixed $\xi$ (meaning a fixed “nonequilibrium
temperature”). Notably, a nonzero $\bar{\theta}_{L^{2}}$ in the steady state
indicates global rotation of the spins, consistent with the identification of
the model as an active oscillator model.

We numerically calculate $\sigma(N_{0})$ by

$\sigma(N_{0})=\sqrt{\langle(N-N_{0})^{2}\rangle}$ (15)

in the steady state, where $N$ is the total number of particles in an open
square box of increasing linear length $l_{s}<<L$, and $N_{0}\equiv c_{0}\
l_{s}^{2}$ is the average number of particles within that box;
$\langle...\rangle_{t}$ denotes averaged over time.

All the averages are performed in the steady states.

###  IV.2 Steady states in the agent-based model

Our simulations clearly reveal the existence of distinct ordered and
disordered states. The most spectacular results from our simulations are in
fact manifested in the movies MOVIE1, MOVIE2, MOVIE3 and MOVIE4 in the
Supplemental Material sup . Movies MOVIE1 and MOVIE2 respectively represent
the temporal evolution of the spin configurations in the ordered steady states
and states without any ordering, whereas the corresponding movies showing the
time evolution of the density distributions are MOVIE3 (ordered state) and
MOVIE4 (no order).

To motivate the reader further and to give him/her an bird’s eye view of the
ensuing results, we also refer to the snapshots in Fig. 2(a), (b), (d), and
(e) (ordered state), and (c) and (f) (disordered state). In the spin
configuration plots (Fig. 2(a)-(c)), an arrow represents the local coarse-
grained spin orientation at a time after the steady state has reached, whereas
the colorbars in Fig. 2(d)-(f) represent the local density fluctuation. Fig.
2(a)-(b) show a nearly ordered spin configurations for the given parameter
values which is validated by the corresponding density colormaps shown in Fig.
2(d)-(e) where the uniform color gradient signifies a homogeneous ordered
phase. In Fig. 2(a) and (d), we have shown the plots for $L=64$ wheres Fig.
2(b) and (e) are for $L=256$. These figures clearly demonstrate that ordering
persists even with a 4 fold increase in $L$, as expected in ordered states. In
contrast, Fig. 2(c) and (f) signifies the disordered states (for the parameter
values mentioned) as evident from the disorderly orientations of the spin
vectors in (c) and from the density colormap in (f) which displays wide
fluctuations in the local density, including some empty or near empty regions.
Unsurprisingly, the movies mentioned above directly correspond to the model’s
configuration snapshots. For instance, movies MOVIE1 and MOVIE3 correspond to
the snapshots in Fig. 2(a) and (d), which are clearly phase-ordered states,
together with small density fluctuations around the mean density. In contrast,
movies MOVIE2 and MOVIE4 correspond to the snapshots in Fig. 2(c) and (f),
which clearly lack any phase order, and have large density fluctuations.

Figure 2: Representative snapshots showing ordered [(a) and (b)], and
disordered (c) spin configurations, and ordered [(d) and (e)], and disordered
(f) density distributions. The model parameters are mentioned on top of each
figure. For better visual clarity, we have shown the phase configurations
taken from a small part of each system (a $16^{2}$ section is shown); the
density distributions are shown over full system sizes. Subplots (a) and (b)
show ordered configurations for $L=64,\,256$, showing the persistence of order
even after 4 fold increase in $L$; subplot (c) shows clear lack of any order.
Similarly, subplots (d) and (e) show small density fluctuations about the mean
density $c_{0}=5$ in the snapshots for $L=64,\,256$, whereas in subplot (f)
much bigger density fluctuations can be seen, signifying a disordered state.
Other parameters: $c_{0}=5$, $\xi=0.1$, and $p_{\theta}=p_{c}=0.5$.

The ordered and disordered states can be further distinguished by considering
the probability distribution $P(\theta)$ of the orientation
$\theta\in[0,2\pi]$ and the probability distribution $P(c)$ of the density
$c$. The corresponding plots are shown in Fig. 3. In a perfectly ordered
uniform state, $P(\theta)\propto\delta(\theta-\theta_{0})$ and
$P(c)\propto\delta(c-c_{0})$, where $\theta_{0}$ is an arbitrary reference
state $\in[0,2\pi]$. Nonetheless, it is reasonable to expect that in the
ordered states, both $P(\theta)$ and $P(c)$ will continue to be sufficiently
narrow and will have peaks around the respective mean values of $\theta$ and
$c$, as shown in Fig. 3(a) and (b), respectively. In Fig. 3, we observe that
$P(\theta)$ and $P(c)$ are sharply peaked around angle $\theta_{0}\approx
330^{\circ}$ and at the mean density $c_{0}=5$. In general, the mean density,
as obtained from and related to the peak of $P(c)$, should depend on the
simulation parameter for average density. While the specific value of
$\theta_{0}$ (calculated in the rotating frame) is of no physical
significance, the narrow distribution does imply an orientationally ordered
state, accompanied by a nearly smooth distribution of the particle number,
corresponding to $P(c)$ being sharply peaked around $c=c_{0}$. Quite
obviously, the distributions, $P(\theta)$ and $P(c)$, would be sufficiently
broad with larger variance in disordered states as also demonstrated in Fig. 3
where we observe that the phase distribution is flat, indicating all values of
the phase is equally probable, along with a wide distribution of the density
implying large fluctuations of the density. Note that we have also observed
another form of disordered phase in which the density distribution has two
peaks, one near zero and another smaller peak at a density far greater than
the mean density $c_{0}$ (see Fig. 6 below). Two kinds of disordered phases
are discussed in detail later.

Figure 3: Plots of (a) $P(\theta)$ versus $\theta$ and (b) $P(c)$ versus $c$
in the ordered (red) and disordered (blue) states. As argued, in the ordered
state [corresponding to the snapshots in Fig. 2(a) and (d)], both $P(\theta)$
and $P(c)$ are narrow, peaked around their mean values, whereas in the
disordered state [corresponding to the snapshots in Fig. 2(c) and (f)],
$P(\theta)$ is nearly flat, and $P(c)$ is much broader than for the ordered
state (see text).

####  IV.2.1 Ordered states

Are all the ordered states in the agent-based model, parametrized by the
simulation parameters $g_{1},\,g_{2},\,c_{s},\,c_{s2}$ statistically
equivalent? That is to ask: Are all of them characterized by the same scaling
exponents? We now set out to address this generic question by quantitatively
studying these ordered states below by using the agent-based model.

The natures of order in the NESS for the different parameters are investigated
in Fig. 4. Fig. 4(a) and Fig. 4(b) demonstrate instantaneous spin
configurations, and Fig. 4(c) and Fig. 4(d) show the corresponding densities
on a $32^{2}$ segment of a $256^{2}$ lattice for better visualization. The
system exhibits SQLRO/MNF for the active coefficients (a & c) $g_{1}=1.0$,
$g_{2}=0.02$ and WQLRO/GNF for (b & d) $g_{1}=2.0$, $g_{2}=0.03$ respectively,
the quantitative validation of which is provided below.

We have numerically calculated $\langle\theta^{2}\rangle$ from (13) and
extracted its dependence on $\ln L$ and plotted $\langle\theta^{2}\rangle/\ln
L$ versus $\ln L$ in Fig. 5(a). Since in equilibrium
$\langle\theta^{2}\rangle$ scales as $\ln L$ meaning QLRO, we investigate if
the nonlinear effects alter the $\ln L$-dependence of
$\langle\theta^{2}\rangle$ in the active problem. In fact, we detect two
distinctly different kinds of dependence, both in the phase-ordered steady
states: Noting that for QLRO or in the equilibrium limit,
$\langle\theta^{2}\rangle/\ln L$ is a constant independent of $L$, we find
$\langle\theta^{2}\rangle$ necessarily increases with $\ln L$. While our
numerical results emphatically rule out SRO, it is also not LRO, since
$\langle\theta^{2}\rangle$ grows with $\ln L$ (hence with $L$). This behavior
is then reminiscent of QLRO, being a “middle ground” between SRO and LRO.
However, a careful observation of Fig. 5(a) reveals that depending upon the
parameters, there are actually two qualitatively distinctly different scaling
behaviors - in one case $\langle\theta^{2}\rangle$ grows slower than $\ln L$,
and in the other case faster than $\ln L$. We name them stronger than QLRO or
SQLRO and weaker than QLRO or WQLRO, respectively. However, the system
eventually becomes disordered for large enough values of $\xi$. For a detailed
picture of the nature of ordering as a function of the control parameters, we
refer to Table 1.

In parallel to calculating $\langle\theta^{2}\rangle$, we further investigate
the number fluctuations $\sigma(N_{0})$ as defined in (15), and extract its
dependence on $N_{0}$. In equilibrium systems away from any critical point
$\sigma(N_{0})\sim\sqrt{N_{0}}$. To see if active effects can change this
result, we plot $\sigma(N_{0})/\sqrt{N_{0}}$ versus $\sqrt{N_{0}}$ in Fig.
5(b). This, just like $\langle\theta^{2}\rangle$ as a function of $\ln L$,
displays two qualitatively different kind of scaling with $\sqrt{N_{0}}$. In
one case, $\sigma(N_{0})$ grows slower than $\sqrt{N_{0}}$ for model
parameters for which $\langle\theta^{2}\rangle$ shows SQLRO. This case clearly
implies miniscule (with respect to its equilibrium counterpart) number
fluctuations or MNF. For other choices, $\sigma(N_{0})$ grows faster than
$\sqrt{N_{0}}$ for model parameters for which $\langle\theta^{2}\rangle$ shows
WQLRO. This case gives giant (again with respect to its equilibrium
counterpart) number fluctuations or GNF. Our investigation confirms that MNF
and GNF in Fig. 5(b) respectively correspond to SQLRO and WQLRO of Fig. 5(a),
and therefore establishes the SQLRO $\Leftrightarrow$ MNF and WQLRO
$\Leftrightarrow$ GNF correspondence in the agent-based model.

Figure 4: NESS snapshots showing (a-b) the ordered spin configurations and
(c-d) the corresponding instantaneous density profiles in a $32^{2}$ segment
of a system with $L=256$, shown for two different combinations of active
coefficients $g_{1}$ and $g_{2}$: (a & c) $g_{1}=1.0$, $g_{2}=0.02$ and (b &
d) $g_{1}=2.0$, $g_{2}=0.03$ respectively. Parameters: $c_{0}=5$, $\xi=0.1$,
$p_{\theta}=p_{c}=0.5$, $c_{s}=1$, and $c_{s2}=1$. Figure 5:  Plots showing
scaling in the ordered phase obtained from the agent-based simulations: (a)
$\langle\theta^{2}\rangle/\ln L$ versus $\ln L$ showing SQLRO and WQLRO, and
(b) $\sigma(N_{0})/\sqrt{N_{0}}$ versus $\ln N_{0}$ showing MNF and GNF
($L=128$) with parameters : $c_{0}=5$, $\xi=0.1$, $p_{\theta}=p_{c}=0.5$,
$c_{s}=1$, $c_{s2}=1$, (SQLRO/MNF): $g_{1}=1.0$, $g_{2}=0.02$ and (WQLRO/GNF):
$g_{1}=2.0$, $g_{2}=0.03$ respectively. The red broken horizontal lines denote
the behavior of plots in the equilibrium limit, or the (a) QLRO and (b) NNF
correspondence associated with the plots. The SQLRO-MNF and WQLRO-GNF
correspondences are clearly established in the agent-based model (see text).
Figure 6: Plot of $P(c)$ versus $c$ showing two distinct types of disordered
phase obtained in the agent-based model, referred as rough and aggregate
phases. Parameters: $L=64$, $c_{0}=5,\,\xi=0.1$, $p_{\theta}=p_{c}=0.5$,
(‘Rough’): $g_{1}=4.0$, $g_{2}=0.01$, $c_{s}=1$, and $c_{s2}=1$, and
(‘Aggregate’): $g_{1}=-1.0$, $g_{2}=0.01$, $c_{s}=-1$, and $c_{s2}=-1$,
respectively. Types of Order
---
$g_{1}$ |  $g_{2}$ |  Order
$1.0$ |  |  SQLRO
$1.5$ |  $0.01$ |  SQLRO
$2.0$ |  |  SQLRO
$1.0$ |  |  SQLRO
$1.5$ |  $0.02$ |  SQLRO
$2.0$ |  |  QLRO∗
$1.0$ |  |  SQLRO
$1.5$ |  $0.03$ |  QLRO∗
$2.0$ |  |  WQLRO
|  $0.01$ |
$4.0$ |  $0.02$ |  Rough phase
|  $0.03$ |
Table 1: The natures of ordering obtained for different values of the
parameters $g_{1},\;\text{and}\;g_{2}$ respectively. The values of the
coefficients $c_{s}$ and $c_{s2}$ are kept $1$. Other parameters: $c_{0}=5$,
$\xi=0.1$, and $p_{\theta}=p_{c}=0.5$. Increasing the magnitude of $g_{1}$ can
result in a transition from ordered to disordered rough phase in our
simulations (see text).

####  IV.2.2 Disordered states

We now use our agent-based model to study the nature of the disordered states
and characterize them. Surprisingly, our simulation studies reveal the
existence of two distinct types of disordered states, characterized by their
distinctive density distributions (see Fig. 6). In one case, the probability
distribution of the number density is still peaked about $c_{0}$, but is much
wider (and also asymmetric about $c_{0}$). Clearly, a broad density
distribution peaked at $c_{0}$ is found in our simulations implies lack of
stable cluster formations of any preferred sizes. We name this rough phase. In
the second case, that we call as the aggregate phase, it is not peaked about
$c_{0}$; instead, a large peak close to zero, and a smaller peak at a density
much larger than $c_{0}$ are observed deep inside the aggregate phases.

Visual inspection of movies of the disordered states (see Supplemental
Material sup ) reveals that the life times of the density clusters deep inside
the aggregate phases clearly exceed the simulations times. In the rough phase,
however, no persistent clusters of preferred sizes are visible in NESS. See
movies MOVIE2 and MOVIE4 for the time-evolution of the phase and density
configurations in the rough phase, and the corresponding movies MOVIE5 and
MOVIE6 for the aggregate phase in Ref. sup . Macroscopic distinctions between
the rough and aggregate phase steady states are clearly visible from these
movies.

A more detailed analysis of the different kind of disordered states (see Fig.
27, Fig. 28 and Fig. 29), changes in the nature of the disordered states
depending upon $c_{s}$ (see Fig. 30 and Fig. 31), and the existence of the
rough phases for $g_{2}=0$ as obtained from our simulations (see Fig. 32 and
Fig. 33), are discussed in details in Appendix O.

We find that the magnitudes and signs of the coefficients in the agent-based
model determine the nature of the disordered states. Depending upon the sign
of these model parameters ($g_{2}$, $c_{s}$, and $c_{s2}$) the nearly phase-
ordered system can get disordered corresponding to the aggregate phases or the
rough phases in our simulations, which are manifested by disordered spin
morphology and the formation of clusters of particles of typical sizes in NESS
(see Fig. 27; Appendix O.1). Further, we also note that it is thus possible to
switch between the ordered and rough phases by tuning the values of the active
coefficients. For instance, it is found that increasing only the magnitude of
$g_{1}$ can results in a transition from ordered to rough phase as noted in
Table 1 (also see Fig. 2). See Table 2 for a general classification of the
steady states, parametrized by the signs of the model parameters.

Nature of ordering as parametrized by the model parameters
---
$g_{1}$ |  $g_{2}$ |  $c_{s}$ |  $c_{s2}$ |  Observed phase
$+ve$ |  $+ve$ |  $+ve$ |  $+ve$ |  Ordered phase
$+ve$ |  $+ve$ |  $+ve$ |  $-ve$ |  Aggregate phase
$+ve$ |  $+ve$ |  $-ve$ |  $+ve$ |  Ordered phase
$+ve$ |  $+ve$ |  $-ve$ |  $-ve$ |  Aggregate phase
$-ve$ |  $+ve$ |  $+ve$ |  $+ve$ |  Ordered phase
$-ve$ |  $+ve$ |  $+ve$ |  $-ve$ |  Aggregate phase
$-ve$ |  $+ve$ |  $-ve$ |  $+ve$ |  Ordered phase
$-ve$ |  $+ve$ |  $-ve$ |  $-ve$ |  Aggregate phase
$+ve$ |  $-ve$ |  $+ve$ |  $+ve$ |  Aggregate phase
$+ve$ |  $-ve$ |  $+ve$ |  $-ve$ |  Ordered phase
$+ve$ |  $-ve$ |  $-ve$ |  $+ve$ |  Aggregate phase
$+ve$ |  $-ve$ |  $-ve$ |  $-ve$ |  Ordered phase
$-ve$ |  $-ve$ |  $+ve$ |  $+ve$ |  Aggregate phase
$-ve$ |  $-ve$ |  $+ve$ |  $-ve$ |  Ordered phase
$-ve$ |  $-ve$ |  $-ve$ |  $+ve$ |  Aggregate phase
$-ve$ |  $-ve$ |  $-ve$ |  $-ve$ |  Ordered phase
Table 2: Nature of order (or lack thereof) in the NESS in our agent-based
model, depending upon the possible signs of the active coefficients. For the
phases written in bold, the absolute values of the model parameters correspond
to the ordered phases of TABLE 1. The existence of the aggregate phase does
not depend on the magnitudes of the parameter values, rather it depends on the
sign of the active coefficients.

In the equilibrium limit of the model, we find the phase variance to display
the $\langle\theta^{2}\rangle\sim\ln L$ behavior at sufficiently low noises,
which is the hallmark of the two-dimensional equilibrium XY model (see Fig.
12; Appendix C).

###  IV.3 Summary of the numerical results

We now summarize our results obtained from the Monte-Carlo simulations of the
agent-based models.

(i) The steady states in the agent-based model can be (a) phase-ordered, or
(b) disordered.

(ii) In the phase-ordered states, the variance of the phase fluctuations grows
with the system size $L$. Depending upon the model parameters that determine
the various active dynamical processes, $\langle\theta^{2}\rangle$ can grow
faster or slower than $\ln L$. Similarly, the standard deviation of the number
fluctuations $\sigma(N_{0})$ can grow with $N_{0}$ faster or slower than
$\sqrt{N_{0}}$. Our numerical results reveal a surprising and hitherto unknown
correspondence between the phase and density fluctuations - the model exhibits
the simultaneous occurrence of SQLRO (WQLRO) with MNF (GNF).

(iii) The disordered states are found to be of two distinct types: (a)
aggregate phase characterized by stable density clusters of distinctive size,
and a probability distribution of the density that is not peaked at $c_{0}$,
the mean density, but is bimodal with one peak close to zero, and another at a
density much larger than $c_{0}$, (b) rough phase with no stable density
clusters, corresponding to a density probability distribution that is peaked
around $c_{0}$, but much broader than those for the phase-ordered states.

##  V Derivation of the Hydrodynamic theory

We have demonstrated the existence of stable phase-ordered states mediated by
a mobile conserved species density in our agent-based model. It should be
remembered that this model is a particular model in the general class of
models with many parameters arising from, e.g., more complex forms for the
phase different-dependent and density dependent local rotations, and phase
difference-dependent particle hopping in the update rules. One may
additionally consider replacing the simple spin alignment rules by other
alignment rules involving for instance slightly longer ranged (longer than
nearest neighbor) attractions. If we include all of these, we can make a
microscopic model with practically as many parameters as we may want. However,
so long as all these new variants of the original update rules are all “short
ranged”, i.e., fall off sufficiently fast as a function of the separation, and
retain all the basic invariances of the original update rules, they should
ultimately show the same long wavelength properties as the original model. A
hydrodynamic theory is a particularly powerful tool to establish that, which
ultimately classify the universality class(es). Such a theory necessarily has
a finite number of relevant parameters for a given class of microscopic models
having arbitrary (and practically unlimited) number of parameters. A
hydrodynamic theory for the agent-based model should also provides firmer
footing to the plethora of remarkable results obtained from the agent-based
simulations, by validating them analytically. To that end, we set up the
hydrodynamic theory for the stable phase-ordered states, and calculate the
scaling properties of fluctuations in those states. We also study the
stability conditions of the phase-ordered states, which are expected to give
us clues about the nature of the disordered states found in the agent-based
model.

###  V.1 Identifying the hydrodynamic variables

As in our agent-based model, we consider an active system consisting of a
collection of mobile particles of fixed total number on a substrate, with each
lattice site containing one XY spin. Thus each site is assigned an “XY spin
vector” of fixed length. In order to formulate hydrodynamics, we must coarse-
grain this particle-based description into a position-dependent two-component
spin vector field, or equivalently a complex scalar field $Z({\bf r},t)$,
characterized by an amplitude and a phase. In principle, the coarse-grain
procedure should lead to fluctuations in both the amplitude and phase of $Z$.
However, since we are considering a nearly phase-ordered collection of XY
spins, and not at any continuous transition between the orientationally
ordered and disordered phases, the spin amplitude fluctuations are massive, or
have relaxation times that remain finite in the long wavelength limit. This is
due to the fact that these amplitude fluctuations, unlike the fluctuations in
the phases, are not Goldstone modes associated with spontaneous breaking of
continuous rotation invariance of the spin space. Thus, in the absence of any
symmetry arguments for the amplitude dynamics to be slow, we expect it to be
generically fast. Hence, in the spirit of hydrodynamics, we ignore the
amplitude fluctuations, and set $Z^{2}=1$ without any loss of generality. In
contrast, invariance under a global rotation of the XY spins (equivalently, a
rotation of the chosen reference state) implies that the broken symmetry phase
fluctuations are slow variables with relaxation time-scales diverging in the
long wavelength limit. Furthermore, since we do not consider birth or death of
the mobile species, their number density must be a conserved quantity,
corresponding to it being a slow variable as well. Furthermore, due to the
friction from the substrate, there is no momentum conservation. Therefore, in
the long wavelength limit, we neglect the amplitude fluctuations and retain
phase $\theta({\bf x},t)$, and number density $c({\bf x},t)$ as the
hydrodynamic variables.

In another different physical realization of the system, we could imagine the
system to have mobile particles of fixed numbers, each carrying an XY spins on
a $2d$ substrate, which has the same hydrodynamic equations for the phase-
ordered states, owing to the identical number of broken symmetry mode and
conserved variable (one each), and have the same symmetries.

###  V.2 Hydrodynamic equations of motion

We now systematically formulate the hydrodynamic equations of motion of the
phase $\theta({\bf x},t)$, and number density $c({\bf x},t)$. For a driven
system as ours, in the absence of any free energy functional to derive the
equation of motions, we must write down the equations of motion by appealing
to the general symmetries and conservation laws. Here, symmetries not only
include the symmetries of the underlying microscopic dynamics, they should
also include the symmetries of the states as well. This in particular means it
depends on the specific symmetries broken in the ordered states. The relevant
symmetries are as follows:

(i) For a collection of phase-ordered XY spins, the relevant broken symmetry
is the rotational invariance in the order parameter space 111This is analogous
to the breaking of the rotational invariance in the “spin” space of the
classical spin models, e.g., the XY or the Heisenberg model, in the low
temperature ferromagnetic phase. This is not to be confused with the
rotational invariance of the physical space, which is assumed to be isotropic
(i.e., unbroken) throughout this article.

(ii) Further, the system should be invariant under a translation and rotation
in the physical space and a constant shift of the phase (equivalently, under
rotation of all the spins, or of the reference state by an arbitrary amount).

These symmetry considerations dictate that the equation of motion for $\theta$
should have the general form

$\partial_{t}\theta=\kappa\nabla^{2}\theta+\frac{\lambda}{2}({\bm{\nabla}}\theta)^{2}+\Omega(c)+f_{\theta},$
(16)

where, $\kappa$ is the stiffness, $\Omega(c)$ is a general function of the
density, and we have retained the lowest order nonlinear terms in $\theta$ and
spatial gradients, in the spirit of hydrodynamics assuming small phase
fluctuations in a nearly phase-ordered state. The $\kappa\nabla^{2}\theta$
represents the standard spin relaxation term of equilibrium origin. The
remaining two terms on the rhs of (16) represent active rotation in response
to the local phase difference and number density. None of these terms can be
obtained from a free energy functional. If we set $c=c_{0}=const.$, i.e.,
consider uniform density, then (16), the equation of motion of $\theta$,
reduces to the well-known Kardar-Parisi-Zhang equation Kardar  _et al._
(1986), after absorbing a constant term $\Omega(c_{0})$ by going to a rotating
frame.

The equation for the number density $c$ should follow a conservation law and
have the form

$\partial_{t}c=D\nabla^{2}c+\lambda_{0}{\bm{\nabla}}\cdot(\tilde{\Omega}(c){\bm{\nabla}}\theta)+f_{c}.$
(17)

Equation (17) corresponds to a particle current

${\bf
J}_{c}=-D{\bm{\nabla}}c-\lambda_{0}\tilde{\Omega}(c){\bm{\nabla}}\theta.$ (18)

Here, $D$ is the diffusivity of the mobile species, and $\lambda_{0}$ is a
coupling constant that can be positive or negative. Again we have truncated up
to the lowest order in $\theta$ and spatial gradients; $\tilde{\Omega}(c)$ is
another general functions of $c$; in general $\Omega(c)$ and
$\tilde{\Omega}(c)$ are two different independent nonlinear functions of $c$.
The $\lambda_{0}\tilde{\Omega}(c){\bm{\nabla}}\theta$-term describes particle
current in response to a local phase difference; it may enhance or suppress
the diffusive current $D{\bm{\nabla}}c$. The two signs of $\lambda_{0}$ for a
given $\tilde{\Omega}(c)$ are reminiscent of contractile and extensile active
matters Marchetti  _et al._ (2013). Indeed, we shall see below that the two
possible signs of $\tilde{\lambda}_{0}$ corresponds to two different physical
states of the system. Gaussian noises $f_{\theta}$ and $f_{c}$ are white and
conserved noises, respectively, which are added phenomenologically to model
the inherent underlying stochastic dynamics. These are the nonequilibrium
analogs of thermal noises in equilibrium models. Their variances are
independent of the damping coefficients in a driven system as ours in the
absence of any Fluctuations-Dissipation Theorem (FDT) Chaikin and Lubensky
(1995). We assume the variances to have the forms as given in (4) and (5)
above.

Equations (16) and (17) can be obtained from a complex Ginzburg-Landau
equation (CGLE) coupled with a density field $c$; see Appendix F. The
advection by a flow field is irrelevant in the RG sense is discussed in
Appendix E.

Now write $c({\bf x},t)$ as a sum of the mean concentration $c_{0}$ and local
fluctuations around it, i.e., $c({\bf x},t)=c_{0}+\delta c({\bf x},t)$. We
expand the functions $\Omega(c)$ and $\tilde{\Omega}(c)$ about $c_{0}$. With
this, the hydrodynamic equations retaining up to the lowest order
nonlinearities take the form

$\displaystyle\frac{\partial\theta}{\partial
t}=\kappa\nabla^{2}\theta+\Omega_{1}\delta
c+\frac{\lambda}{2}({\bm{\nabla}}\theta)^{2}+\Omega_{2}(\delta
c)^{2}+f_{\theta},$ (19a) $\displaystyle\frac{\partial\delta c}{\partial
t}=\lambda_{0}\tilde{\Omega}_{0}\nabla^{2}\theta+D\nabla^{2}\delta
c+\lambda_{0}\tilde{\Omega}_{1}{\bm{\nabla}}\cdot(\delta
c{\bm{\nabla}}\theta)+f_{c}.$ (19b)

Here, $\Omega_{1}\equiv\partial\Omega/\partial
c|_{c=c_{0}},\,\Omega_{2}\equiv\partial^{2}\Omega/\partial^{2}c|_{c=c_{0}},\,\tilde{\Omega}_{0}\equiv\tilde{\Omega}|_{c=c_{0}},\,\tilde{\Omega}_{1}\equiv\partial\tilde{\Omega}/\partial
c|_{c=c_{0}}$. Without any loss of generality, we set $\Omega_{1}>0$;
$\tilde{\Omega}_{0}$ can be positive or negative. Using the fact that the
phase $\theta$, an angle, is dimensionless, and $c$ has the dimension of
surface (areal) density, we can fix the dimensions of the various model
parameters in (19a) and (19b). Parameters $\kappa$ and $D$ have the dimensions
of diffusivity, $\lambda_{0}\tilde{\Omega}_{0}$ has the dimensions of
frequency, $\lambda,\,\Omega_{1}$ and $\lambda_{0}\tilde{\Omega}_{1}$ have the
dimensions of diffusivity, and $\Omega_{2}$ has the dimension of
hyperdiffusivity. Lastly, a term $\Omega(c_{0})$ in the hydrodynamic equation
for $\theta$ has been removed by going to a rotating frame: We shift phase
$\theta({\bf r},t)$ by $\Omega(c_{0})$, $\theta({\bf
r},t)\rightarrow\theta({\bf r},t)+\Omega(c_{0})t$, which removes the term
$\Omega(c_{0})$ from (16), that gets generated upon expanding $\Omega(c)$
about $c=c_{0}$. While we have derived the hydrodynamic equations (19a) and
(19b) from symmetry principles and conservation laws, these equations could
also be derived from a general Smoluchowski-like Equation for a basic kinetic
model corresponding to this system, and then by calculating the relevant
moments of the distributions. We should however keep in mind that the
orientational order in the present problem is in the order parameter space,
distinct from the physical space; see Ref. sai (2015); Saintillan and Shelley
(2015) for detailed discussions.

We note that the noises $f_{\theta}$ and $f_{c}$ in Eqs. (16) and (17), or
(19a) and (19b) have their origins in the stochasticity of the agent-based
model. For example, noise $f_{\theta}$ is generated upon coarse-graining the
microscopic phase update rules, which are the stochastic updates of the phase
$\theta$, characterized by $p_{\theta}$ and the randomness in the sites to be
updated in a given Monte-Carlo step, and the “Vicsek”-like noise $\xi$ 222 We
make a technical point here for interested readers. It is expected that such a
coarse-graining procedure should also generated multiplicative noises which
will multiply the different terms in the hydrodynamic equations (as well as
the subleading terms neglected in the hydrodynamic limit). However, since the
microscopic random processes are all short range, $f_{\theta}$ too is short
range. This then gives that all the multiplicative noise terms are irrelevant
in the scaling sense. This consideration leaves $f_{\theta}$ as the only
relevant noise in (19a). Similarly, the noise $f_{c}$ is generated upon
coarse-graining the microscopic density update rules characterized by $p_{c}$,
which include random hopping to a neighboring site with a lower density and
the randomness in selecting the target site in the phase difference-dependent
current. Again, coarse-graining should generate multiplicative noises, all of
which are irrelevant in the scaling sense for reasons similar to the above.
This leaves $f_{c}$ as the only relevant noise in (19b). We have set the noise
crosscorrelations to zero, keeping only the auto-correlations (4) and (5) non-
zero. The model at hand being a driven model, the noise variances and the
damping coefficients are mutually independent. This allows us to choose the
noise variances independently. Whether or not noise crosscorrelations exist
depends upon whether the microscopic stochastic processes affecting the phase
and density updates are crosscorrelated or not. For simplicity, we have chosen
them not to be crosscorrelated.

The hydrodynamic theory for the active XY model in contact with a conserved
species on a substrate developed here is the hydrodynamic theory for active
superfluids on a substrate as well. Hydrodynamic theory for an equilibrium
superfluid has four hydrodynamic variables, viz., phase, density, momentum
density and energy density (equivalently, entropy density). In its driven,
nonequilibrium analogue on a substrate the last two are no longer hydrodynamic
variables, leaving only the phase and density as the two hydrodynamic
variables in the problem. This justifies the applicability of Eqs. (19a) and
(19b) for an active superfluid on a substrate.

For a growing surface or a moving biomembrane, $\theta({\bf x},t)$ is the
local height measured in the Monge gauge with respect to an arbitrary base
plane. Due to the arbitrariness of the location of the base plane, the
dynamics of the surface fluctuations is unaffected by the position of the base
plane. This in turn implies that the dynamics must be invariant under a
constant shift $\theta({\bf x},t)\rightarrow\theta({\bf x},t)+const.$, which
is the same symmetry as the phase. As in a given dimension all systems having
the same symmetry are described by the same hydrodynamic theory, Eqs. (19a)
and (19b) are the hydrodynamic equations for a fluctuating surface with a
conserved species on it. Indeed, due to these symmetry reasons these equations
also serve as the hydrodynamic theory of an active fluid membrane with a
finite tension but without momentum conservation that contains an active
species in it Cagnetta  _et al._ (2020), as mentioned earlier. Nonetheless,
there is a crucial difference between the active XY model problem and the
surface dynamics problem. In the former case, $\theta({\bf x},t)$ is a compact
variable, being confined between 0 and $2\pi$ Chaikin and Lubensky (1995). In
contrast, for a surface, $\theta({\bf x},t)$ is non-compact or unbounded. This
has the consequence that there can be topological defects in the spin model,
which have no analog in the height version of the problem. However, since we
are considering smoothly varying systems, this distinction between the spin
and the height versions of the problem does not arise.

Notice that the density fluctuations cease to become hydrodynamic variables in
two distinct limits: (i) In the immobile limit with conserved but strictly
uniform density, which is achieved when the spins are rigidly grafted on the
lattice points with no other conserved species present, and (ii) particle non-
conservation via birth-death events Altman  _et al._ (2015). In each of these
cases, the phase $\theta$ satisfies the pure KPZ equation Altman  _et al._
(2015).

Although we motivate this study by principally considering active XY spins
grafted on the lattice rigidly together with a diffusing conserved density
trapped near the surface, the hydrodynamic equations (19a) and (19b) apply to
any system having the same symmetries, i.e., invariance under a translation
and rotation in the physical space and a constant shift of the phase
(equivalently, under rotation of all the spins, or of the reference state by
an arbitrary amount). We may consider, for instance, a collection of
diffusively mobile XY spins. This can happen when the mobile particles
themselves carry a XY spin. In the phase-ordered states of this version of the
model that has the same symmetry as above, we again have one broken symmetry
mode (the phase), and one conserved density (of the mobile spins). The
hydrodynamic equations for this model are same as (19a) and (19b). Consider
further an active system consisting of immobile XY spins grafted on lattice
points and an incompressible asymmetric binary fluid on a substrate. An
asymmetric binary fluid is described by a scalar order parameter $\phi$ that
does not have the Ising symmetry, having a non-zero mean. Due to the
conservation of the molecules of the two components that make up the binary
fluid, $\phi$ follows a conserved dynamics. Further, overall
incompressiblility implies that the local mean density is a constant. Lastly,
the presence of a substrate ensures non-conservation of momentum, and the
latter drops out of the hydrodynamics theory. This leaves us with the phase
$\theta$ and the order parameter $\phi$ as the only two hydrodynamic
variables. The rotational symmetry of the phase, absence of the Ising symmetry
of $\phi$ and the conservation of the latter ensures that the coupled dynamics
of $\theta$ and $\phi$ are given by Eqs. (19a) and (19b); see Refs. Safran
(2003); Cates and Tjhung (2018). Interestingly, one-dimensional versions of
Eqs. (19a) and (19b) map onto the hydrodynamic equations for a sedimenting
one-dimensional crystal Lahiri and Ramaswamy (1997) due to symmetry reasons.
Equations (19a) and (19b) also apply to a chiral active hexatic on a substrate
again due to symmetry reasons Maitra  _et al._ (2020).

##  VI Scaling and correlations in the linear theory

It is instructive to first focus on the linearized hydrodynamics [set
$\lambda=\Omega_{2}=\tilde{\Omega}_{1}=0$ in Eqs. (19a) and (19b)], and study
their properties. It is convenient to scale $\delta c$ by $\alpha$, which
brings the linearized equations to the forms

$\displaystyle\frac{\partial\theta}{\partial
t}=\kappa\nabla^{2}\theta+\frac{\Omega_{1}}{\alpha}\delta c+f_{\theta},$ (20a)
$\displaystyle\frac{\partial\delta c}{\partial
t}=\alpha\lambda_{0}\tilde{\Omega}_{0}\nabla^{2}\theta+D\nabla^{2}\delta
c+f_{c},$ (20b)

where we have absorbed a factor of $\alpha$ in $f_{c}$. We now fix $\alpha$ by
setting

$\Omega_{1}/\alpha=\alpha|\lambda_{0}\tilde{\Omega}_{0}|=V_{0}$ (21)

without any loss of generality. Let us now analyze the mode structure of the
coupled linear equations (20a) and (20b). By Fourier transforming (20a) and
(20b) in space and time we get the dispersion relation

$\omega=\frac{i}{2}\left[-(\kappa+D)k^{2}\pm\sqrt{(\kappa-D)^{2}k^{4}-4\lambda_{0}\tilde{\Omega}_{0}\Omega_{1}k^{2}}\right].$
(22)

There are two distinct cases, delineated by the sign of the product
$\lambda_{0}\tilde{\Omega}_{0}\Omega_{1}$. We first consider
$\lambda_{0}\tilde{\Omega}_{0}\Omega_{1}>0$. We are interested in the long
wavelength limit. We then get the dispersion relation for Fourier frequency
$\omega$ as

$\omega=\pm V_{0}k-i\Gamma k^{2},$ (23)

for sufficiently small wavevector $k$, a wavevector regime for which
$V_{0}^{2}\gg\tilde{\Gamma}^{2}k^{2}/4,\;\Gamma=(\kappa+D)/2$ and further
$\tilde{\Gamma}=\kappa-D$. This corresponds to a linearly stable system with
underdamped propagating waves with speed $V_{0}$; see Fig. 7.

Figure 7: Schematic phase diagram in the $\tilde{\Gamma}-|V_{0}|$ plane for a
given wavevector $k$. Linearly stable and unstable regions are marked. The
magnitude of the slopes of the lines demarcating the two is $k/2$. We focus on
the linearly stable region (see text).

In other words, the underdamped propagating waves will be visible provided the
linear system size $L$ is large enough to ensure $V_{0}\gg\pi\tilde{\Gamma}/L$
for a fixed $\tilde{\Gamma}$; $k_{min}=2\pi/L$ is the smallest wavevector. See
Fig. 8.

Figure 8: Phase diagram in the $1/L-|V_{0}|$ plane. Linearly stable and
unstable regions are marked. Slope of the boundary line between the linearly
stable and unstable regions is $\pi\tilde{\Gamma}$. We focus on the linearly
stable region in the hydrodynamic theory (see text).

We have found that the sign of the product of the coefficients of the linear
cross coupling terms in (20a) and (20b) controls the linear stability. This
can be understood as follows. Assume $\Omega_{1}/\alpha>0$ without any loss of
generality. Now consider a local region with excess particles, i.e., $\delta
c>0$, surrounded by regions of deficit particles ($\delta c<0$). This will
“speed up” the local rotation of the spins, $\partial\theta/\partial t>0$ in
the local patch where $\delta c>0$, surrounded by regions with
$\partial\theta/\partial t<0$, implying slowing down there. This results into
the orientation of the spins in the excess particle region to be “locally
ahead” of the neighboring spins, giving a negative Laplacian of the phase
locally. This will enhance $\delta c$ locally (i.e., more excess particles
locally) if $\lambda_{0}\tilde{\Omega}_{0}\Omega_{1}<0$, which in turn
increases $\partial\theta/\partial t$ further, eventually creating a run away
process leading to the linear instability discussed above. On the other hand,
if $\lambda_{0}\tilde{\Omega}_{0}\Omega_{1}>0$, an initial local positive
$\delta c$ will reduce it via the linear $\nabla^{2}\theta$-term in (20b),
ruling out any linear instability. See Ref. Lahiri and Ramaswamy (1997);
Lahiri  _et al._ (2000) for a similar mechanism in a one-dimensional
sedimenting lattice. In the surface growth version of the model, this linear
instability implies that local excess or deficit of particles (i.e., $\delta
c\neq 0$) enhances the local curvature; see Ref. Banerjee and Basu (2018) for
a similar mechanism of instability in an active membrane.

The correlation functions of $\theta({\bf k},\omega)$ and $\delta c({\bf
k},\omega)$ in the linearly stable case in the linear theory can be calculated
from (20a) and (20b) in a straightforward manner:

$\displaystyle\langle|\theta({\bf
k},\omega)|^{2}\rangle=\frac{2D_{\theta}(\omega^{2}+D^{2}k^{4})+2V_{0}^{2}D_{c}k^{2}}{(V_{0}^{2}k^{2}-\omega^{2})^{2}+\omega^{2}k^{4}(\kappa+D)^{2}},$
(24) $\displaystyle\langle|\delta c({\bf
k},\omega)|^{2}\rangle=\frac{2D_{\theta}V_{0}^{2}k^{4}+2D_{c}k^{2}(\omega^{2}+\kappa^{2}k^{4})}{(V_{0}^{2}k^{2}-\omega^{2})^{2}+\omega^{2}k^{4}(\kappa+D)^{2}},$
(25) $\displaystyle\langle\theta(-{\bf k},-\omega)\delta c({\bf
k},\omega)\rangle\equiv\langle\delta c(-{\bf k},-\omega)\theta({\bf
k},\omega)\rangle^{*}=$
$\displaystyle~{}~{}~{}~{}\frac{2V_{0}k^{2}\left[-D_{\theta}(-i\omega+Dk^{2})+D_{c}(i\omega+\kappa
k^{2})\right]}{(V_{0}^{2}k^{2}-\omega^{2})^{2}+\omega^{2}k^{4}(\kappa+D)^{2}}.$
(26)

Unsurprisingly, both $\langle|\theta({\bf k},\omega)|^{2}\rangle$ and
$\langle|\delta c({\bf k},\omega)|^{2}\rangle$ are real and positive definite,
but the cross-correlation function $\langle\theta(-{\bf k},-\omega)\delta
c({\bf k},\omega)\rangle$ is in general complex. All the correlators
$\langle|\theta({\bf k},\omega)|^{2}\rangle,\,\langle|\delta c({\bf
k},\omega)|^{2}\rangle$ and $\langle\theta(-{\bf k},-\omega)\delta c({\bf
k},\omega)\rangle$ are peaked at $\omega=\pm V_{0}k$ in the long wavevlength
limit, with the width of each peak scaling with $k$ as $k^{2}$ (again in the
long wavelength limit); the latter implies a dynamic exponent $z=2$ in the
linear theory. This further implies that in the limit $k\rightarrow 0$, the
ratio of the distance between the two peaks ($=2V_{0}k$) and their widths
($=2\Gamma k^{2}$) diverges: define Quality factor $\varphi=V_{0}/(\Gamma k)$
as the ratio of the two, which diverges for $k\rightarrow 0$. Although Eqs.
(19a) and (19b) have their individual damping coefficients, respectively,
$\kappa$ and $D$, in the linear coupled system of (20a) and (20b) it is their
sum $\kappa+D\equiv 2\Gamma$ that acts as the effective damping. Thus for
stability reasons, we must have $\Gamma>0$; individual positivity of $\kappa$
and $D$ are formally not necessary for stability, although since $\kappa$ and
$D$ are physically identified with the microscopic spin stiffness and particle
diffusivity, respectively, they are expected to be individually positive.

We now calculate the equal-time correlation functions of $\theta({\bf
x},t),\,\delta c({\bf x},t)$. We start from the definition

$C^{0}_{\theta\theta}(k)\equiv\langle|\theta({\bf
k},t)|^{2}\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\langle|\theta({\bf
k},\omega)|^{2}\rangle,$ (27)

where a superscript “0” refers to this being a linear theory result. We now
use (24) above. By applying the residue theorem, we get for the equal-time
phase autocorrelation function in the Fourier space in the linear theory

$C^{0}_{\theta\theta}(k)\approx\frac{D_{\theta}+D_{c}}{2\Gamma k^{2}},$ (28)

see Appendix G.2. Correlator (28) implies

$\Delta_{\theta}^{0}\equiv\langle\theta^{2}({\bf
x},t)\rangle=\frac{D_{\theta}+D_{c}}{4\pi\Gamma}\ln\left(\frac{L}{a_{0}}\right),$
(29)

corresponding to QLRO. As mentioned earlier, (29) holds for the variance of
the surface fluctuations in the EW equation Barabási  _et al._ (1995) as well.

It is useful to compare the linear theory result in (28) with the
corresponding result for the equilibrium XY model in its QLRO phase at
temperature $T$ as given in (1) above: We note that the $k$-dependence is
unchanged from its equilibrium counterpart. In fact (1) and (29) have the same
$\ln L$-dependence. Direct comparison with (29) allows us to define an
effective temperature $T_{eff}^{\theta}$, an effective spring constant
$\kappa_{eff}$ and a dimensionless or reduced temperature ${\cal T}_{\theta}$
for the phase fluctuations by

$\displaystyle T_{eff}^{\theta}$ $\displaystyle\equiv D_{\theta}+D_{c},$ (30a)
$\displaystyle\kappa_{eff}$ $\displaystyle\equiv 2\Gamma,$ (30b)
$\displaystyle{\cal T}_{\theta}$
$\displaystyle\equiv\frac{T^{\theta}_{eff}}{\kappa_{eff}}=\frac{D_{\theta}+D_{c}}{2\Gamma}.$
(30c)

Equation (29) allows us to calculate the Debye-Waller factor (DWF) $W$, which
for an equilibrium system tells us about the depression of the order parameter
due to thermal fluctuations in the ordered phase; here, it would reveal the
corresponding depression of the order parameter due to the noises, which
generalize the role of temperature in a driven nonequilibrium system. This is
defined as Chaikin and Lubensky (1995)

${\cal R}e\langle\exp(i\theta)\rangle=\exp(-\langle\theta^{2}({\bf
x},t)\rangle/2)\equiv\exp(-W)$ (31)

Thus we find

$W\equiv\frac{1}{2}\langle\theta^{2}({\bf
x},t)\rangle=\frac{D_{\theta}+D_{c}}{8\pi\Gamma}\ln\left(\frac{L}{a_{0}}\right),$
(32)

such that the order parameter

${\cal
R}e\langle\exp(i\theta)\rangle=\exp[-W]=\left(\frac{L}{a_{0}}\right)^{-\eta_{0}},$
(33)

where

$\eta_{0}=\frac{D_{\theta}+D_{c}}{8\pi\Gamma}$ (34)

is model-dependent. Thus, the order parameter decays algebraically with $L$ as
$L^{-\eta_{0}}$, eventually vanishing in thermodynamic limit
$L\rightarrow\infty$. Defining the spin variable as a complex function $Z({\bf
r},t)\equiv Z_{0}\exp\left[i\theta({\bf x},t)\right]$ and noting that deep in
the phase-ordered state, amplitude $Z_{0}\approx const.$, we can calculate the
spin-spin equal-time correlation function in the real space in the linear
theory:

$\displaystyle C_{ZZ}^{0}({r})$ $\displaystyle\equiv$
$\displaystyle\langle\cos[\theta({\bf x},t)-\theta({\bf
x^{\prime}},t)]\rangle$ (35) $\displaystyle=$ $\displaystyle{\cal
R}e\langle\exp[i\\{\theta({\bf x},t)-\theta({\bf x^{\prime}},t)\\}]\rangle$
$\displaystyle=$ $\displaystyle\exp[-\langle\\{\theta({\bf x},t)-\theta({\bf
x^{\prime}},t)\\}^{2}\rangle/2],$

where $r=|\bf x-x^{\prime}|$. As shown in Appendix I,

$C_{\theta\theta}^{0}(r)\equiv\langle\left[\theta({\bf x},t)-\theta({\bf
x^{\prime}},t)\right]^{2}\rangle\approx\frac{D_{\theta}+D_{c}}{2\pi\Gamma}\ln(r/a_{0})$
(36)

for large $r$. This gives

$C_{ZZ}^{0}(r)\approx\exp\left[-\frac{D_{\theta}+D_{c}}{4\pi\Gamma}\ln(r/a_{0})\right]\equiv\left(r/a_{0}\right)^{-\tilde{\gamma}_{0}}$
(37)

for large $r$, where $\tilde{\gamma}_{0}=(D_{\theta}+D_{c})/(4\pi\Gamma)$.
Thus the equal-time spin correlation function decays algebraically with
distance, eventually vanishing at infinite distance. The above forms of $W$,
$C^{0}_{\theta\theta}(r)$ and $C_{ZZ}^{0}(r)$ are the hallmarks of QLRO, and
are identical in form with their counterparts for the $2d$ XY model in its low
temperature QLRO phase.

We now focus on the number density fluctuations and calculate the equal-time
density correlator defined as:

$\langle|\delta c({\bf
k},t)|^{2}\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\langle|\delta
c({\bf k},\omega)|^{2}\rangle.$ (38)

We now use (25) above and integrate over $\omega$, see Appendix G.2. This
gives in the linear theory

$C_{cc}^{0}(k)\equiv\langle|\delta c({\bf
k},t)|^{2}\rangle\approx\frac{D_{\theta}+D_{c}}{2\Gamma}$ (39)

in the long wavelength limit, that is independent of $k$. It is useful to
compare (39) with the equilibrium density auto correlation function of a non-
critical system, which at temperature $T$ has the form

$C_{cc}(k)_{eq}\equiv\langle|\delta c({\bf
k},t)|^{2}\rangle|_{eq}=\frac{T}{\chi},$ (40)

which is independent of $k$. Here, $\chi$ is the effective equilibrium
susceptibility that is finite in the long wavelength limit away from any
critical point. Comparing (39) with (40), we find

$\displaystyle T_{eff}^{c}\equiv D_{\theta}+D_{c},$ (41a)
$\displaystyle\chi\equiv 2\Gamma,$ (41b)

for the effective temperature $T_{eff}^{c}$ and effective susceptibility
$\chi$, respectively, in the long wavelength limit. The equal-time real space
density auto correlation function

$C_{cc}^{0}(r)\equiv\langle\delta c({\bf x},t)\delta c({\bf
x^{\prime}},t)\rangle$ (42)

is the inverse Fourier transform of $C_{cc}^{0}(k)$. For the convenience of
evaluating this inverse Fourier transform, we now include a hyperdiffusion
term $-\zeta\nabla^{4}\delta c$ in (19b) to get

$\frac{\partial\delta c}{\partial
t}=\lambda_{0}\tilde{\Omega}_{0}\nabla^{2}\theta+D\nabla^{2}\delta
c-\zeta\nabla^{4}\delta
c+\lambda_{0}\tilde{\Omega}_{1}{\bm{\nabla}}\cdot(\delta
c{\bm{\nabla}}\theta)+f_{c}.$ (43)

This gives

$\displaystyle\langle\delta c({\bf x},t)\delta c({\bf
x^{\prime}},t)\rangle\approx\frac{D_{c}+D_{\theta}}{2\Gamma}\int_{k=0}^{\infty}\frac{d^{2}k}{(2\pi)^{2}}\frac{\exp(i{\bf
k\cdot r})}{1+\zeta^{2}k^{2}}$
$\displaystyle\approx\frac{D_{c}+D_{\theta}}{4\Gamma}\exp(-r/\zeta)\frac{1}{\sqrt{r}},$
(44)

for $r\gg\zeta$; a convergence factor $1/(1+\zeta^{2}k^{2})$ has been
introduced to ensure convergence of the integral, where we have used the
asymptotic form of the Bessel function for large $r$. Here, $\zeta$ is a
length scale that sets the scale of the short range interaction; for
$r\gg\zeta$, $\langle\delta c({\bf x},t)\delta c({\bf x^{\prime}},t)\rangle$
decreases rapidly.

Lastly, we consider the cross-correlation function $\langle\delta
c({\bf-k},-\omega)\theta({\bf k},\omega)\rangle$. We start from (26) above and
use the definition

$\langle\delta c({\bf-k},t)\theta({\bf
k},t)\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\langle\delta
c({\bf-k},-\omega)\theta({\bf k},\omega)\rangle.$ (45)

Staying within the linear theory and again integrate over $\omega$, we find

$C_{\times}^{0}(k)\equiv\langle\delta c({\bf-k},t)\theta({\bf
k},t)\rangle=\frac{D_{c}\kappa-D_{\theta}D}{2\Gamma V_{0}};$ (46)

see Appendix G.2 for a detailed derivation. This immediately gives

$C_{\times}^{0}(r)\equiv\langle\delta c({\bf x},t)\theta({\bf
x^{\prime}},t)\rangle\approx\frac{D_{c}\kappa-D_{\theta}D}{4\pi\Gamma
V_{0}}\exp(-r/\zeta)\frac{1}{\sqrt{r}}$ (47)

for $r\gg\zeta$. We define a dimensionless ratio $\varpi$ that will be useful
in our subsequent discussions. It is defined in the linear theory as

$\varpi=\frac{C_{cc}^{0}(k)}{C_{\times}^{0}(k)}$ (48)

that is independent of $k$ and is an ${\cal O}(1)$ number when all the model
parameters are themselves ${\cal O}(1)$ numbers.

At this stage, it is interesting to note that we can extract the same
effective temperature from both $\langle\theta({\bf x},t)^{2}\rangle$ and
$\langle(\delta c({\bf x},t))^{2}\rangle$: $T_{eff}^{\theta}=T_{eff}^{c}\equiv
T_{eff}$, apparently implying an effective equilibrium-like behavior for the
linearly coupled system of $\theta$ and $\delta c$ at a temperature $T_{eff}$.
Nonetheless, a non-zero cross-correlation function $\langle\delta c({\bf
x},t)\theta({\bf x^{\prime}},t)\rangle$ breaks the condition of thermal
equilibrium, or the Fluctuation-Dissipation Theorem (FDT) Chaikin and Lubensky
(1995); see Appendix J.

Let us further examine the consequence of (44). Let us consider a finite area
$\tilde{L}\times\tilde{L}$ and study the density fluctuations in it. If we
consider two such small areas of equal size with a linear size $\tilde{L}$,
one centered at ${\bf x}$ and the other at ${\bf x}^{\prime}$, we can
calculate the variance of the number fluctuations in each of the areas and how
that variance scale with the area. Now writing ${\bf
r}=(x-x^{\prime},\,y-y^{\prime})$, we get for the variance of the number
fluctuations $\delta N$ about the mean $N_{0}$ in an area
$\tilde{L}\times\tilde{L}$

$\displaystyle\langle(\delta
N)^{2}\rangle=\int_{0}^{\tilde{L}}dx\int_{0}^{\tilde{L}}dy\int_{0}^{\tilde{L}}dx^{\prime}\int_{0}^{\tilde{L}}dy^{\prime}\langle\delta
c({\bf x},t)\delta c({\bf x^{\prime}},t)\rangle$ (49) $\displaystyle\approx$
$\displaystyle\int_{0}^{\tilde{L}}dx\int_{0}^{\tilde{L}}dy\int_{0}^{\tilde{L}}dx^{\prime}\int_{0}^{\tilde{L}}dy^{\prime}\frac{D_{c}+D_{\theta}}{4\Gamma}\exp(-r/\zeta)\frac{1}{\sqrt{r}}\,d\phi$
$\displaystyle\approx$
$\displaystyle\tilde{L}^{2}\frac{D_{c}+D_{\theta}}{4\Gamma},$

assuming $\tilde{L}\gg\zeta$, where $\phi$ is the polar angle. Thus
$\langle(\delta N)^{2}\rangle$ scales with $\tilde{L}$ as $\tilde{L}^{2}$.
Since the mean number $N_{0}$ in an area $\tilde{L}\times\tilde{L}$ scales
with $\tilde{L}^{2}$, we find that $\langle(\delta N)^{2}\rangle$ scales as
$N_{0}$, giving the standard deviation of the number fluctuations
$\sigma(N_{0})\propto\sqrt{N_{0}}$, corresponding to normal number fluctuation
(NNF), as expected in an equilibrium system with short range interactions away
from any critical points.

On the other hand, for $\lambda_{0}\tilde{\Omega}_{0}\Omega_{1}<0$, the
dispersion relation by using (22) becomes

$\omega=\pm iV_{0}k,$ (50)

to the lowest order in $k$. Thus one of the modes grows with a growth rate
that scales with $k$ corresponding to linear instability. These instabilities
are static.

We thus conclude that in the linearized theory, a nearly phase-ordered state,
characterized by $\langle\theta({\bf x},t)^{2}\rangle$ and $\langle(\delta
c({\bf x}.t)^{2}\rangle$, is either statistically indistinguishable from the
$2d$ XY model in its QLRO phase and ordinary diffusing particles in a lattice-
gas system, respectively, with both being at temperature $T_{eff}$, in so far
as the phase fluctuations and the density fluctuations are concerned, or the
system is linearly unstable implying absence of any ordering. In the next
Section, we focus on the linearly stable phase-ordered states and ask whether
these are also nonlinearly stable or not. If these phase-ordered states are
indeed nonlinearly stable, then we ask: what are their universal scaling
properties? We address this issue systematically below.

##  VII Nonlinear effects

To know whether nonlinear effects preserve or destabilize the linearly stable
uniform states in $2d$, we start from the hydrodynamic equations that contain
the lowest order nonlinear terms in fields and gradients:

$\displaystyle\partial_{t}\theta=\kappa\nabla^{2}\theta+\frac{\lambda}{2}({\bm{\nabla}}\theta)^{2}+V_{0}\delta
c+\Omega_{2}(\delta c)^{2}+f_{\theta},$ (51a) $\displaystyle\partial_{t}\delta
c=D\nabla^{2}\delta
C+V_{0}\nabla^{2}\theta+\lambda_{1}{\bm{\nabla}}\cdot(\delta
c{\bm{\nabla}}\theta)+f_{c},$ (51b)

where, we have retained only the terms lowest order in fields and gradients;
$\lambda_{1}\equiv\lambda_{0}\tilde{\Omega}_{1}$. Noises $f_{\theta}$ and
$f_{c}$ are still assumed to be zero-mean Gaussian-distributed with variances
given by (4) and (5); we ignore any cross-correlations between $f_{\theta}$
and $f_{c}$ for simplicity. We are interested in the scaling properties in the
long wavelength limit, by which we specifically mean wavevector $k$ satisfying
$k\ll\frac{V_{0}}{\kappa-D}$, such that there are underdamped propagating
waves in the wavevector regime of interest. Identifying ${\bm{\nabla}}\theta$
with a “velocity” ${\bf v}_{s}\equiv{\bm{\nabla}}\theta$, the resulting
equations have structural similarities with the coupled Burgers (BMHD)
equations; see Refs. Basu  _et al._ (1999); Basu and Frey (2004, 2009) and
Appendix L.

Equations (51a) and (51b) are invariant under ${\bf x}\rightarrow-{\bf
x},\,\theta({\bf x},t)\rightarrow\theta(-{\bf x},t),\,\delta c({\bf
x},t)\rightarrow\delta c(-{\bf x},t)$. Further, if $\lambda=\lambda_{1}$, the
model equations (51a) and (51b) admit invariance under a pseudo-Galilean
transformation Forster  _et al._ (1977); Barabási  _et al._ (1995); Basu and
Frey (2004):

${\bf x}\rightarrow{\bf x}+{\bf
U}_{0}t,\,{\bm{\nabla}}\theta\rightarrow{\bm{\nabla}}\theta+{\bf
U}_{0},\,\frac{\partial}{\partial t}\rightarrow\frac{\partial}{\partial
t}+\lambda{\bf U}_{0}\cdot{\bm{\nabla}},$ (52)

where $U_{0}$ is a constant. In general however there is no physical
requirement of imposing the Galilean invariance on the model equations (51a)
and (51b). Hence, in this model there are no restrictions on the magnitudes
and signs of $\lambda,\,\lambda_{1}$; in fact, they can be of the same or
opposite signs. In the limit of immobile particles, i.e., if we set
$V_{0}=\Omega_{2}=0$, Eq. (51a) reduces to the well-known Kardar-Parisi-Zhang
(KPZ) equation Kardar  _et al._ (1986); Barabási  _et al._ (1995). It is known
that at $2d$ the KPZ equation only has a rough phase, which for the phase
variable $\theta({\bf x},t)$ means destruction of the phase-ordered state. We
will now study whether the additional nonlinear terms and the propagating
modes can alter this picture, and make stable order possible.

###  VII.1 Renormalization group analysis

The presence of the nonlinear terms in (51a) and (51b) preclude any exact
enumeration of the correlation functions. Instead, perturbative approaches are
needed. Naïve perturbation theory, however, produces diverging corrections to
the model parameters, similar to the $2d$ KPZ equation. These divergences are
systematically handled within the dynamic RG framework Forster  _et al._
(1977); Hohenberg and Halperin (1977); Bausch  _et al._ (1976); Barabási  _et
al._ (1995), which is conveniently implemented by using a path integral
description, equivalent to and constructed from Eqs. (51a) and (51b), in terms
of the fields $\theta({\bf x},t)$ and $\delta c({\bf x},t)$, and their dynamic
conjugate fields, respectively, $\hat{\theta}({\bf x},t)$ and $\hat{c}({\bf
x},t)$ Bausch  _et al._ (1976), after averaging over the noise distributions
with variances given by (4) and (5).

The momentum shell RG procedure consists of integrating over the short
wavelength Fourier modes of $\theta({\bf r},t)$, $\delta c({\bf r},t)$,
$\hat{\theta}({\bf r},t)$ and $\hat{c}({\bf r},t)$, followed by rescaling of
lengths and times Forster  _et al._ (1977); Hohenberg and Halperin (1977);
Barabási  _et al._ (1995). In particular, we follow the usual approach of
initially restricting the wavevectors to be within a bounded circular
Brillouin zone: $|{\bf k}|<\Lambda$. However, the precise value of the upper
cutoff $\Lambda$ has no effect on our final results. The fields $\theta({\bf
r},t)$, $\delta c({\bf r},t)$, $\hat{\theta}({\bf r},t)$ and $\hat{c}^{(}{\bf
r},t)$ are separated into the high and low wave vector parts $\theta({\bf
r},t)=\theta^{<}({\bf r},t)+\theta^{>}({\bf r},t)$, $\delta c({\bf
r},t)=\delta c^{<}({\bf r},t)+\delta c^{>}({\bf r},t)$, $\hat{\theta}({\bf
r},t)=\hat{\theta}^{<}({\bf r},t)+\hat{\theta}^{>}({\bf r},t)$, and
$\hat{c}({\bf r},t)=\hat{c}^{<}({\bf r},t)+\hat{c}^{>}({\bf r},t)$ where
$\theta^{>}({\bf r},t)$, $\delta c^{>}({\bf r},t)$, $\hat{\theta}^{>}({\bf
r},t)$ and $\hat{c}^{>}({\bf r},t)$ have support in the large wave vector
(short wavelength) range $\Lambda e^{-l}<|{\bf k}|<\Lambda$, while
$\theta^{<}({\bf r},t)$, $\delta c^{<}({\bf r},t)$, $\hat{\theta}^{<}({\bf
r},t)$ and $\hat{c}^{<}({\bf r},t)$ have support in the small wave vector
(long wavelength) range $|{\bf k}|<e^{-l}\Lambda$; $b\equiv e^{l}>1$. We then
integrate out $\theta^{>}({\bf r},t)$, $\delta c^{>}({\bf r},t)$,
$\hat{\theta}^{>}({\bf r},t)$ and $\hat{c}^{>}({\bf r},t)$ perturbatively in
the anhamornic coupling $\lambda,\,\Omega_{2}$ and $\lambda_{1}$; as usual,
this resulting perturbation theory of $\theta^{<}({\bf r},t)$, $\delta
c^{<}({\bf r},t)$, $\hat{\theta}^{<}({\bf r},t)$ and $\hat{c}^{<}({\bf r},t)$
can be represented by Feynman graphs, with the order of perturbation theory
reflected by the number of loops in the graphs we consider. The Feynman
vertices representing the anharmonic couplings are illustrated in Fig. 15 in
Appendix G.1. After this perturbative step, we rescale lengths, with ${\bf
r}\rightarrow{\bf r}^{\prime}e^{l}$, which restores the UV cutoff back to
$\Lambda$, together with rescaling of time $t\rightarrow t^{\prime}e^{zl}$
(equivalently in the Fourier space, the momentum and frequency are rescaled as
${\bf k}\rightarrow{\bf k}^{\prime}/b$ and
$\omega\rightarrow\omega^{\prime}/b^{z}$, respectively), where $z$ is dynamic
exponent. This is then followed by rescaling the long wave length part of the
fields that we define in the Fourier space for calculational convenience. We
fix the scale factors by demanding that certain harmonic part of the action
functional (161) in Appendix G.1 do not scale. In $d$-space dimensions, we
have

$\displaystyle\theta({\bf k},\omega)=b^{\chi_{\theta}}\theta^{\prime}({\bf
k^{\prime}},\omega^{\prime}),\,\delta c({\bf k},\omega)=b^{\chi_{c}}\delta
c^{\prime}({\bf k^{\prime}},\omega^{\prime}),$ $\displaystyle\hat{\theta}({\bf
k},\omega)=b^{\hat{\chi}_{\theta}}\hat{\theta}^{\prime}({\bf
k^{\prime}},\omega^{\prime}),\,\hat{c}({\bf
k},\omega)=b^{\hat{\chi}_{c}}\hat{c}^{\prime}({\bf
k^{\prime}},\omega^{\prime}).$ (53)

This gives the exponents
$\chi_{\theta}=4,\chi_{c}=3,\hat{\chi}_{\theta}=2,\hat{\chi}_{c}=3$ at $2d$;
see also Appendix H.2.

We restrict ourselves here to a one-loop perturbative RG calculation. See
Appendix H for the detailed calculations.

Before we proceed further, let us quickly revisit the dynamic RG calculation
for the $2d$ KPZ equation Tang  _et al._ (1990); Barabási  _et al._ (1995);
Täuber (2014); Kardar  _et al._ (1986). Our model reduces to the KPZ equation
if $\delta c=0=V_{0}$. One defines an effective dimensionless coupling
constant
$\tilde{g}=\frac{\lambda^{2}D_{\theta}}{4\kappa^{3}}\frac{k_{d}}{(2\pi)^{d}}$
in $d$ dimensions; $k_{d}$ is the area of a unit hypersphere in $d$ dimension.
The RG flow equation for the coupling $\tilde{g}$ at $2d$ (which is the lower
critical dimension of the KPZ equation) reads

$\frac{d\tilde{g}}{dl}=\tilde{g}^{2}.$ (54)

This can be easily integrated to give a scale-dependent $\tilde{g}(l)$

$\tilde{g}(l)=\frac{\tilde{g}(0)}{1-\tilde{g}_{0}l},$ (55)

where $\tilde{g}(l=0)$ is a constant of integration, and is the value of
$\tilde{g}$ at the smallest scale or its unrenormalized value; $l$ is the
logarithm of a length-scale, and is the “RG time”. Clearly at the RG fixed
point, $d\tilde{g}(l)/dl=0$, giving $\tilde{g}=0$ is the only fixed point,
which is an unstable fixed point. In fact, as $l$ approaches a finite
nonuniversal value determined by the microscopic or unrenormalized parameters,
$\tilde{g}(l)$ diverges at a finite $l$, signaling breakdown of the
perturbation theory. In terms of a height field that satisfies the $2d$ KPZ
equation, this is interpreted as the existence of a rough phase, inaccessible
to perturbative calculations, or for the phase variable $\theta$, the non-
existence of any phase ordered state at $2d$ in the thermodynamic limit.
Notice that the flow of the coupling $\tilde{g}$ to infinity (as a result of
$d\tilde{g}/dl>0$) is crucial to this breakdown of perturbation theory and the
consequent lack of any phase-ordered states. Flow equation (54) shows that at
$2d$, coupling $\tilde{g}$ is marginally relevant (consistent with the lack of
a smooth phase in the $2d$ KPZ equation, or the lack of an ordered state at
$2d$ here); $\tilde{g}=0$, which is the only fixed point of (54), is
marginally unstable.

We calculate the relevant one-loop corrections to all the model parameters of
our model equations (51a) and (51b); see Appendix H.1 for details. It turns
out that only $\kappa,D,D_{\theta},D_{c}$ have relevant fluctuation-
corrections at the one-loop order due to the nonlinear terms. These
corrections are $\log$-divergent at $2d$. However, there are no diverging
fluctuation corrections to the wave speed $V_{0}$ or the coupling constants
$\lambda,\,\lambda_{1}$ and $\Omega_{2}$ at the one-loop order.

Rescaling space, time and the fields in the standard way as discussed above
(see also Appendix H.1), we get the RG recursion relations for the various
model parameters at the one-loop order. We define two dimensionless parameters

$\mu_{1}=\frac{\Omega_{2}}{\lambda},\;\,\mu_{2}=\frac{\lambda_{1}}{\lambda}.$
(56)

We also define

$\displaystyle\overline{D}=D_{\theta}+D_{c},\,\tilde{D}=D_{\theta}-D_{c}.$
(57)

We find

$\displaystyle\frac{d\tilde{D}}{dl}=A(\mu_{1},\mu_{2})g\tilde{D},$ (58a)
$\displaystyle\frac{d\overline{D}}{dl}=G(\mu_{1},\mu_{2})g\overline{D},$ (58b)

where

$\displaystyle A(\mu_{1},\mu_{2})=\mu_{1}^{2}-\mu_{2}^{2}/2+\mu_{1}+1/4,$
(59a) $\displaystyle
G(\mu_{1},\mu_{2})=\mu_{1}^{2}+\mu_{2}^{2}/2+\mu_{1}+1/4.$ (59b)

Both $A$ and $G$ are dimensionless ratios; $A$ can be of any sign, but $G$ is
non-negative by construction. We further define the dimensionless coupling
constant

$g=\frac{\lambda^{2}\overline{D}}{8\Gamma^{3}}\frac{k_{d}}{(2\pi)^{d}}$ (60)

in $d$-dimensions, that is the analogue of the dimensionless coupling constant
$\tilde{g}$ in the KPZ problem.

By using the rescaling factors defined in (53), we find the scaling of $g$ as

$g^{\prime}=b^{2-d}g,$ (61)

showing that at $2d$, $g$ does not scale at all under naïve rescaling of
space, time and the long wavelength parts of the fields. This means $2d$ is
the critical dimension, just as $2d$ is the critical dimension of the KPZ
equation. We further define

${\cal
N}\equiv\frac{\tilde{D}}{\overline{D}}=\frac{D_{\theta}-D_{c}}{D_{\theta}+D_{c}},$
(62)

which is a dimensionless ratio. We find

$\frac{d{\cal N}}{dl}=g\left[A(\mu_{1},\mu_{2})-G(\mu_{1},\mu_{2}){\cal
N}\right].$ (63)

Thus, the flow of $\cal N$ is independent of all other parameters (apart from
$\mu_{1}$ and $\mu_{2}$). At the RG fixed point

${\cal N}^{*}=\frac{A(\mu_{1},\mu_{2})}{G(\mu_{1},\mu_{2})},$ (64)

which is linearly stable. Here and below a “superscript ∗” represents a RG
fixed point value. The RG recursion equations for $\Gamma,\tilde{\Gamma}$ read
(we set ${\cal N}={\cal N}^{*}$ in the long RG time limit)

$\displaystyle\frac{d\Gamma}{dl}=\Gamma\left[g\left\\{B(\mu_{1},\mu_{2})-\frac{C(\mu_{1},\mu_{2}){\cal
N}^{*}}{4}+\left(\frac{\mu_{2}^{2}-1}{8}+\frac{C(\mu_{1},\mu_{2}){\cal
N}^{*}}{4}\right)\frac{\tilde{\Gamma}}{\Gamma}\right\\}+\tilde{g_{1}}\frac{C(\mu_{1},\mu_{2})}{2}\frac{\tilde{\Gamma}}{\Gamma}\right],$
(65a)
$\displaystyle\frac{d\tilde{\Gamma}}{dl}=\Gamma\left[g\left\\{F(\mu_{1},\mu_{2})-\frac{H(\mu_{1},\mu_{2}){\cal
N}^{*}}{4}+\left(\frac{4\mu_{1}\mu_{2}-\mu_{2}^{2}-1}{8}+\frac{H(\mu_{1},\mu_{2}){\cal
N}^{*}}{4}\right)\frac{\tilde{\Gamma}}{\Gamma}\right\\}+\tilde{g_{1}}\frac{H(\mu_{1},\mu_{2})}{2}\frac{\tilde{\Gamma}}{\Gamma}\right],$
(65b)

where

$\displaystyle\tilde{g_{1}}=\frac{G(\mu_{1},\mu_{2})-A(\mu_{1},\mu_{2})}{2G(\mu_{1},\mu_{2})}g.$
(66)

$\displaystyle
B(\mu_{1},\mu_{2})=\frac{7\mu_{1}\mu_{2}}{2}+\frac{9\mu_{2}^{2}}{8}+\mu_{1}+\mu_{2}+\frac{1}{8},$
(67a) $\displaystyle
C(\mu_{1},\mu_{2})=\mu_{1}\mu_{2}-\mu_{2}^{2}/2+\mu_{2}/2-\mu_{1},$ (67b)
$\displaystyle
H(\mu_{1},\mu_{2})=\mu_{2}^{2}/2-\mu_{1}\mu_{2}-\mu_{1}+\mu_{2}/2,$ (67c)
$\displaystyle
F(\mu_{1},\mu_{2})=-\mu_{1}\mu_{2}-\frac{9\mu_{2}^{2}}{8}+\mu_{1}+\mu_{2}+\frac{1}{8}.$
(67d)

Parameters $A,G,B,C,H$ and $F$ depend only on $\mu_{1}(l=0)$ and
$\mu_{2}(l=0)$, and hence are all marginal as well as dimensionless.

Since the nonlinear coupling constants $\lambda,\,\Omega_{2}$ and
$\lambda_{1}$ do not receive any relevant fluctuation corrections themselves
at the one-loop order (see Appendix H.1), the dimensionless parameters
$\mu_{1}$ and $\mu_{2}$ too do not receive any diverging fluctuation
corrections at the one-loop order, i.e., $\mu_{1}$ and $\mu_{2}$ are constants
along the RG flow. This implies

$\frac{d\mu_{1}}{dl}=0=\frac{d\mu_{2}}{dl}.$ (68)

Thus, if the other renormalized model parameters and the scaling exponents
depend upon $\mu_{1}$ and $\mu_{2}$, they will actually depend upon
$\mu_{1}(l=0)$ and $\mu_{2}(l=0)$, i.e., the initial or unrenormalized values
of $\mu_{1}$ and $\mu_{2}$, which are the tuning parameters in this theory as
we see below.

Furthermore, we obtain

$\displaystyle\frac{d}{dl}\left(\frac{\tilde{\Gamma}}{\Gamma}\right)$
$\displaystyle=$ $\displaystyle
g\left[F-\frac{AH}{4G}+\frac{\tilde{\Gamma}}{\Gamma}\left\\{-B+\frac{AC}{4G}+\frac{H}{4}+\frac{4\mu_{1}\mu_{2}-\mu_{2}^{2}-1}{8}\right\\}-\left(\frac{\tilde{\Gamma}}{\Gamma}\right)^{2}\left(\frac{\mu_{2}^{2}-1}{8}+\frac{C}{4}\right)\right]$
(69) $\displaystyle=$ $\displaystyle
g\left[a+m\left(\frac{\tilde{\Gamma}}{\Gamma}\right)+\overline{C}\left(\frac{\tilde{\Gamma}}{\Gamma}\right)^{2}\right].$

where

$\displaystyle a=F-AH/(4G),$ (70a)
$\displaystyle\overline{C}=-C/4-(\mu_{2}^{2}-1)/8,$ (70b) $\displaystyle
m=-B+AC/(4G)+H/4$ $\displaystyle+(4\mu_{1}\mu_{2}-\mu_{2}^{2}-1)/8$ (70c)

are yet another set of dimensionless constants. The ratio
$\tilde{\Gamma}/\Gamma$ is yet another dimensionless number. In the spirit of
hydrodynamics and scaling, we assume scaling solutions

$\Gamma(l)=\Gamma(0)l^{\eta_{2}},\;\tilde{\Gamma}(l)=\tilde{\Gamma}(0)l^{\eta_{2}},$
(71)

in the long wavelength limit in a putative stable ordered phase, which satisfy
the flow equations (65a) and (65b) for some yet unknown scaling exponent
$\eta_{2}$. Thence,
$\psi\equiv\tilde{\Gamma}(l)/\Gamma(l)=\tilde{\Gamma}_{0}/\Gamma_{0}$ is a
constant at the RG fixed point. By using (69), we find that $\psi$ follows the
equation

$\psi=\frac{1}{2\overline{C}}\left[-m\pm\sqrt{m^{2}-4a\overline{C}}\right],$
(72)

For a physically meaningful solution, $\psi$ must be real and also be in the
range between -1 and 1, since $|\tilde{\Gamma}|\leq\Gamma$ by definition 333
We have assumed the same scaling exponents $\eta_{2}$ for both $\Gamma(l)$ and
$\tilde{\Gamma}(l)$. If we instead assume different scaling exponents given by
$\Gamma(l)=\Gamma_{0}l^{\eta_{a}}$ and
$\tilde{\Gamma}(l)=\tilde{\Gamma}_{0}l^{\eta_{b}}$, then if
$\eta_{a}>\eta_{b}$, then $\psi=0$, else if $\eta_{a}<\eta_{b}$, then $\psi$
diverges, which is unphysical. We note that the condition for $\psi$ to be
real is $m^{2}>4a\overline{C}$, together with a choice for $\psi$ in eq (72)
to have stable values is
$\frac{1}{2\overline{C}}\left[-m-\sqrt{m^{2}-4a\overline{C}}\right]$ ; see
Fig. 9 for regions in the $\mu_{1}-\mu_{2}$ plane where ${\cal I}m\,\psi=0$.

This is to be supplemented by the RG flow equation for $g$ itself, which
reads:

$\frac{dg}{dl}=g^{2}\left[G-3\left(B-\frac{AC}{4G}\right)+3\psi\overline{C}\right]=-\Delta_{1}g^{2},$
(73)

where

$\Delta_{1}\equiv-\left[G-3\left(B-\frac{AC}{4G}\right)+3\psi\overline{C}\right].$
(74)

If $\Delta_{1}>0$, $dg/dl<0$, $g(l)$ decreases as the “RG time” $l$ increases,
paving the way for possible phase ordered states. Assuming positive
$\Delta_{1}$, solving (73) we get

$g(l)=\frac{g_{0}}{g_{0}\Delta_{1}l+1}\approx\frac{1}{\Delta_{1}l}$ (75)

for large RG time $l$. Here, $g_{0}\equiv g(l=0)$ is the “initial” or
unrenormalized value of $g$. Here and below, we use a suffix “0” to represent
quantities evaluated at the “initial RG time” ($l=0$), or the corresponding
unrenormalized, microscopic values. In contrast, if $\Delta_{1}<0$, then the
solution of (73) gives

$g(l)=\frac{g_{0}}{-g_{0}|\Delta_{1}|l+1}.$ (76)

Evidently, as $l\rightarrow 1/\Delta_{1}g_{0}$, a non-universal value
controlled by the unrenormalized model parameters, $g(l)\rightarrow\infty$,
signaling breakdown of the perturbation theory, which is reminiscent of the
flow of $\tilde{g}(l)$ in the $2d$ KPZ problem. This in turn implies loss of
the phase-ordered state. We do not discuss this further here, and below we
focus on $\Delta_{1}>0$. Therefore, with $\Delta_{1}>0$ in (73), coupling $g$
is marginally irrelevant, and the fixed point $g^{*}=0$ is marginally stable,
whereas for $\Delta_{1}<0$ coupling $g$ is marginally relevant and $g^{*}=0$
is a marginally unstable fixed point; see Fig. 10 for regions in the
$\mu_{1}-\mu_{2}$ plane with $\Delta_{1}>0$. However, $\Delta_{1}>0$ is not
the necessary and sufficient condition for the existence of phase-ordered
states, as we discuss below. Notice that in principle, one could define
several other dimensionless coupling constants by using $\Omega_{2}$ and
$\lambda_{1}$ in place of $\lambda$ in the definition of $g$ in (60), all of
which will be proportional to $g$, the proportionality constants being
combinations of $\mu_{1},\,\mu_{2}$. Thus, their RG flow equations can be
obtained straightforwardly from the RG flow equation (73) of $g$, using the
fact that $\mu_{1}$ and $\mu_{2}$ are marginal. We however choose to work in
terms of $g,\,\mu_{1},\,\mu_{2}$ for reasons of convenience and simplicity of
explanations; see also Ref. Haldar and Basu (2022) for similar technical
issues.

We thus conclude that we must have $\Delta_{1}>0$ together with a real
$-1\leq\psi\leq 1$ for the existence of stable uniform phase-ordered steady
states, and for which the above scaling solutions hold good 444 What is the
nature of the steady state when $\psi$ is complex? From (65a) and (65b), both
$\Gamma(l)$ and $\tilde{\Gamma}(l)$ are complex when $\psi$ is complex. This
means (23) will now have a higher order in $k$ imaginary term of the form
$i\tilde{\Gamma}^{2}k^{3}/V_{0}$, which triggers linear instability at
intermediate wavevectors. Of course, these instabilities would be cutoff by
even higher order linear damping terms of equilibrium origin, neglected in the
hydrodynamic approach. These can potentially lead to patterned steady states..

Now substitute (75) in (58b) to get

$\frac{d\overline{D}}{dl}=\frac{G(\mu_{1},\mu_{2})}{\Delta_{1}(\mu_{1},\mu_{2})l}\overline{D},$
(77)

solving which we find

$\overline{D}(l)=\overline{D}(0)l^{\eta_{1}},$ (78)

with

$\eta_{1}\equiv\frac{G(\mu_{1},\mu_{2})}{\Delta_{1}(\mu_{1},\mu_{2})}>0,$ (79)

since $G(\mu_{1},\mu_{2})>0$ by construction and
$\Delta_{1}(\mu_{1},\mu_{2})>0$ for the stability. Thus
$\overline{D}(l)/\overline{D}(0)$ necessarily diverges in the entire stable
region of the phase space. Similarly, substituting for ${\cal N}^{*}$ and
$\psi=\tilde{\Gamma}(l)/\Gamma(l)$ in (65a) we obtain

$\displaystyle\frac{d\Gamma}{dl}$ $\displaystyle=$
$\displaystyle\Gamma\frac{1}{\Delta_{1}l}\left[B-\frac{AC}{4G}-\overline{C}\psi\right].$
(80)

Solving Eq. (80) we get

$\Gamma(l)=\Gamma(0)l^{\eta_{2}},$ (81)

with

$\displaystyle\eta_{2}$
$\displaystyle\equiv\frac{1}{\Delta_{1}}\left[B-\frac{AC}{4G}-\overline{C}\psi\right]=\frac{1}{3}\left[1+\frac{G}{\Delta_{1}}\right]$
(82) $\displaystyle=\frac{1}{3}(1+\eta_{1}).$

Equation (82) shows that $\eta_{2}>0$. This implies that
$\Gamma(l)\gg\Gamma(0)$ in the renormalized theory in the long wavelength
limit. This further means that the dynamics in the long wavelength limit of
the renormalized theory is faster (albeit logarithmically) than that in the
linear or noninteracting theory; see also Eq. (95) below in this context. As
mentioned above, $\overline{D}(0)$ and $\overline{\Gamma}(0)$ are to be
identified with the respective unrenormalized or bare parameters
$\overline{D}$ and $\Gamma$.

Clearly, the scaling exponents $\eta_{1},\eta_{2}$ that characterize the
stable, ordered phase are parametrized by the real values of ${\cal N}^{*}$
and $\psi$, together with the constraint $\Delta_{1}>0$, which in turn are
functions of $\mu_{1},\,\mu_{2}$.

Figure 9: Regions with $-1\leq\psi\leq 1$ in the $\mu_{1}-\mu_{2}$ plane:
Light colored regions have bounded values of $\psi$. In the remaining region,
$\psi$ has unphysical solutions. Figure 10: Regions with $\Delta_{1}>0$ in the
$\mu_{1}-\mu_{2}$ plane: Light colored regions have $\Delta_{1}>0$; in the
remaining region $\Delta_{1}<0$.

As we have argued above, in the spirit of scaling solutions

$\displaystyle\tilde{D}(l)=\tilde{D}(0)l^{\eta_{1}(\mu_{1},\mu_{2})},$ (83a)
$\displaystyle\tilde{\Gamma}(l)=\tilde{\Gamma}(0)l^{\eta_{2}(\mu_{1},\mu_{2})};$
(83b)

here, $\tilde{D}(0)\equiv(D_{\theta}(l=0)-D_{c}(l=0))$. The above results for
scale-dependent $\overline{D}(l)$ and $\Gamma(l)$ imply wavevector-dependent
renormalized parameters

$\displaystyle\overline{D}(k)\equiv
D_{\theta}(k)+D_{c}(k)=\overline{D}(0)\left[\ln\left(\frac{\Lambda}{k}\right)\right]^{\eta_{1}},$
(84a)
$\displaystyle\Gamma(k)\equiv\left[\kappa(k)+D(k)\right]/2=\Gamma(0)\left[\ln\left(\frac{\Lambda}{k}\right)\right]^{\eta_{2}}.$
(84b)

Assuming finite ${\cal N}^{*}$ and $\psi$, we further find that

$\tilde{D}(k)\approx\tilde{D}_{0}\left[\ln\left(\frac{\Lambda}{k}\right)\right]^{\eta_{1}},\,\tilde{\Gamma}(k)\approx\tilde{\Gamma}_{0}\left[\ln\left(\frac{\Lambda}{k}\right)\right]^{\eta_{2}},$
(85)

in the limit of small wavevector.

As an explicit example, let us work out the case with ${\cal N}^{*}=0$ and
$\psi=0$ (these conditions imply $D_{\theta}^{*}=D_{c}^{*}$ and
$\kappa^{*}=D^{*}$ at the RG fixed point). These conditions are satisfied at
isolated points in the $\mu_{1}-\mu_{2}$ plane. These are calculated from the
simultaneous solutions of ${\cal N}^{*}=0$ and $\psi=0$. The condition ${\cal
N}^{*}=0$ is satisfied by setting $A(\mu_{1},\mu_{2})=0$. On the other hand in
eq. (69), if

$F(\mu_{1},\mu_{2})=0,\;{\rm
and}\;-B+H/4+(4\mu_{1}\mu_{2}-\mu_{2}^{2}-1)/8=m<0,$ (86)

then the solution $\psi=0$ gives a stable RG fixed point. These conditions are
used to find

$\mu_{1}=\frac{\sqrt{2}-1}{2},\mu_{2}=1.$ (87)

We now use the RG flow equations (58b) and (65a) to obtain

$\displaystyle\frac{d\overline{D}}{dl}=\frac{2\sqrt{2}}{27-2\sqrt{2}}\frac{\overline{D}}{l},$
(88a)
$\displaystyle\frac{d\Gamma}{dl}=\frac{9}{27-2\sqrt{2}}\frac{\Gamma}{l}.$
(88b)

Solving, we obtain

$\eta_{1}\approx 0.117,\eta_{2}\approx 0.372$ (89)

in this special case. We next find the scaling exponents given any general
${\cal N}^{*}\neq 0$ and any real $|\psi|\leq 1$, maintaining $\Delta_{1}>0$.
We present the results in a tabular form; see Table 4.

The existence of continuously varying scaling exponents found above depends
crucially on the nonrenormalization of the coupling constants, and hence
$\mu_{1},\,\mu_{2}$ at the one-loop order. In the absence of any symmetries,
there is no guarantee that $\mu_{1},\mu_{2}$ should remain unrenormalized at
higher loop orders. We now argue that these actually do not matter. This is so
because they will involve higher powers of $g$, and will, hence, vanish like
$1/l$ to a power greater than 1. Therefore, their integrals over $l$ from zero
to infinity will be finite, so they will not change the anomalous behavior of
$D$ and $\kappa$. Similarly, they cannot make any divergent contribution to
$\mu_{1}$ and $\mu_{2}$. Therefore, our results, even though they were only
derived to one loop order, are, in fact, asymptotically exact. This
automatically implies that the continuous variation of the scaling exponents,
making them nonuniversal, is also asymptotically exact Toner .

Wavevector-dependent $\overline{D}(k)$ and $\Gamma(k)$ may be used to
calculate the renormalized correlation functions in the stable regimes of the
theory, which we discuss in the next Section VII.2.

###  VII.2 Scaling of the correlation functions and order

Since in the renormalized theory, coupling constant $g(k)\rightarrow 0$ as
$k\rightarrow 0$, the renormalized theory is effectively linear, albeit with
renormalized model parameter. This consideration gives the renormalized
correlation functions straight forwardly.

We start by constructing the renormalized correlation functions in the Fourier
space as a function of $\bf k$ and $\omega$. We get

$\displaystyle\langle|\theta({\bf