# Power-law correlation in the homogeneous disordered state of anisotropically
self-propelled systems

Kyosuke Adachi  RIKEN Interdisciplinary Theoretical and Mathematical Sciences
Program (iTHEMS), 2-1 Hirosawa, Wako 351-0198, Japan  Nonequilibrium Physics
of Living Matter RIKEN Hakubi Research Team, RIKEN Center for Biosystems
Dynamics Research (BDR), 2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe 650-0047,
Japan  Hiroyoshi Nakano  Institute for Solid State Physics, University of
Tokyo, 5-1-5, Kashiwanoha, Kashiwa 277-8581, Japan

###### Abstract

Self-propelled particles display unique collective phenomena, due to the
intrinsic coupling of density and polarity. For instance, the giant number
fluctuation appears in the orientationally ordered state, and the motility-
induced phase separation appears in systems with repulsion. Effects of strong
noise typically lead to a homogeneous disordered state, in which the coupling
of density and polarity can still play a significant role. Here, we study
universal properties of the homogeneous disordered state in two-dimensional
systems with uniaxially anisotropic self-propulsion. Using hydrodynamic
arguments, we propose that the density correlation and polarity correlation
generically exhibit power-law decay with distinct exponents ($-2$ and $-4$,
respectively) through the coupling of density and polarity. Simulations of
self-propelled lattice gas models indeed show the predicted power-law
correlations, regardless of whether the interaction type is repulsion or
alignment. Further, by mapping the model to a two-component boson system and
employing non-Hermitian perturbation theory, we obtain the analytical
expression for the structure factors, the Fourier transform of the correlation
functions. This reveals that even the first order of the interaction strength
induces the power-law correlations.

##  I Introduction

Active matter, a crowd of self-propelled entities, displays a variety of
collective behaviors [1, 2, 3]. Using minimal particle models of self-
propelled elements [4] and hydrodynamic arguments [5, 6], extensive studies
have been performed to elucidate the fundamental mechanisms and universal
properties of collective phenomena in active matter systems [7]. An intrinsic
feature of active matter is the interplay of density and polarity, which
emerges from the polar self-propelled motion of each particle. The coupling of
density and polarity produces distinctive fluctuation properties and
spontaneous inhomogeneity, such as the giant number fluctuation in the
orientationally ordered state [7], motility-induced phase separation (MIPS)
[8], and (micro)phase separation in flocking [9, 10, 11] or active nematics
[12]. Naturally, the density-polarity coupling is expected to play a
significant role even in the homogeneous disordered state, which appears when
noise predominates over the effects of interactions.

According to the studies on externally driven many-body systems [13, 14, 15,
16], violation of the detailed balance induces the long-range density
correlation with power-law decay, as long as spatial anisotropy exists in
dynamics. Recently, based on an analogy between externally driven force and
anisotropic self-propulsion, we have proposed example models of active matter
that show the same type of long-range density correlation. Specifically, we
have focused on the homogeneous disordered state of uniaxially self-propelled
particles with repulsive interactions [17, 18]. Our results raise fundamental
questions: (i) whether the power-law correlation is universal and independent
of the interaction type (e.g., repulsion or alignment), and (ii) whether the
anisotropic self-propulsion can lead to unique properties that have no
counterparts in externally driven systems. Broadly, the effects of spatial
anisotropy on collective behaviors of active matter have been studied for the
flocking transition [19, 20, 21], active nematics [22], and MIPS [23, 24].

Stochastic many-body systems, including active matter models, are generically
mapped to and can be analyzed as non-Hermitian quantum systems [25, 26, 27]. A
prototypical example is the correspondence between the asymmetric simple
exclusion process and the XXZ quantum spin chain [28, 29, 30, 31]. More
recently, extensive studies on non-Hermitian systems [32] have facilitated
mutual interactions between the classical and quantum domains, as exemplified
by topologically protected edge modes [33, 34]. Given such developments,
mapping active matter models to quantum systems not only imports analytical
tools from quantum theory but also potentially offers valuable insights into
non-Hermitian physics.

In this paper, we study universal properties of the homogeneous disordered
state in active matter models with uniaxial self-propulsion, such as variants
of the active lattice gas model for MIPS [35, 36, 37, 38, 39, 40] and the
active Ising model for flocking [9, 10]. From hydrodynamic arguments, we
propose that the density correlation and polarity correlation generically
exhibit power-law decay with distinct exponents ($-2$ and $-4$, respectively)
through the coupling of density and polarity inherent to active matter.
Performing simulations of lattice gas models with repulsive or aligning
interaction, we numerically confirm that the predicted power-law correlations
appear regardless of the interaction type. For these models, employing the
mapping to a quantum system and non-Hermitian perturbation theory, we
analytically find that even the first order of the interaction strength
induces the power-law correlations.

##  II Hydrodynamic argument

We consider a homogeneous disordered state of self-propelled particles in two
dimensions, with the direction of self-propulsion restricted along the $x$
axis. We assume that the particles flip their direction and move by diffusion
or self-propulsion with interactions such as repulsion or alignment. The
macroscopic dynamics is expected to have the same form regardless of the
interaction type and may be described 111Equation (1) is a simpler version of
Eq. (2) in Ref. [17] [or Eq. (S-21) in Ref. [40]], which is derived for a
repulsive active lattice gas model. by a linear fluctuating hydrodynamic
equation:

$\partial_{t}\rho_{\pm}=\mp
v\partial_{x}\rho_{\pm}+\bar{D}\bm{\nabla}^{2}\rho_{\pm}+\eta+\bar{\gamma}(\rho_{\mp}-\rho_{\pm})\pm\xi,$
(1)

where $v,\bar{D},\bar{\gamma}>0$. Here, $\rho_{\pm}(\bm{r},t)$ is the coarse-
grained density field for particles that are self-propelled in the $\pm x$
direction. The first term on the right-hand side represents advection induced
by self-propulsion with velocity $\pm v$. The second and third terms describe
stochastic diffusion in a standard way as used in the so-called Model B for
dynamics of conserved density [42, 43]. $\eta$ is a conserved Gaussian noise
with $\braket{\eta(\bm{r},t)}=0$,
$\braket{\eta(\bm{r},t)\eta(\bm{r}^{\prime},t^{\prime})}=-2\Delta\bm{\nabla}^{2}\delta(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime})$,
and $\Delta>0$. The last two terms denote flipping (i.e., rotation) of the
direction of self-propulsion with rate $\bar{\gamma}$. $\xi$ is a non-
conserved Gaussian noise with $\braket{\xi(\bm{r},t)}=0$,
$\braket{\xi(\bm{r},t)\xi(\bm{r}^{\prime},t^{\prime})}=2\Delta^{\prime}\delta(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime})$,
and $\Delta^{\prime}>0$. The parameters $\bar{D}$ and $\bar{\gamma}$ quantify
the strengths of diffusion and flipping, respectively, the values of which
depend on the detail of the interaction. The main conclusion will not change
if we consider more complex noise correlations or further anisotropic terms
(e.g., by replacing $\bar{D}\bm{\nabla}^{2}$ with
$\bar{D}_{x}{\partial_{x}}^{2}+\bar{D}_{y}{\partial_{y}}^{2}$).

Using the total density field
$\rho(\bm{r},t):=\rho_{+}(\bm{r},t)+\rho_{-}(\bm{r},t)$ and the polarization
density field $w(\bm{r},t):=\rho_{+}(\bm{r},t)-\rho_{-}(\bm{r},t)$ [1], we can
rewrite Eq. (1) as

$\displaystyle\partial_{t}\rho$
$\displaystyle=-v\partial_{x}w+\bar{D}\bm{\nabla}^{2}\rho+2\eta$ (2)
$\displaystyle\partial_{t}w$
$\displaystyle=-v\partial_{x}\rho+\bar{D}\bm{\nabla}^{2}w-2\bar{\gamma}w+2\xi.$
(3)

Though we can directly solve these linear equations, we here use the adiabatic
approximation [44] to reach the essential results in a simple way. We notice
that $\rho(\bm{r},t)$ is a slow variable since the total particle number is
conserved, while $w(\bm{r},t)$ is a fast variable with decay rate
$2\bar{\gamma}$. Thus, we use the adiabatic approximation as

$w\simeq-\frac{v}{2\bar{\gamma}}\partial_{x}\rho+\frac{1}{\bar{\gamma}}\xi,$
(4)

where we neglect higher-order gradient terms by considering long-wavelength
fluctuations. This equation indicates the coupling of density and polarity at
the hydrodynamic level. Substituting this into Eq. (2), we obtain

$\displaystyle\partial_{t}\rho$
$\displaystyle\simeq\left[\left(\bar{D}+\frac{v^{2}}{2\bar{\gamma}}\right){\partial_{x}}^{2}+\bar{D}{\partial_{y}}^{2}\right]\rho+\zeta$
(5) $\displaystyle\zeta$
$\displaystyle:=2\eta-\frac{v}{\bar{\gamma}}\partial_{x}\xi$ (6)

Due to the uniaxial self-propulsion with speed $v$, the obtained diffusion
coefficient is anisotropic [i.e., $\bar{D}+v^{2}/(2\bar{\gamma})$ and
$\bar{D}$ in the $x$ and $y$ directions, respectively], and the noise
amplitude is also anisotropic, as seen from
$\braket{\zeta(\bm{r},t)\zeta(\bm{r}^{\prime},t^{\prime})}=-2[(4\Delta+v^{2}\Delta^{\prime}/\bar{\gamma}^{2}){\partial_{x}}^{2}+4\Delta{\partial_{y}}^{2}]\delta(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime})$.

Solving Eq. (5), we can obtain the steady-state density structure factor
$\bar{S}_{d}(\bm{k}):=\Omega^{-1}\lim_{t\to\infty}\braket{}{\delta\rho(\bm{k},t)}{{}^{2}}$,
where $\delta\rho(\bm{k},t):=\int
d^{2}\bm{r}e^{-i\bm{k}\cdot\bm{r}}[\rho(\bm{r},t)-\bar{\rho}_{0}]$ with the
average density $\bar{\rho}_{0}$, and $\Omega$ is the total area of the
system. For $\bm{k}\simeq\bm{0}$, we obtain

$\bar{S}_{d}(\bm{k})\simeq\frac{(4\Delta+v^{2}\Delta^{\prime}/\bar{\gamma}^{2}){k_{x}}^{2}+4\Delta{k_{y}}^{2}}{[\bar{D}+v^{2}/(2\bar{\gamma})]{k_{x}}^{2}+\bar{D}{k_{y}}^{2}}.$
(7)

We notice that $\bar{S}_{d}(\bm{k})$ is generically discontinuous at
$\bm{k}\to\bm{0}$, as characterized by

$\lim_{k_{x}\to 0}\bar{S}_{d}(k_{x},0)-\lim_{k_{y}\to
0}\bar{S}_{d}(0,k_{y})=\frac{v^{2}\Delta^{\prime}}{{\bar{\gamma}}^{2}\bar{D}}\frac{1-2\bar{\gamma}\Delta/(\bar{D}\Delta^{\prime})}{1+v^{2}/(2\bar{\gamma}\bar{D})}\neq
0,$ (8)

except for the special case with $2\bar{\gamma}\Delta=\bar{D}\Delta^{\prime}$
where the fluctuation-dissipation relation happens to hold in Eq. (5).

From the discontinuity of $\bar{S}_{d}(\bm{k})$, the steady-state density
correlation function,
$\bar{C}_{d}(\bm{r}):=\lim_{t\to\infty}\braket{\rho(\bm{0},t)\rho(\bm{r},t)}-{\bar{\rho}_{0}}^{2}=\int
d^{2}\bm{k}e^{i\bm{k}\cdot\bm{r}}\bar{S}_{d}(\bm{k})$, shows power-law decay
for $|\bm{r}|\to\infty$ [15]:

$\begin{split}\bar{C}_{d}(x,0)&\sim x^{-2}\\\ \bar{C}_{d}(0,y)&\sim
y^{-2}.\end{split}$ (9)

Note that the fluctuating hydrodynamic equation similar to Eq. (5) has been
widely used to explain long-range density correlation of the form of Eq. (9)
in anisotropic nonequilibrium systems such as driven lattice gas models [15,
16]. The authors have observed the same type of power-law density correlation
in an active lattice gas model [17] and active Brownian particles [18] with
uniaxial self-propulsion.

We further consider the steady-state polarity correlation function,
$\bar{C}_{p}(\bm{r}):=\lim_{t\to\infty}\braket{w(\bm{0},t)w(\bm{r},t)}$. Since
$\bar{C}_{p}(\bm{r})\sim{\partial_{x}}^{2}\bar{C}_{d}(\bm{r})$ from Eq. (4),
we obtain the power-law polarity correlation for $|\bm{r}|\to\infty$:

$\begin{split}\bar{C}_{p}(x,0)&\sim x^{-4}\\\ \bar{C}_{p}(0,y)&\sim
y^{-4},\end{split}$ (10)

where the exponent is decreased by $2$ compared to the density correlation.
Since we do not assume the interaction type in the above discussion, the
power-law correlations of density [Eq. (9)] and polarity [Eq. (10)] are
expected to appear in many-particle systems with uniaxial self-propulsion
regardless of the detail of the interaction.

##  III Numetical simulations

###  III.1 Lattice gas model with uniaxial self-propulsion

To confirm the prediction [Eqs. (9) and (10)] for uniaxially self-propelled
particles with distinct interaction types, we consider a lattice gas model in
two dimensions with lattice size $L\times L$, total particle number $N$, and
periodic boundary conditions in both directions (Fig. 1). Each particle has
“spin” $s\in\\{+,-\\}$ as an internal state, which corresponds to the
direction of self-propulsion along the $x$ axis. The particle configuration is
specified by $\bm{n}:=\\{n_{i,s}\\}$, where $n_{i,s}$ is the number of
particles with spin $s$ at site $i$. A particle with spin $s$ follows two
kinds of dynamics: (i) the particle hops one site right with rate
$(1+s\varepsilon)D(\bm{n})$, left with rate $(1-s\varepsilon)D(\bm{n})$, and
up or down with rate $D(\bm{n})$ for each, and (ii) the spin flips to $-s$
with rate $\gamma(\bm{n})$. Here, $\varepsilon\in[0,1]$ is a parameter for the
strength of self-propulsion, where $\varepsilon=0$ corresponds to no self-
propulsion. The interactions between particles are represented by the
functional forms of $D(\bm{n})$ and $\gamma(\bm{n})$.

We consider repulsive and aligning interactions within Glauber-type dynamics:

$\begin{split}D(\bm{n})&=D_{0}\frac{2}{1+e^{\Delta
E_{\mathrm{hop}}(\bm{n})}}\\\ \gamma(\bm{n})&=\gamma_{0}\frac{2}{1+e^{\Delta
E_{\mathrm{flip}}(\bm{n})}},\end{split}$ (11)

where $D_{0},\gamma_{0}>0$, and the functional forms of $\Delta
E_{\mathrm{hop}}(\bm{n})$ and $\Delta E_{\mathrm{flip}}(\bm{n})$ specify the
interactions. For hopping from a site $i$ to a neighboring site $j$, we
consider a repulsive interaction:

$\Delta E_{\mathrm{hop}}(\bm{n})=U(n_{j}-n_{i}+1),$ (12)

where $U\geq 0$ is the strength of repulsion, and $n_{i}:=\sum_{s}n_{i,s}$ is
the local density. For spin flipping from $s$ to $-s$, we consider an aligning
interaction:

$\Delta E_{\mathrm{flip}}(\bm{n})=J(sm_{i}-1),$ (13)

where $J\geq 0$ is the strength of alignment, and $m_{i}:=\sum_{s}sn_{i,s}$ is
the local polarity.

Figure 1: Lattice gas model with uniaxial self-propulsion. We consider $N$
particles with spin $+$ (red, right arrow) or $-$ (blue, left arrow) on a
square lattice with size $L\times L$. Each particle stochastically (i) hops to
the neighboring site at a rate favoring the direction of its spin or (ii)
flips its spin. The hopping rate $D(\bm{n})$ and flipping rate
$\gamma(\bm{n})$ depend on the particle configuration $\bm{n}$, representing
the repulsive and aligning interactions between particles, respectively [see
Eq. (11)]. Figure 2: Interaction-dependent phase behavior. (a, b) With
repulsive interaction $U$, (a) the homogeneous disordered steady state appears
for weak $U$, and (b) the system shows MIPS for strong $U$. (c, d) With
aligning interaction $J$, (c) the homogeneous disordered steady state appears
for weak $J$, and (d) the system shows flocking for strong $J$. In each
typical snapshot, the red, blue, and purple dots mean that the site is
occupied by one or more particles with the local polarity $m_{i}$ positive,
negative, and zero, respectively (see also Fig. 1). The arrow in (d) indicates
the direction of the flocking motion. The used parameters are (a, b)
$D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=0.01$, $L=200$, $N=2\times 10^{4}$,
and (a) $U=0.1$ or (b) $U=5$; (c, d) $D_{0}=1$, $\varepsilon=0.8$,
$\gamma_{0}=2$, $L=200$, $N=1.6\times 10^{5}$, and (c) $J=0.1$ or (d) $J=1$.

In Fig. 2, we show the qualitative phase behavior of our model. In the system
with pure repulsion (i.e., $U>0$ and $J=0$), MIPS appears for large $U$ [Fig.
2(b)], as observed in similar models [38, 17, 40]. The cluster configuration
in MIPS reflects the anisotropy of self-propulsion. In the system with pure
alignment (i.e., $U=0$ and $J>0$), flocking appears for large $J$ [Fig. 2(d)],
as observed in the active Ising model [9, 10]. The cluster configuration in
flocking reflects the anisotropy of self-propulsion. When $U$ and $J$ are
small enough, noise overcomes the effects of interactions, and the steady
state is homogeneous and disordered [Figs. 2(a) and (c)]. In the following, we
focus on the collective properties of these homogeneous disordered states.

For the lattice gas model, we define the density correlation function:

$C_{d}(\bm{r}_{j}):=\frac{1}{L^{2}}\sum_{i}\braket{n_{i}n_{i+j}}-{\rho_{0}}^{2},$
(14)

where $\bm{r}_{j}$ is the coordinate of site $j$, $\braket{\cdots}$ is the
steady-state ensemble average, and $\rho_{0}:=N/L^{2}$. Similarly, we define
the polarity correlation function:

$C_{p}(\bm{r}_{j}):=\frac{1}{L^{2}}\sum_{i}\braket{m_{i}m_{i+j}}.$ (15)

We also define the corresponding density and polarity structure factors:

$S_{d}(\bm{k}):=\sum_{j}e^{-i\bm{k}\cdot\bm{r}_{j}}C_{d}(\bm{r}_{j})$ (16)

and

$S_{p}(\bm{k}):=\sum_{j}e^{-i\bm{k}\cdot\bm{r}_{j}}C_{p}(\bm{r}_{j}),$ (17)

where $\bm{k}:=(2\pi n_{x}/L,2\pi n_{y}/L)\in[-\pi,\pi)\times[-\pi,\pi)$ with
$n_{x},n_{y}\in\mathbb{Z}$. Note that $S_{d}(\bm{0})=0$ from the particle
number conservation (i.e., $\sum_{i}n_{i}=N$).

Figure 3: Power-law correlations in the system with repulsion and uniaxial
self-propulsion. (a, b) Heat maps of the structure factors for (a) density
[$S_{d}(\bm{k})$] and (b) polarity [$S_{p}(\bm{k})$]. (c, d) Correlation
functions along the $x$ axis (blue) and $y$ axis (orange) for (c) density
[$C_{d}(\bm{r})$] and (d) polarity [$C_{p}(\bm{r})$]. The inset of (d) is an
expanded view, suggesting that $C_{p}(\bm{r})$ finally decays from the
positive side in both $x$ and $y$ directions. (e, f) Log-log plots of the (e)
density and (f) polarity correlation functions, which correspond to (c) and
(d), respectively. In (e), the two lines with light color suggest the
correlation function after subtracting $-\rho_{0}/L^{2}$, which comes from a
finite-size effect [see Eq. (67)]. The black dashed lines in (e) and (f) are
algebraic functions, indicating that $C_{d}(\bm{r})\sim r^{-2}$ and
$C_{p}(\bm{r})\sim r^{-4}$. The used parameters are $D_{0}=1$,
$\varepsilon=0.8$, $\gamma_{0}=0.5$, $U=0.05$, $L=200$, and $N=1.6\times
10^{6}$. Figure 4: Power-law correlations in the system with alignment and
uniaxial self-propulsion. We plot (a, b) the structure factors and (c-f) the
correlation functions in the same way as shown in Fig. 3. In particular, (e)
and (f) suggest that $C_{d}(\bm{r})\sim r^{-2}$ and $C_{p}(\bm{r})\sim
r^{-4}$. The used parameters are $D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=2$,
and $J=0.025$, $L=200$, and $N=1.6\times 10^{6}$.

###  III.2 Power-law correlations for repulsive interaction

We first consider the system with purely repulsive interaction (i.e., $U>0$
and $J=0$). We set $D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=0.5$, $U=0.05$,
$L=200$, and $N=1.6\times 10^{6}$, where the steady state is homogeneous and
disordered. We selected a high density, $\rho_{0}=N/L^{2}=40$, to clearly see
the asymptotic behavior of the correlation functions with a limited number of
samples. See Appendix A for the simulation method and Appendix B for the
numerical sampling of configurations in the steady state.

In Figs. 3(a) and (b), we show the heat maps of the obtained density structure
factor $S_{d}(\bm{k})$ and polarity structure factor $S_{p}(\bm{k})$,
respectively. $S_{d}(\bm{k})$ is discontinuous at $\bm{k}\to\bm{0}$, as
expected from the hydrodynamic argument in Sec. II. Correspondingly, the
density correlation function $C_{d}(\bm{r})$ exhibits the power-law decay with
exponent $-2$ in the $x$ and $y$ directions [Figs. 3(c) and (e)], which is
consistent with Eq. (9). Lastly, the polarity correlation function
$C_{p}(\bm{r})$ shows the power-law decay with exponent $-4$ [Figs. 3(d) and
(f)] as predicted in Eq. (10).

###  III.3 Power-law correlations for aligning interaction

We next consider the system with purely aligning interaction (i.e., $U=0$ and
$J>0$). We set $D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=2$, and $J=0.025$,
$L=200$, and $N=1.6\times 10^{6}$ (corresponding to $\rho_{0}=40$), where the
steady state is homogeneous and disordered. See Appendix A for the simulation
method and Appendix B for the numerical sampling of configurations in the
steady state.

In Fig. 4, we arrange the figures in the same order as Fig. 3, the case of
repulsive interaction. The observed density structure factor $S_{d}(\bm{k})$
[Fig. 4(a)] and density correlation function $C_{d}(\bm{r})$ [Figs. 4(c) and
(e)] are qualitatively similar to those in the repulsive case [Figs. 3(a),
(c), and (e)]. Particularly, we observe the power-law decay with exponent $-2$
in the $x$ and $y$ directions, as predicted by the hydrodynamic argument [Eq.
(9)]. The polarity structure factor $S_{p}(\bm{k})$ and correlation function
$C_{p}(\bm{r})$ [Figs. 4(b) and (d)] exhibit distinct behaviors from those for
the repulsive case [Figs. 3(b) and (d)]. Specifically, polarity shows positive
correlations in all directions [Fig. 4(d)] due to the aligning interaction,
while the repulsive interaction leads to the negative polarity correlation at
a short distance [Fig. 3(d)]. Despite such difference in the short-distance
property, the long-distance polarity correlation function $C_{p}(\bm{r})$
exhibits the power-law decay with exponent $-4$ in the $x$ and $y$ directions
[Fig. 4(f)], which is also consistent with the hydrodynamic argument [Eq.
(10)].

##  IV Perturbation theory for interactions

To show that very weak interaction, regardless of repulsion or alignment, is
sufficient for the power-law correlations observed in the lattice gas model
(Sec. III), we analytically derive the correlation functions up to the first
order in the interaction strength. For the systematic derivation of the
distribution of particle configurations, we map our model (Fig. 1) to a non-
Hermitian quantum model by the Doi-Peliti method [25, 26, 27] and apply
perturbation theory [45]. With this approach, we can obtain analytical
expressions of the structure factors without taking the hydrodynamic limit.
Readers who are not interested in the detailed derivation can skip to the main
results in Sec. IV.4.

###  IV.1 Doi-Peliti method

The master equation for our model is written as

$\frac{d}{dt}P(\bm{n},t)=\sum_{\bm{n}^{\prime}}W(\bm{n},\bm{n}^{\prime})P(\bm{n}^{\prime},t),$
(18)

where $P(\bm{n},t)$ is the time-dependent distribution of the configuration
$\bm{n}$, and $W(\bm{n},\bm{n}^{\prime})$ is the transition rate from
$\bm{n}^{\prime}$ to $\bm{n}$. The steady-state distribution
$P_{\mathrm{st}}(\bm{n})$ should satisfy

$\sum_{\bm{n}^{\prime}}W(\bm{n},\bm{n}^{\prime})P_{\mathrm{st}}(\bm{n}^{\prime})=0.$
(19)

Regarding $W(\bm{n},\bm{n}^{\prime})$ as a matrix element of $W$, we can think
of $P_{\mathrm{st}}(\bm{n})$ as a component of the eigenvector of $W$
corresponding to the eigenvalue of zero.

We introduce a boson Fock space spanned by Fock bases
$\\{\ket{\bm{n}}\\}_{\bm{n}}$. Each Fock basis $\ket{\bm{n}}$
($=\ket{\\{n_{i,s}\\}}$) is defined using the vacuum state $\ket{0}$ as

$\ket{\bm{n}}:=\prod_{i,s}\left(\hat{a}_{i,s}^{\dagger}\right)^{n_{i,s}}\ket{0},$
(20)

and the corresponding left Fock basis is defined as

${}_{L}\\!\bra{\bm{n}}:=\bra{0}\prod_{i,s}\frac{1}{n_{i,s}!}\left(\hat{a}_{i,s}\right)^{n_{i,s}}.$
(21)

Here, $\hat{a}_{i,s}^{({\dagger})}$ is the annihilation (creation) operator of
a boson with spin $s$ at site $i$, and we assume the commutation relation:
$[\hat{a}_{i,s},\hat{a}_{j,s^{\prime}}^{\dagger}]=\delta_{i,j}\delta_{s,s^{\prime}}$
and
$[\hat{a}_{i,s},\hat{a}_{j,s^{\prime}}]=[\hat{a}_{i,s}^{\dagger},\hat{a}_{j,s^{\prime}}^{\dagger}]=0$.
The number operator of a boson with spin $s$ at site $i$ is given as
$\hat{n}_{i,s}:=\hat{a}_{i,s}^{\dagger}\hat{a}_{i,s}$. Note that
$\ket{\bm{n}}$ and ${}_{L}{\bra{\bm{n}}}$ satisfy biorthogonality:
${}_{L}\\!\braket{\bm{n}}{\bm{n}^{\prime}}=\delta_{\bm{n},\bm{n}^{\prime}}$.

We define the Fourier transformation of the operator $\hat{a}_{j,s}$ for
wavevector $\bm{k}=(2\pi n_{x}/L,2\pi n_{y}/L)$ with
$n_{x},n_{y}\in\mathbb{Z}$ and $k_{x},k_{y}\in[-\pi,\pi)$:

$\hat{a}_{\bm{k},s}:=\frac{1}{L}\sum_{j}e^{-i\bm{k}\cdot\bm{r}_{j}}\hat{a}_{j,s},$
(22)

where $\bm{r}_{j}$ is the coordinate of site $j$. The operators satisfy the
boson commutation relation:
$[\hat{a}_{\bm{k},s},\hat{a}_{\bm{k}^{\prime},s^{\prime}}^{\dagger}]=\delta_{\bm{k},\bm{k}^{\prime}}\delta_{s,s^{\prime}}$
and
$[\hat{a}_{\bm{k},s},\hat{a}_{\bm{k}^{\prime},s^{\prime}}]=[\hat{a}_{\bm{k},s}^{\dagger},\hat{a}_{\bm{k}^{\prime},s^{\prime}}^{\dagger}]=0$.
The inverse transformation is given as

$\hat{a}_{j,s}=\frac{1}{L}\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{r}_{j}}\hat{a}_{\bm{k},s}.$
(23)

We define the time-dependent state vector:

$\ket{\psi(t)}:=\sum_{\bm{n}}P(\bm{n},t)\ket{\bm{n}}.$ (24)

Then, the master equation [Eq. (18)] is equivalent to the following imaginary-
time Schrödinger equation [27]:

$\frac{d}{dt}\ket{\psi(t)}=-\hat{H}\ket{\psi(t)},$ (25)

where the matrix element of the pseudo-Hamiltonian $\hat{H}$ is defined by

${}_{L}\\!\braket{\bm{n}}{\hat{H}}{\bm{n}^{\prime}}:=-W(\bm{n},\bm{n}^{\prime}).$
(26)

Thus, the steady-state vector,

$\ket{\psi_{\mathrm{st}}}:=\sum_{\bm{n}}P_{\mathrm{st}}(\bm{n})\ket{\bm{n}},$
(27)

is the eigenvector of $\hat{H}$ corresponding to the eigenvalue of zero. Since
$P_{\mathrm{st}}(\bm{n})={}_{L}\\!\braket{\bm{n}}{\psi_{\mathrm{st}}}$, the
calculation of the steady-state distribution $P_{\mathrm{st}}(\bm{n})$ is
equivalent to the eigenvalue problem of $\hat{H}$ (i.e., calculation of
$\ket{\psi_{\mathrm{st}}}$).

The steady-state average of a configuration-dependent physical quantity
$A(\bm{n})$ can be expressed [27] as

$\braket{A(\bm{n})}=\sum_{\bm{n}}A(\bm{n})P_{\mathrm{st}}(\bm{n})=\braket{\mathcal{P}}{A({\hat{\bm{n}}})}{\psi_{\mathrm{st}}},$
(28)

where $\hat{\bm{n}}:=\\{\hat{n}_{i,s}\\}$ and

$\bra{\mathcal{P}}:=\bra{0}e^{\sum_{i,s}\hat{a}_{i,s}}.$ (29)

Note that $\sum_{\bm{n}}P_{\mathrm{st}}(\bm{n})=1$ is equivalent to
$\braket{\mathcal{P}}{\psi_{\mathrm{st}}}=1$, as obtained from Eq. (28) by
setting $A(\bm{n})=1$.

Using this formulation, we can write the density correlation function [Eq.
(14)] as

$\displaystyle C_{d}(\bm{r}_{j})$
$\displaystyle=\frac{1}{L^{2}}\sum_{i}\braket{n_{i}n_{i+j}}-{\rho_{0}}^{2}$
$\displaystyle=\frac{1}{L^{2}}\sum_{i}\braket{\mathcal{P}}{\hat{n}_{i}\hat{n}_{i+j}}{\psi_{\mathrm{st}}}-{\rho_{0}}^{2}$
$\displaystyle=\frac{1}{L^{2}}\sum_{i}\sum_{s,s^{\prime}}\braket{\mathcal{P}}{\hat{a}_{i,s}\hat{a}_{i+j,s^{\prime}}}{\psi_{\mathrm{st}}}+\rho_{0}\delta_{\bm{r}_{j},\bm{0}}-{\rho_{0}}^{2},$
(30)

where $\hat{n}_{i}:=\sum_{s}\hat{n}_{i,s}$, $\rho_{0}=N/L^{2}$, and we used
$\bra{\mathcal{P}}\hat{a}_{i,s}^{\dagger}=\bra{\mathcal{P}}$. Similarly, the
polarity correlation function [Eq. (15)] is obtained as

$\displaystyle C_{p}(\bm{r}_{j})$
$\displaystyle=\frac{1}{L^{2}}\sum_{i}\braket{m_{i}m_{i+j}}$
$\displaystyle=\frac{1}{L^{2}}\sum_{i}\braket{\mathcal{P}}{\hat{m}_{i}\hat{m}_{i+j}}{\psi_{\mathrm{st}}}$
$\displaystyle=\frac{1}{L^{2}}\sum_{i}\sum_{s,s^{\prime}}ss^{\prime}\braket{\mathcal{P}}{\hat{a}_{i,s}\hat{a}_{i+j,s^{\prime}}}{\psi_{\mathrm{st}}}+\rho_{0}\delta_{\bm{r}_{j},\bm{0}},$
(31)

where $\hat{m}_{i}:=\sum_{s}s\hat{n}_{i,s}$.

We also obtain the corresponding density and polarity structure factors as

$\displaystyle S_{d}(\bm{k})$
$\displaystyle=\sum_{j}e^{-i\bm{k}\cdot\bm{r}_{j}}C_{d}(\bm{r}_{j})$
$\displaystyle=\sum_{s,s^{\prime}}\braket{\mathcal{P}}{\hat{a}_{\bm{k},s}\hat{a}_{-{\bm{k},s^{\prime}}}}{\psi_{\mathrm{st}}}+\rho_{0}-{\rho_{0}}^{2}L^{2}\delta_{\bm{k},\bm{0}}$
(32)

and

$\displaystyle S_{p}(\bm{k})$
$\displaystyle=\sum_{j}e^{-i\bm{k}\cdot\bm{r}_{j}}C_{p}(\bm{r}_{j})$
$\displaystyle=\sum_{s,s^{\prime}}ss^{\prime}\braket{\mathcal{P}}{\hat{a}_{\bm{k},s}\hat{a}_{-{\bm{k},s^{\prime}}}}{\psi_{\mathrm{st}}}+\rho_{0},$
(33)

respectively.

###  IV.2 Pseudo-Hamiltonian

We divide $\hat{H}$ [Eq. (26)] into two parts:

$\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{int}},$ (34)

where $\hat{H}_{0}$ is the unperturbed (i.e., non-interacting) part, and
$\hat{H}_{\mathrm{int}}$ is the perturbation (i.e., interaction) part.

$\hat{H}_{0}$ is obtained as

$\displaystyle\hat{H}_{0}$
$\displaystyle=-D_{0}\sum_{\braket{i,j},s}\hat{a}_{j,s}^{\dagger}\hat{a}_{i,s}-\varepsilon
D_{0}\sum_{i,s}s(\hat{a}_{i+\hat{x},s}^{\dagger}\hat{a}_{i,s}-\hat{a}_{i-\hat{x},s}^{\dagger}\hat{a}_{i,s})$
$\displaystyle\ \ \ \
-\gamma_{0}\sum_{i,s}\hat{a}_{i,-s}^{\dagger}\hat{a}_{i,s}+(4D_{0}+\gamma_{0})N,$
(35)

where $\sum_{\braket{i,j}}$ denotes the summation over all pairs of
neighboring sites $(i,j)$, with the pairs $(i,j)$ and $(j,i)$ being treated as
distinct, and $\hat{x}$ is the unit translation along the $x$ axis. The first
and third terms are Hermitian and represent the symmetric hopping and spin
flipping, respectively; the second term is anti-Hermitian and represents self-
propulsion; the last term is the diagonal part of $\hat{H}_{0}$, which is a
constant due to the conservation of the total particle number $N$.

$\hat{H}_{\mathrm{int}}$ is obtained as

$\displaystyle\hat{H}_{\mathrm{int}}$
$\displaystyle=D_{0}\sum_{\braket{i,j},s}(\hat{a}_{j,s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})\tanh\left[\frac{U(\hat{n}_{j}-\hat{n}_{i}+1)}{2}\right]$
$\displaystyle\ \ \ \ +\varepsilon
D_{0}\sum_{i,s}s(\hat{a}_{i+\hat{x},s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})\tanh\left[\frac{U(\hat{n}_{i+\hat{x}}-\hat{n}_{i}+1)}{2}\right]$
$\displaystyle\ \ \ \ -\varepsilon
D_{0}\sum_{i,s}s(\hat{a}_{i-\hat{x},s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})\tanh\left[\frac{U(\hat{n}_{i-\hat{x}}-\hat{n}_{i}+1)}{2}\right]$
$\displaystyle\ \ \ \
+\gamma_{0}\sum_{i,s}(\hat{a}_{i,-s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})\tanh\left[\frac{J(s\hat{m}_{i}-1)}{2}\right].$
(36)

Here, the first three lines and the last line represent the effects of
repulsion and alignment [see Eqs. (11)-(13)], respectively. Note that
$\hat{H}_{\mathrm{int}}$ is generically non-Hermitian.

###  IV.3 Non-Hermitian perturbation theory

Considering small $U$ and $J$, we expand $\hat{H}_{\mathrm{int}}$ [Eq. (36)]
up to the first order in terms of $U$ and $J$:

$\displaystyle\hat{H}_{\mathrm{int}}\simeq\hat{H}_{\mathrm{int}}^{\prime}$
$\displaystyle:=\frac{UD_{0}}{2}\sum_{\braket{i,j},s}(\hat{a}_{j,s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})(\hat{n}_{j}-\hat{n}_{i}+1)$
$\displaystyle\ \ \ \ +\frac{\varepsilon
UD_{0}}{2}\sum_{i,s}s(\hat{a}_{i+\hat{x},s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})(\hat{n}_{i+\hat{x}}-\hat{n}_{i}+1)$
$\displaystyle\ \ \ \ -\frac{\varepsilon
UD_{0}}{2}\sum_{i,s}s(\hat{a}_{i-\hat{x},s}^{\dagger}\hat{a}_{i,s}-\hat{n}_{i,s})(\hat{n}_{i-\hat{x}}-\hat{n}_{i}+1)$
$\displaystyle\ \ \ \
+\frac{J\gamma_{0}}{2}\sum_{i,s}(a_{i,-s}^{\dagger}a_{i,s}-\hat{n}_{i,s})(s\hat{m}_{i}-1).$
(37)

Then, we expand the steady-state vector $\ket{\psi_{\mathrm{st}}}$ up to the
first order in terms of $\hat{H}_{\mathrm{int}}^{\prime}$:

$\ket{\psi_{\mathrm{st}}}\simeq\ket{\psi_{\mathrm{st}}^{(0)}}+\ket{\psi_{\mathrm{st}}^{(1)}},$
(38)

where $\ket{\psi_{\mathrm{st}}^{(0)}}$ is the steady-state vector of
$\hat{H}_{0}$, and $\ket{\psi_{\mathrm{st}}^{(1)}}$ is the first-order
correction. Applying perturbation theory for non-Hermitian systems [45], we
obtain

$\ket{\psi_{\mathrm{st}}^{(1)}}=-\sum_{n\neq
0}\frac{{}_{L}\\!\braket{\psi^{(0)}(n)}{\hat{H}_{\mathrm{int}}^{\prime}}{\psi_{\mathrm{st}}^{(0)}}}{E^{(0)}(n)}\ket{\psi^{(0)}(n)}.$
(39)

Here, $\ket{\psi^{(0)}(n)}$, ${}_{L}\\!\bra{\psi^{(0)}(n)}$, and $E^{(0)}(n)$
are the eigenvector, left eigenvector, and eigenvalue of $\hat{H}_{0}$
specified by $n$, respectively, and $n=0$ corresponds to the steady state
[e.g., $\ket{\psi^{(0)}(0)}=\ket{\psi_{\mathrm{st}}^{(0)}}$]. We assume
biorthogonality:
${}_{L}\\!\braket{\psi^{(0)}(n)}{\psi^{(0)}(m)}=\delta_{n,m}$.

To calculate $\ket{\psi_{\mathrm{st}}^{(1)}}$ using Eq. (39), we need to
obtain $\ket{\psi_{\mathrm{st}}^{\mathrm{(0)}}}$ and
$\\{\ket{\psi^{(0)}(n)},{}_{L}\\!\bra{\psi^{(0)}(n)},E^{(0)}(n)\\}$ for all
$n$ by solving the eigenvalue problem for $\hat{H}_{0}$. In the following, we
employ linear transformations of the annihilation and creation operators that
preserve the boson commutation relation.

From Eqs. (23) and (35), we can express $\hat{H}_{0}$ as

$\hat{H}_{0}=\sum_{\bm{k}}\begin{pmatrix}\hat{a}_{\bm{k},+}^{\dagger}&\hat{a}_{\bm{k},-}^{\dagger}\end{pmatrix}h_{\bm{k}}\begin{pmatrix}\hat{a}_{\bm{k},+}\\\
\hat{a}_{\bm{k},-}\end{pmatrix},$ (40)

where

$h_{\bm{k}}=\begin{pmatrix}2D_{0}(2-\cos k_{x}-\cos k_{y}+i\varepsilon\sin
k_{x})+\gamma_{0}&-\gamma_{0}\\\ -\gamma_{0}&2D_{0}(2-\cos k_{x}-\cos
k_{y}-i\varepsilon\sin k_{x})+\gamma_{0}\end{pmatrix}.$ (41)

Solving the eigenvalue problem of $2\times 2$ matrix $h_{\bm{k}}$, we obtain
the eigenvector $\bm{u}_{\bm{k},\alpha}$, left eigenvector
$\bm{v}_{\bm{k},\alpha}^{T}$, and eigenvalue $\epsilon_{\bm{k},\alpha}$
($\alpha\in\\{0,1\\}$,
$h_{\bm{k}}\bm{u}_{\bm{k},\alpha}=\epsilon_{\bm{k},\alpha}\bm{u}_{\bm{k},\alpha}$,
and
$\bm{v}_{\bm{k},\alpha}^{T}h_{\bm{k}}=\bm{v}_{\bm{k},\alpha}^{T}\epsilon_{\bm{k},\alpha}$):

$\displaystyle\bm{u}_{\bm{k},0}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
\sqrt{1-{f_{\bm{k}}}^{2}}+if_{\bm{k}}\end{pmatrix}$ (42)
$\displaystyle\bm{u}_{\bm{k},1}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
-\sqrt{1-{f_{\bm{k}}}^{2}}+if_{\bm{k}}\end{pmatrix},$ (43)
$\displaystyle\bm{v}_{\bm{k},0}$
$\displaystyle=\frac{1}{\sqrt{2}}\frac{1}{\sqrt{1-{f_{\bm{k}}}^{2}}}\begin{pmatrix}\sqrt{1-{f_{\bm{k}}}^{2}}-if_{\bm{k}}\\\
1\end{pmatrix}$ (44) $\displaystyle\bm{v}_{\bm{k},1}$
$\displaystyle=\frac{1}{\sqrt{2}}\frac{1}{\sqrt{1-{f_{\bm{k}}}^{2}}}\begin{pmatrix}\sqrt{1-{f_{\bm{k}}}^{2}}+if_{\bm{k}}\\\
-1\end{pmatrix},$ (45)

and

$\displaystyle\epsilon_{\bm{k},0}$ $\displaystyle=2D_{0}(2-\cos k_{x}-\cos
k_{y})+\gamma_{0}(1-\sqrt{1-{f_{\bm{k}}}^{2}})$ (46)
$\displaystyle\epsilon_{\bm{k},1}$ $\displaystyle=2D_{0}(2-\cos k_{x}-\cos
k_{y})+\gamma_{0}(1+\sqrt{1-{f_{\bm{k}}}^{2}}),$ (47)

where $f_{\bm{k}}:=(2\varepsilon D_{0}/\gamma_{0})\sin k_{x}$ and
$\sqrt{1-{f_{\bm{k}}}^{2}}=i\sqrt{{f_{\bm{k}}}^{2}-1}$ for $|f_{\bm{k}}|>1$.
Note that $\bm{u}_{\bm{k},\alpha}$ and $\bm{v}_{\bm{k},\alpha}^{T}$ satisfy
biorthogonality:
$\bm{v}_{\bm{k},\alpha}^{T}\bm{u}_{\bm{k}^{\prime},\alpha^{\prime}}=\delta_{\bm{k},\bm{k}^{\prime}}\delta_{\alpha,\alpha^{\prime}}$.

We define a linear transformation of $\hat{a}_{\bm{k},s}$ and
$\hat{a}_{\bm{k},s}^{\dagger}$ as

$\displaystyle\hat{b}_{\bm{k},\alpha}$
$\displaystyle:=\bm{v}_{\bm{k},\alpha}^{T}\begin{pmatrix}\hat{a}_{\bm{k},+}\\\
\hat{a}_{\bm{k},-}\end{pmatrix}$ (48)
$\displaystyle\hat{\bar{b}}_{\bm{k},\alpha}$
$\displaystyle:=\begin{pmatrix}\hat{a}_{\bm{k},+}^{\dagger}&\hat{a}_{\bm{k},-}^{\dagger}\end{pmatrix}\bm{u}_{\bm{k},\alpha}.$
(49)

From the biorthogonality of $\bm{u}_{\bm{k},\alpha}$ and
$\bm{v}_{\bm{k},\alpha}^{T}$, the inverse transformation is given as

$\displaystyle\begin{pmatrix}\hat{a}_{\bm{k},+}\\\
\hat{a}_{\bm{k},-}\end{pmatrix}=\sum_{\alpha}\bm{u}_{\bm{k},\alpha}\hat{b}_{\bm{k},\alpha}$
(50)
$\displaystyle\begin{pmatrix}\hat{a}_{\bm{k},+}^{\dagger}&\hat{a}_{\bm{k},-}^{\dagger}\end{pmatrix}=\sum_{\alpha}\hat{\bar{b}}_{\bm{k},\alpha}\bm{v}_{\bm{k},\alpha}^{T}.$
(51)

Note that $\hat{\bar{b}}_{\bm{k},\alpha}\neq\hat{b}_{\bm{k},\alpha}^{\dagger}$
in general, except the case with no self-propulsion (i.e., $\varepsilon=0$),
where $f_{\bm{k}}=0$ and $\bm{u}_{\bm{k},\alpha}=\bm{v}_{\bm{k},\alpha}$.
However, since the commutation relation holds (i.e.,
$[\hat{b}_{\bm{k},\alpha},\hat{\bar{b}}_{\bm{k}^{\prime},\alpha^{\prime}}]=\delta_{\bm{k},\bm{k}^{\prime}}\delta_{\alpha,\alpha^{\prime}}$
and
$[\hat{b}_{\bm{k},\alpha},\hat{b}_{\bm{k}^{\prime},\alpha^{\prime}}]=[\hat{\bar{b}}_{\bm{k},\alpha},\hat{\bar{b}}_{\bm{k}^{\prime},\alpha^{\prime}}]=0$),
$\hat{b}_{\bm{k},\alpha}$, $\hat{\bar{b}}_{\bm{k},\alpha}$, and
$\hat{\bar{b}}_{\bm{k},\alpha}\hat{b}_{\bm{k},\alpha}$ are the annihilation,
creation, and number operators for a quasiparticle specified by
$(\bm{k},\alpha)$.

Using the spectral decomposition
$h_{\bm{k}}=\sum_{\alpha}\epsilon_{\bm{k},\alpha}\bm{u}_{\bm{k},\alpha}\bm{v}_{\bm{k},\alpha}^{T}$
in Eq. (40), we obtain

$\hat{H}_{0}=\sum_{\bm{k},\alpha}\epsilon_{\bm{k},\alpha}\hat{\bar{b}}_{\bm{k},\alpha}\hat{b}_{\bm{k},\alpha}.$
(52)

Thus, any eigenvector of $\hat{H}_{0}$ is specified by the set of occupancies
for all quasiparticle states, which is denoted by $\\{n_{\bm{k},\alpha}\\}$.
The eigenvector $\ket{\psi^{(0)}(\\{n_{\bm{k},\alpha}\\})}$ is given as

$\ket{\psi^{(0)}(\\{n_{\bm{k},\alpha}\\})}=\prod_{\bm{k},\alpha}\left(\hat{\bar{b}}_{\bm{k},\alpha}\right)^{n_{\bm{k},\alpha}}\ket{0}.$
(53)

The corresponding left eigenvector
${}_{L}\\!\bra{\psi^{(0)}(\\{n_{\bm{k},\alpha}\\})}$ and eigenvalue
$E^{(0)}(\\{n_{\bm{k},\alpha}\\})$ are

${}_{L}\\!\bra{\psi^{(0)}(\\{n_{\bm{k},\alpha}\\})}=\bra{0}\prod_{\bm{k},\alpha}\frac{1}{n_{\bm{k},\alpha}!}\left(\hat{b}_{\bm{k},\alpha}\right)^{n_{\bm{k},\alpha}}$
(54)

and

$E^{(0)}(\\{n_{\bm{k},\alpha}\\})=\sum_{\bm{k},\alpha}\epsilon_{\bm{k},\alpha}n_{\bm{k},\alpha},$
(55)

respectively. Note that $\sum_{\bm{k},\alpha}n_{\bm{k},\alpha}=N$ since
$\sum_{\bm{k},\alpha}\hat{\bar{b}}_{\bm{k},\alpha}\hat{b}_{\bm{k},\alpha}=\sum_{i,s}\hat{a}_{i,s}^{\dagger}\hat{a}_{i,s}=N$
is the fixed total particle number.

From Eqs. (53) and (55), we find that the steady-state vector of
$\hat{H}_{0}$, which corresponds to the eigenvalue of zero [i.e.,
$E^{(0)}(\\{n_{\bm{k},\alpha}\\})=0$], is obtained by taking $n_{\bm{0},0}=N$
and $n_{\bm{k},\alpha}=0$ for all $(\bm{k},\alpha)\neq(\bm{0},0)$:

$\ket{\psi_{\mathrm{st}}^{(0)}}=\frac{1}{(\sqrt{2}L)^{N}}\left(\hat{\bar{b}}_{\bm{0},0}\right)^{N}\ket{0},$
(56)

where the normalization factor $(\sqrt{2}L)^{-N}$ ensures that the total
probability for the configuration is $1$ (i.e.,
$\braket{\mathcal{P}}{\psi_{\mathrm{st}}^{(0)}}=\sum_{\bm{n}}{}_{L}\\!\braket{\bm{n}}{\psi_{\mathrm{st}}^{(0)}}=1$)
[see the notes after Eq. (29)]. This steady-state vector represents the random
distribution of $N$ particles with no spatial or spin correlations since
$\hat{\bar{b}}_{\bm{0},0}\propto\sum_{i,s}\hat{a}_{i,s}^{\dagger}$, which
follows from Eqs. (22), (42), and (49). All the other eigenvectors are
specified by the set of occupancies $\\{n_{\bm{k},\alpha}\\}$ with
$n_{\bm{0},0}<N$. Thus, $\ket{\psi_{\mathrm{st}}^{(1)}}$ in Eq. (39) is
expressed as

$\ket{\psi_{\mathrm{st}}^{(1)}}=-\sideset{}{{}^{\prime}}{\sum}_{\\{n_{\bm{k},\alpha}\\}}\frac{{}_{L}\\!\braket{\psi^{(0)}(\\{n_{\bm{k},\alpha}\\})}{\hat{H}_{\mathrm{int}}^{\prime}}{\psi_{\mathrm{st}}^{(0)}}}{E^{(0)}(\\{n_{\bm{k},\alpha}\\})}\ket{\psi^{(0)}(\\{n_{\bm{k},\alpha}\\})},$
(57)

where $\sum_{\\{n_{\bm{k},\alpha}\\}}^{\prime}$ denotes the summation over all
sets of occupancies $\\{n_{\bm{k},\alpha}\\}$ satisfying
$\sum_{\bm{k},\alpha}n_{\bm{k},\alpha}=N$ and $n_{\bm{0},0}<N$.

To reduce Eq. (57) to a more explicit formula, we rewrite
$\hat{H}_{\mathrm{int}}^{\prime}$ [Eq. (37)] using $\hat{b}_{\bm{k},\alpha}$
and $\hat{\bar{b}}_{\bm{k},\alpha}$. From Eqs. (23), (37), (50), and (51), we
obtain

$\displaystyle\hat{H}_{\mathrm{int}}^{\prime}$
$\displaystyle=\frac{1}{L^{2}}\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\bm{q}\\\
s_{1},s_{2},s_{3},s_{4}\end{subarray}}\\!\\!\\!\\!\\!V_{s_{1},s_{2},s_{3},s_{4}}(\bm{k},\bm{k}^{\prime},\bm{q})\hat{a}_{\bm{q}+\bm{k},s_{1}}^{\dagger}\hat{a}_{\bm{q}-\bm{k},s_{2}}^{\dagger}\hat{a}_{\bm{q}-\bm{k}^{\prime},s_{4}}\hat{a}_{\bm{q}+\bm{k}^{\prime},s_{3}}$
$\displaystyle=\frac{1}{L^{2}}\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\bm{q}\\\
\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\end{subarray}}\\!\\!\\!\\!\\!\\!\tilde{V}_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}}(\bm{k},\bm{k}^{\prime},\bm{q})\hat{\bar{b}}_{\bm{q}+\bm{k},\alpha_{1}}\hat{\bar{b}}_{\bm{q}-\bm{k},\alpha_{2}}\hat{b}_{\bm{q}-\bm{k}^{\prime},\alpha_{4}}\hat{b}_{\bm{q}+\bm{k}^{\prime},\alpha_{3}},$
(58)

where

$\displaystyle V_{s_{1},s_{2},s_{3},s_{4}}(\bm{k},\bm{k}^{\prime},\bm{q})$
$\displaystyle=UD_{0}\delta_{s_{1},s_{3}}\delta_{s_{2},s_{4}}[2-\cos(q_{x}+k_{x})+\cos(q_{x}+k_{x}^{\prime})-\cos(k_{x}-k_{x}^{\prime})-\cos(q_{y}+k_{y})+\cos(q_{y}+k_{y}^{\prime})-\cos(k_{y}-k_{y}^{\prime})]$
$\displaystyle\ \ \ \ +is_{3}\varepsilon
UD_{0}\delta_{s_{1},s_{3}}\delta_{s_{2},s_{4}}[\sin(q_{x}+k_{x})-\sin(q_{x}+k_{x}^{\prime})-\sin(k_{x}-k_{x}^{\prime})]+\frac{J\gamma_{0}}{2}s_{3}s_{4}(\delta_{s_{1},-s_{3}}-\delta_{s_{1},s_{3}})\delta_{s_{2},s_{4}}$
(59)
$\displaystyle\tilde{V}_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}}(\bm{k},\bm{k}^{\prime},\bm{q})$
$\displaystyle=\sum_{s_{1},s_{2},s_{3},s_{4}}V_{s_{1},s_{2},s_{3},s_{4}}(\bm{k},\bm{k}^{\prime},\bm{q})v_{\bm{q}+\bm{k},\alpha_{1},s_{1}}v_{\bm{q}-\bm{k},\alpha_{2},s_{2}}u_{\bm{q}-\bm{k}^{\prime},\alpha_{4},s_{4}}u_{\bm{q}+\bm{k}^{\prime},\alpha_{3},s_{3}}.$
(60)

Note that $u_{\bm{k},\alpha,s}$ and $v_{\bm{k},\alpha,s}$ denote the
components of $\bm{u}_{\bm{k},\alpha}$ and $\bm{v}_{\bm{k},\alpha}$,
respectively.

According to the forms of $\ket{\psi_{\mathrm{st}}^{(0)}}$ [Eq. (56)] and
$\hat{H}_{\mathrm{int}}^{\prime}$ [Eq. (58)], the first-order perturbation
[Eq. (57)] causes two-quasiparticle excitations in the language of quantum
theory. From Eqs. (53)-(58), we obtain

$\ket{\psi_{\mathrm{st}}^{(1)}}=-\frac{N(N-1)}{(\sqrt{2}L)^{N}L^{2}}\sideset{}{{}^{\prime}}{\sum}_{\bm{k},\alpha_{1},\alpha_{2}}\frac{\tilde{V}_{\alpha_{1},\alpha_{2},0,0}(\bm{k},\bm{0},\bm{0})}{\epsilon_{\bm{k},\alpha_{1}}+\epsilon_{-\bm{k},\alpha_{2}}}\ket{\bm{k},\alpha_{1};-\bm{k},\alpha_{2}}.$
(61)

Here, $\sum_{\bm{k},\alpha_{1},\alpha_{2}}^{\prime}$ denotes the summation
over all sets of
$(\bm{k},\alpha_{1},\alpha_{2})\notin\\{(\bm{0},0,0),(\bm{0},0,1),(\bm{0},1,0)\\}$,
and

$\ket{\bm{k},\alpha_{1};-\bm{k},\alpha_{2}}:=\hat{\bar{b}}_{\bm{k},\alpha_{1}}\hat{\bar{b}}_{-\bm{k},\alpha_{2}}\left(\hat{\bar{b}}_{\bm{0},0}\right)^{N-2}\ket{0},$
(62)

which corresponds to the excitation of two quasiparticles specified by
$(\bm{k},\alpha_{1})$ and $(-\bm{k},\alpha_{2})$ from the unperturbed steady
state $\ket{\psi_{\mathrm{st}}^{(0)}}$. For completeness, we give the
expression of $\tilde{V}_{\alpha_{1},\alpha_{2},0,0}(\bm{k},\bm{0},\bm{0})$:

$\displaystyle\tilde{V}_{\alpha_{1},\alpha_{2},0,0}(\bm{k},\bm{0},\bm{0})$
$\displaystyle=\sum_{s_{1},s_{2}}\left[UD_{0}(2-\cos k_{x}-\cos
k_{y})-\frac{J\gamma_{0}}{2}s_{1}s_{2}\right]$ $\displaystyle\ \ \ \ \ \ \ \ \
\ \ \times v_{\bm{k},\alpha_{1},s_{1}}v_{-\bm{k},\alpha_{2},s_{2}}.$ (63)

We now have the analytical expression of the steady-state vector
$\ket{\psi_{\mathrm{st}}}$ up to the first order of the interaction parameters
$U$ and $J$ [Eqs. (56) and (61)]. Since $\ket{\psi_{\mathrm{st}}}$ is
equivalent to the steady-state distribution through
$P_{\mathrm{st}}(\bm{n})={}_{L}\\!\braket{\bm{n}}{\psi_{\mathrm{st}}}$, we can
calculate physical quantities using the obtained formula, as we explain in the
next section.

###  IV.4 Correlation functions and structure factors

Using the expression of
$\ket{\psi_{\mathrm{st}}}\simeq\ket{\psi_{\mathrm{st}}^{(0)}}+\ket{\psi_{\mathrm{st}}^{(1)}}$
obtained by perturbation theory [Eqs. (56) and (61)], we expand the density
structure factor $S_{d}(\bm{k})$ [Eq. (32)] up to the first order in $U$ and
$J$:

$S_{d}(\bm{k})\simeq S_{d}^{(0)}(\bm{k})+S_{d}^{(1)}(\bm{k}).$ (64)

For $\bm{k}\neq\bm{0}$, the unperturbed structure factor $S_{d}^{(0)}(\bm{k})$
is given as

$S_{d}^{(0)}(\bm{k})=\rho_{0},$ (65)

and the first-order correction $S_{d}^{(1)}(\bm{k})$ is obtained as

$S_{d}^{(1)}(\bm{k})=-\frac{N(N-1)}{L^{4}}\sum_{s,s^{\prime}}\sum_{\alpha,\alpha^{\prime}}\frac{u_{\bm{k},\alpha,s}u_{-\bm{k},\alpha^{\prime},s^{\prime}}\tilde{V}_{\alpha,\alpha^{\prime},0,0}(\bm{k},\bm{0},\bm{0})}{\epsilon_{\bm{k},\alpha}+\epsilon_{-\bm{k},\alpha^{\prime}}}.$
(66)

See Eqs. (42)-(47) and (63) for the expressions of $\bm{u}_{\bm{k},\alpha}$,
$\tilde{V}_{\alpha,\alpha^{\prime},0,0}(\bm{k},\bm{0},\bm{0})$, and
$\epsilon_{\bm{k},\alpha}$.

For $\bm{k}=\bm{0}$, we can check $S_{d}^{(0)}(\bm{0})=S_{d}^{(1)}(\bm{0})=0$,
consistent with $S_{d}(\bm{0})=0$, which generally holds from the particle
number conservation (i.e., $\sum_{i}n_{i}=N$). Note that $S_{d}^{(0)}(\bm{k})$
[$=\rho_{0}(1-\delta_{\bm{k},\bm{0}})$] leads to the unperturbed density
correlation function:

$C_{d}^{(0)}(\bm{r}):=\frac{1}{L^{2}}\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}S_{d}(\bm{k})=\rho_{0}\delta_{\bm{r},\bm{0}}-\frac{\rho_{0}}{L^{2}},$
(67)

which includes a term $-\rho_{0}/L^{2}$ due to the finite-size effect. In
Figs. 3(e) and 4(e), the two lines with light color represent the observed
density correlation function after subtracting this term (i.e.,
$C_{d}(\bm{r})+\rho_{0}/L^{2}$).

Similarly, we expand the polarity structure factor as

$S_{p}(\bm{k})\simeq S_{p}^{(0)}(\bm{k})+S_{p}^{(1)}(\bm{k}),$ (68)

where

$S_{p}^{(0)}(\bm{k})=\rho_{0},$ (69)

and

$S_{p}^{(1)}(\bm{k})=-\frac{N(N-1)}{L^{4}}\sum_{s,s^{\prime}}\sum_{\alpha,\alpha^{\prime}}ss^{\prime}\frac{u_{\bm{k},\alpha,s}u_{-\bm{k},\alpha^{\prime},s^{\prime}}\tilde{V}_{\alpha,\alpha^{\prime},0,0}(\bm{k},\bm{0},\bm{0})}{\epsilon_{\bm{k},\alpha}+\epsilon_{-\bm{k},\alpha^{\prime}}}.$
(70)

We examine the long-wavelength behaviors of $S_{d}(\bm{k})$ and
$S_{p}(\bm{k})$, which determine the long-distance behaviors of the
correlation functions $C_{d}(\bm{r})$ and $C_{p}(\bm{r})$, respectively. The
terms with $(\alpha,\alpha^{\prime})=(0,0)$ in $S_{d}^{(1)}(\bm{k})$ [Eq.
(66)] and $S_{p}^{(1)}(\bm{k})$ [Eq. (70)] can be non-analytic at
$\bm{k}\to\bm{0}$, while the other terms with
$(\alpha,\alpha^{\prime})\neq(0,0)$ are analytic. Noticing that
$\tilde{V}_{0,0,0,0}(\bm{k},\bm{0},\bm{0})\simeq
UD_{0}\bm{k}^{2}-(J\varepsilon^{2}{D_{0}}^{2}/\gamma_{0}){k_{x}}^{2}$ for
$\bm{k}\simeq\bm{0}$, we can expand $S_{d}^{(1)}(\bm{k})$ and
$S_{p}^{(1)}(\bm{k})$ for $\bm{k}\simeq\bm{0}$:

$\displaystyle S_{d}^{(1)}(\bm{k})$
$\displaystyle\simeq-\frac{N(N-1)}{L^{4}}\frac{U\bm{k}^{2}-J(\varepsilon^{2}D_{0}/\gamma_{0}){k_{x}}^{2}}{\bm{k}^{2}+(2\varepsilon^{2}D_{0}/\gamma_{0}){k_{x}}^{2}}+(\mathrm{a.t.})$
(71) $\displaystyle S_{p}^{(1)}(\bm{k})$
$\displaystyle\simeq\left(\frac{\varepsilon
D_{0}}{\gamma_{0}}k_{x}\right)^{2}S_{d}^{(1)}(\bm{k})+(\mathrm{a.t.})$
$\displaystyle=-\frac{N(N-1)}{L^{4}}\left(\frac{\varepsilon
D_{0}}{\gamma_{0}}k_{x}\right)^{2}\frac{U\bm{k}^{2}-J(\varepsilon^{2}D_{0}/\gamma_{0}){k_{x}}^{2}}{\bm{k}^{2}+(2\varepsilon^{2}D_{0}/\gamma_{0}){k_{x}}^{2}}+(\mathrm{a.t.}),$
(72)

where $(\mathrm{a.t.})$ denotes the analytic terms derived from
$(\alpha,\alpha^{\prime})\neq(0,0)$. The first line in Eq. (72) indicates the
coupling of density and polarity at the level of first-order perturbation (see
Sec. II for the hydrodynamic counterpart).

Figure 5: Singularity in structure factors induced by perturbative repulsion.
(a, b) Density structure factor $S_{d}(\bm{k})$ as a function of (a) $k_{x}$
and (b) $k_{y}$, for several fixed values of $k_{y}$ and $k_{x}$,
respectively. (c, d) Polarity structure factor $S_{p}(\bm{k})$, plotted in the
same way as done in (a) and (b). The solid and dashed lines in each panel
represent the values obtained from simulations and those predicted from
perturbation theory. The predicted $S_{d}(\bm{k})$ and $S_{p}(\bm{k})$ are
discontinuous and non-smooth at $\bm{k}\to\bm{0}$, respectively. The used
parameters are $D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=0.5$, $U=0.005$,
$L=100$, and $N=4\times 10^{4}$. Figure 6: Singularity in structure factors
induced by perturbative alignment. In the same way as shown in Fig. 5, we plot
(a, b) the density structure factor $S_{d}(\bm{k})$ and (c, d) the polarity
structure factor $S_{p}(\bm{k})$. The values obtained from simulations (solid
lines) and those predicted from perturbation theory (dashed lines) show good
agreement, suggesting singularities in $S_{d}(\bm{k})$ and $S_{p}(\bm{k})$
similarly to the case of repulsion (Fig. 5). The used parameters are
$D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=2.0$, $J=0.005$, $L=100$, and
$N=4\times 10^{4}$.

From Eq. (71), we find that $S_{d}(\bm{k})$ is generically discontinuous at
$\bm{k}\to\bm{0}$, which is characterized by the difference in the two
limiting values:

$\displaystyle\lim_{k_{x}\to 0}S_{d}(k_{x},0)-\lim_{k_{y}\to 0}S_{d}(0,k_{y})$
$\displaystyle=\frac{N(N-1)}{L^{4}}\frac{(2U+J)\varepsilon^{2}D_{0}/\gamma_{0}}{1+2\varepsilon^{2}D_{0}/\gamma_{0}},$
(73)

up to the first order in the repulsion $U$ ($\geq 0$) and alignment $J$ ($\geq
0$). This formula suggests that the nonzero anisotropic self-propulsion
$\varepsilon$ and interaction $U$ or $J$ are essential, but the interaction
type (i.e., repulsion or alignment) is not essential, to the discontinuity of
$S_{d}(\bm{k})$. In addition, the prefactor $N(N-1)$ suggests that the
discontinuity appears even at the two-particle level (i.e., $N=2$). As in the
case of the hydrodynamic model (see Sec. II), from
$C_{d}(\bm{r})=L^{-2}\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}S_{d}(\bm{k})\simeq\int
d^{2}\bm{k}e^{i\bm{k}\cdot\bm{r}}S_{d}(\bm{k})$, we find that the density
correlation function follows a power law for $|\bm{r}|\to\infty$ reflecting
the discontinuity of $S_{d}(\bm{k})$:

$\begin{split}C_{d}(x,0)&\sim x^{-2}\\\ C_{d}(0,y)&\sim y^{-2}.\end{split}$
(74)

Note that the power-law density correlation in two-particle systems has also
been found in externally driven systems [46, 47].

Focusing on Eq. (72), we see that $S_{p}(\bm{k})$ is continuous but can be
non-smooth at $\bm{k}\to\bm{0}$, which is characterized by

$\displaystyle\lim_{k_{x}\to 0}\frac{\partial^{2}S_{p}(\bm{k})}{{\partial
k_{y}}^{2}}\Bigg{|}_{(k_{x},0)}-\lim_{k_{y}\to
0}\frac{\partial^{2}S_{p}(\bm{k})}{{\partial k_{y}}^{2}}\Bigg{|}_{(0,k_{y})}$
$\displaystyle=-\frac{N(N-1)}{L^{4}}\frac{2(2U+J)\varepsilon^{4}(D_{0}/\gamma_{0})^{3}}{(1+2\varepsilon^{2}D_{0}/\gamma_{0})^{2}},$
(75)

up to the first order in $U$ and $J$. The condition for the non-smoothness of
$S_{p}(\bm{k})$ is the same as that for the discontinuity of $S_{d}(\bm{k})$,
i.e., nonzero anisotropic self-propulsion and interaction (regardless of
repulsion or alignment). From
$C_{p}(\bm{r})=L^{-2}\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}S_{p}(\bm{k})\simeq\int
d^{2}\bm{k}e^{i\bm{k}\cdot\bm{r}}S_{p}(\bm{k})$, we find that the power-law
correlation of polarity appears for $|\bm{r}|\to\infty$ reflecting the non-
smoothness of $S_{p}(\bm{k})$:

$\begin{split}C_{p}(x,0)&\sim x^{-4}\\\ C_{p}(0,y)&\sim y^{-4}.\end{split}$
(76)

###  IV.5 Comparison with numerical results

In the following, we confirm that the obtained formulas [Eqs. (64)-(66) and
Eqs. (68)-(70)] quantitatively reproduce the numerical results when $U$ and
$J$ are small.

We first consider a purely repulsive interaction (i.e., $U>0$ and $J=0$). The
parameters are set as $D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=0.5$,
$U=0.005$, $L=100$, and $N=4\times 10^{4}$ (corresponding to $\rho_{0}=4$). In
Figs. 5(a) and (b), we compare the numerically obtained density structure
factor $S_{d}(\bm{k})$ (solid lines) and the analytical results based on Eqs.
(64)-(66) (dotted lines). Both the (a) $k_{x}$ dependence and (b) $k_{y}$
dependence of $S_{d}(\bm{k})$ are reproduced well by the analytical formulas,
especially when $S_{d}(\bm{k})$ is close to the unperturbed value
$S_{d}^{(0)}(\bm{k})=\rho_{0}=4$. In Figs. 5(c) and (d), we compare the
numerically obtained polarity structure factor $S_{p}(\bm{k})$ (solid lines)
and the analytical results based on Eqs. (68)-(70) (dotted lines), which also
show good agreement.

We next consider a purely aligning interaction (i.e., $U=0$ and $J>0$). The
parameters are $D_{0}=1$, $\varepsilon=0.8$, $\gamma_{0}=2$, $J=0.005$,
$L=100$, and $N=4\times 10^{4}$ (corresponding to $\rho_{0}=4$). In Fig. 6,
the figures are placed in the same order as Fig. 5. In all figures, the
numerical results agree with the analytical predictions with high accuracy.

##  V Discussion and summary

In this paper, we have studied the correlation properties in the homogeneous
disordered state of particle systems with uniaxial self-propulsion. Starting
with hydrodynamic arguments, we have proposed that such systems should
generically show power-law density correlation and polarity correlation with
exponent $-2$ and $-4$, respectively. Performing simulations of a lattice gas
model with uniaxial self-propulsion and repulsion or alignment between
particles, we have found that the predicted power-law correlations indeed
appear regardless of the interaction type. Further, using the Doi-Peliti
method and perturbation theory, we have mapped the model to a two-component
boson system and analyzed the effects of interaction as quasiparticle
excitations. We have analytically obtained the formulas for the density and
polarity structure factors, which have singularities that lead to the power-
law decay of correlation functions, even in the first order of the interaction
strength. Both the hydrodynamic and microscopic formulas [Eqs. (4) and (72)]
suggest that the coupling of density and polarity is essential to the power-
law correlation of polarity.

We have focused on the homogeneous disordered state, which generically appears
when the interaction strength is weak compared to the noise effect [Figs. 2(a)
and (c)]. When the interaction is strong, phase separation or long-range order
can appear through phase transitions such as MIPS and flocking [Figs. 2(b) and
(d)]. It is interesting to examine whether the power-law correlation for
density fluctuation or polarity fluctuation still appears in such
inhomogeneous or ordered states.

According to our results, in anisotropic biological systems, we may observe
the power-law correlations of density or polarity as a universal
nonequilibrium collective phenomenon, regardless of the detail of
interactions. For example, spatial anisotropy can be introduced as an
anisotropic interaction between a cell population and orientationally aligned
substrate [48, 49, 50]. We note that observation of the power law using a few
samples can be hard, given that the correlation functions decay relatively
quickly since the exponents are $-2$ and $-4$ for density and polarity,
respectively [i.e., $C_{d}(r)\sim r^{-2}$ and $C_{p}(r)\sim r^{-4}$].
Nevertheless, it will be interesting to probe nonequilibrium properties
inherent in biological systems by adding spatial anisotropy.

###### Acknowledgements.

We thank Ryusuke Hamazaki, Kyogo Kawaguchi, Hiroshi Noguchi, Yuta Sekino,
Kazuaki Takasan, Takeru Yokota for fruitful discussions. The computations in
this study were performed using the facilities of the Supercomputer Center at
the Institute for Solid State Physics, the University of Tokyo. This work was
supported by JSPS KAKENHI Grant Numbers JP20K14435 (to K.A.) and JP22K13978
(to H.N.).

Figure 7: Relaxation dynamics of structure factors for the model with
repulsion. We plot the time evolution of the long-wavelength components of the
structure factors $S_{d}(\bm{k})$ and $S_{p}(\bm{k})$ averaged over
independent samples for the parameter sets used in (a) Fig. 3 and (b) Fig. 5.
The gray area in each panel suggests the time points used for averaging.
Figure 8: Relaxation dynamics of structure factors for the model with
alignment. We plot the time evolution of the long-wavelength components of the
structure factors $S_{d}(\bm{k})$ and $S_{p}(\bm{k})$ averaged over
independent samples for the parameter sets used in (a) Fig. 4 and (b) Fig. 6.
The gray area in each panel suggests the time points used for averaging.

##  Appendix A Numerical simulation of lattice gas model

For the simulation of the lattice gas model explained in Sec. III.1, we
discretize the time with interval $\Delta t$. For each sample, we generate the
initial configuration of $N$ particles by placing each particle on a randomly
chosen site with a randomly chosen spin. In a single Monte Carlo (MC) step, we
update the configuration as follows.

  1. (1)

We randomly choose a particle.

  2. (2)

The chosen particle (with spin $s$ at site $i$) (i) hops by one site or (ii)
flips the spin with the following probabilities.

    1. (i)

The particle hops one site right with probability
$2(1+s\varepsilon)D_{0}\Delta t/[1+e^{U(\Delta n_{i\to j}+1)}]$, left with
probability $2(1-s\varepsilon)D_{0}\Delta t/[1+e^{U(\Delta n_{i\to j}+1)}]$,
or up or down with probability $2D_{0}\Delta t/[1+e^{U(\Delta n_{i\to j}+1)}]$
for each, where $\Delta n_{i\to j}:=n_{j}-n_{i}$ is the difference in the
local density between the target site $j$ and the departure site $i$.

    2. (ii)

The spin of the particle flips to $-s$ with probability $2\gamma_{0}\Delta
t/[1+e^{J(sm_{i}-1)}]$, where
$m_{i}=\sum_{s^{\prime}}s^{\prime}n_{i,s^{\prime}}$ is the local polarity at
site $i$.

  3. (3)

We repeat the procedures (1) and (2) $N$ times and increment time by $\Delta
t$.

We set $\Delta t=1/(8D_{0}+\gamma_{0})$ for the simulations of the model with
purely repulsive interaction (i.e., $U>0$ and $J=0$) (Figs. 3 and 5) and
$\Delta t=1/(4D_{0}+2\gamma_{0})$ for the simulations of the model with purely
aligning interaction (i.e., $U=0$ and $J>0$) (Figs. 4 and 6).

##  Appendix B Numerical sampling of configurations

In the lattice gas model simulations, we checked the relaxation of the long-
wavelength components of the structure factors, $S_{d}(\bm{k})$ and
$S_{p}(\bm{k})$ for $\bm{k}=(2\pi/L,0)$ and $\bm{k}=(0,2\pi/L)$, which is
slower than that of the short-wavelength components. In Figs. 7(a) and (b), we
plot the relaxation dynamics of these quantities averaged over independent
samples with the parameter sets used for Figs. 3 and 5 (purely repulsive
interaction), respectively. The gray area in each panel suggests the time
points used for averaging. Specifically, we took averages over $18432$
independent samples and $1800$ time points (taken every 100 MC steps after
relaxation with $2\times 10^{4}$ MC steps) to obtain the structure factors and
correlation functions for Fig. 3; $92160$ independent samples and $1800$ time
points (taken every 100 MC steps after relaxation with $2\times 10^{4}$ MC
steps) for Fig. 5.

In Figs. 8(a) and (b), we plot the corresponding relaxation dynamics with the
parameter sets used for Figs. 4 and 6 (purely aligning interaction),
respectively. We took averages over $44800$ independent samples and $800$ time
points (taken every 100 MC steps after relaxation with $4\times 10^{4}$ MC
steps) for Figs. 4; $4096$ independent samples and $50000$ time points (taken
every 100 MC steps after relaxation with $1\times 10^{6}$ MC steps) for Figs.
6.

## References

  * Marchetti  _et al._ [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143 (2013).
  * Needleman and Dogic [2017] D. Needleman and Z. Dogic, Active matter at the interface between materials science and cell biology, Nat. Rev. Mater. 2, 1 (2017).
  * Gompper  _et al._ [2020] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, T. Kiørboe, E. Lauga, W. C. K. Poon, A. DeSimone, S. Muiños-Landin, A. Fischer, N. A. Söker, F. Cichos, R. Kapral, P. Gaspard, M. Ripoll, F. Sagues, A. Doostmohammadi, J. M. Yeomans, I. S. Aranson, C. Bechinger, H. Stark, C. K. Hemelrijk, F. J. Nedelec, T. Sarkar, T. Aryaksama, M. Lacroix, G. Duclos, V. Yashunsky, P. Silberzan, M. Arroyo, and S. Kale, The 2020 motile active matter roadmap, J. Phys. Condens. Matter 32, 193001 (2020).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, I, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75, 1226 (1995).
  * Toner and Tu [1995] J. Toner and Y. Tu, Long-Range order in a Two-Dimensional dynamical XY model: How birds fly together, Phys. Rev. Lett. 75, 4326 (1995).
  * Toner [2012] J. Toner, Reanalysis of the hydrodynamic theory of fluid, polar-ordered flocks, Phys. Rev. E 86, 031918 (2012).
  * Chaté [2020] H. Chaté, Dry aligning dilute active matter, Annu. Rev. Condens. Matter Phys. 11, 189 (2020).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Motility-Induced phase separation, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Solon and Tailleur [2013] A. P. Solon and J. Tailleur, Revisiting the flocking transition using active spins, Phys. Rev. Lett. 111, 078101 (2013).
  * Solon and Tailleur [2015] A. P. Solon and J. Tailleur, Flocking with discrete symmetry: The two-dimensional active ising model, Phys. Rev. E 92, 042119 (2015).
  * Solon  _et al._ [2015] A. P. Solon, H. Chaté, and J. Tailleur, From phase to microphase separation in flocking models: the essential role of nonequilibrium fluctuations, Phys. Rev. Lett. 114, 068101 (2015).
  * Ngo  _et al._ [2014] S. Ngo, A. Peshkov, I. S. Aranson, E. Bertin, F. Ginelli, and H. Chaté, Large-scale chaos and fluctuations in active nematics, Phys. Rev. Lett. 113, 038302 (2014).
  * Garrido  _et al._ [1990] P. L. Garrido, J. L. Lebowitz, C. Maes, and H. Spohn, Long-range correlations for conservative dynamics, Phys. Rev. A 42, 1954 (1990).
  * Dorfman  _et al._ [1994] J. R. Dorfman, T. R. Kirkpatrick, and J. V. Sengers, Generic Long-Range correlations in molecular fluids, Annu. Rev. Phys. Chem. 45, 213 (1994).
  * Schmittmann and Zia [1995] B. Schmittmann and R. K. P. Zia, Statistical mechanics of driven diffusive systems, in  _Phase Transitions and Critical Phenomena_ , Vol. 17 (Academic Press, 1995).
  * Schmittmann and Zia [1998] B. Schmittmann and R. K. P. Zia, Driven diffusive systems. an introduction and recent developments, Phys. Rep. 301, 45 (1998).
  * Adachi  _et al._ [2022] K. Adachi, K. Takasan, and K. Kawaguchi, Activity-induced phase transition in a quantum many-body system, Phys. Rev. Res. 4, 013194 (2022).
  * Nakano and Adachi [2024] H. Nakano and K. Adachi, Universal properties of repulsive self-propelled particles and attractive driven particles, Phys. Rev. Res. 6, 013074 (2024).
  * Brambati  _et al._ [2022] M. Brambati, G. Fava, and F. Ginelli, Signatures of directed and spontaneous flocking, Phys Rev E 106, 024608 (2022).
  * Solon  _et al._ [2022] A. Solon, H. Chaté, J. Toner, and J. Tailleur, Susceptibility of polar flocks to spatial anisotropy, Phys. Rev. Lett. 128, 208004 (2022).
  * Chatterjee  _et al._ [2022] S. Chatterjee, M. Mangeat, and H. Rieger, Polar flocks with discretized directions: The active clock model approaching the vicsek model, EPL 138, 41001 (2022).
  * Mishra  _et al._ [2014] S. Mishra, S. Puri, and S. Ramaswamy, Aspects of the density field in an active nematic, Philos. Trans. R. Soc. A 372, 20130364 (2014).
  * Bröker  _et al._ [2023] S. Bröker, J. Bickmann, M. Te Vrugt, M. E. Cates, and R. Wittkowski, Orientation-Dependent propulsion of active brownian spheres: From Self-Advection to programmable cluster shapes, Phys. Rev. Lett. 131, 168203 (2023).
  * [24] S. Othman, J. Midya, T. Auth, and G. Gompper, Phase behavior and dynamics of active brownian particles in an alignment field, arXiv:2403.02947 .
  * Doi [1976] M. Doi, Second quantization representation for classical many-particle system, J. Phys. A 9, 1465 (1976).
  * Peliti [1985] L. Peliti, Path integral approach to birth-death processes on a lattice, J. Phys. France 46, 1469 (1985).
  * Täuber [2014] U. C. Täuber,  _Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior_ (Cambridge University Press, 2014).
  * Gwa and Spohn [1992] L. H. Gwa and H. Spohn, Six-vertex model, roughened surfaces, and an asymmetric spin hamiltonian, Phys. Rev. Lett. 68, 725 (1992).
  * Sandow [1994] S. Sandow, Partially asymmetric exclusion process with open boundaries, Phys. Rev. E 50, 2660 (1994).
  * Kim [1995] D. Kim, Bethe ansatz solution for crossover scaling functions of the asymmetric XXZ chain and the Kardar-Parisi-Zhang-type growth model, Phys. Rev. E 52, 3512 (1995).
  * Essler and Rittenberg [1996] F. H. L. Essler and V. Rittenberg, Representations of the quadratic algebra and partially asymmetric diffusion with open boundaries, J. Phys. A 29, 3375 (1996).
  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Non-Hermitian physics, Adv. Phys. 69, 249 (2020).
  * Murugan and Vaikuntanathan [2017] A. Murugan and S. Vaikuntanathan, Topologically protected modes in non-equilibrium stochastic systems, Nat. Commun. 8, 13881 (2017).
  * Dasbiswas  _et al._ [2018] K. Dasbiswas, K. K. Mandadapu, and S. Vaikuntanathan, Topological localization in out-of-equilibrium dissipative systems, Proc. Natl. Acad. Sci. U.S.A. 115, E9031 (2018).
  * Thompson  _et al._ [2011] A. G. Thompson, J. Tailleur, M. E. Cates, and R. A. Blythe, Lattice models of nonequilibrium bacterial dynamics, J. Stat. Mech. 2011, P02029 (2011).
  * Peruani  _et al._ [2011] F. Peruani, T. Klauss, A. Deutsch, and A. Voss-Boehme, Traffic jams, gliders, and bands in the quest for collective motion of self-propelled particles, Phys. Rev. Lett. 106, 128101 (2011).
  * Whitelam  _et al._ [2018] S. Whitelam, K. Klymko, and D. Mandal, Phase separation and large deviations of lattice active matter, J. Chem. Phys. 148, 154902 (2018).
  * Kourbane-Houssene  _et al._ [2018] M. Kourbane-Houssene, C. Erignoux, T. Bodineau, and J. Tailleur, Exact hydrodynamic description of active lattice gases, Phys. Rev. Lett. 120, 268003 (2018).
  * Partridge and Lee [2019] B. Partridge and C. F. Lee, Critical Motility-Induced phase separation belongs to the ising universality class, Phys. Rev. Lett. 123, 068002 (2019).
  * [40] R. Mukherjee, S. Saha, T. Sadhu, A. Dhar, and S. Sabhapandit, Hydrodynamics of a hard-core non-polar active lattice gas, arXiv:2405.19984 .
  * Note [1] Equation (1) is a simpler version of Eq. (2) in Ref. [17] [or Eq. (S-21) in Ref. [40]], which is derived for a repulsive active lattice gas model.
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49, 435 (1977).
  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Speck  _et al._ [2015] T. Speck, A. M. Menzel, J. Bialké, and H. Löwen, Dynamical mean-field theory and weakly non-linear analysis for the phase separation of active brownian particles, J. Chem. Phys. 142, 224109 (2015).
  * Sternheim and Walker [1972] M. M. Sternheim and J. F. Walker, Non-Hermitian hamiltonians, decaying states, and perturbation theory, Phys. Rev. C 6, 114 (1972).
  * Sasa [2005] S.-I. Sasa, Long range spatial correlation between two brownian particles under external driving, Physica D 205, 233 (2005).
  * Lefevere and Tasaki [2005] R. Lefevere and H. Tasaki, High-temperature expansion for nonequilibrium steady states in driven lattice gases, Phys. Rev. Lett. 94, 200601 (2005).
  * Luo  _et al._ [2023] Y. Luo, M. Gu, M. Park, X. Fang, Y. Kwon, J. M. Urueña, J. Read de Alaniz, M. E. Helgeson, C. M. Marchetti, and M. T. Valentine, Molecular-scale substrate anisotropy, crowding and division drive collective behaviours in cell monolayers, J. R. Soc. Interface 20, 20230160 (2023).
  * Gu  _et al._ [2023] M. Gu, X. Fang, and Y. Luo, Data-Driven model construction for anisotropic dynamics of active matter, PRX Life 1, 013009 (2023).
  * Skillin  _et al._ [2023] N. P. Skillin, B. E. Kirkpatrick, K. M. Herbert, B. R. Nelson, G. K. Hach, K. A. Günay, R. M. Khan, F. W. DelRio, T. J. White, and K. S. Anseth, Stiffness anisotropy coordinates supracellular contractility driving long-range myotube-ecm alignment, bioRxiv 10.1101/2023.08.08.552197 (2023).


realization due to $d\geq 0$. Accordingly,
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is internally positive
and ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\in{\mathscr{P}}$,
i.e., ${\mathscr{P}}_{\\!{{\rho}}}\subset{\mathscr{P}}$. From this result and
the definition of ${\mathscr{P}}_{\\!(0,1)}$, the last claim is directly
implied.

###  A-B Proof of Theorem 7

Let $\varepsilon>0$ and ${\mathrm{g}}\in{\mathscr{P}}$ with transfer function
$G$. Since ${\mathscr{P}}\subset\ell^{1}$, each element of ${\mathscr{P}}$ is
a stable impulse response, and therefore, we have $r(G)<1$. Let ${{\rho}}$ and
${a}$ be positive real scalars such that ${{\rho}}\in(r(G),1)$ and
${a}<(1-{{\rho}})\varepsilon$. Consider an impulse response
${\mathrm{g}}^{{\scalebox{0.63}{($\varepsilon$)}}}=(g^{{\scalebox{0.63}{($\varepsilon$)}}}_{t})_{t=0}^{\infty}$
with transfer function $G^{{\scalebox{0.63}{($\varepsilon$)}}}$ where, for any
$t\in{\mathbb{Z}}_{+}$, $g_{t}$ is defined as
$g^{{\scalebox{0.63}{($\varepsilon$)}}}_{t}=g_{t}+{a}{{\rho}}^{t}$. For any
$t\in{\mathbb{Z}}_{+}$, one has $g_{t}\geq 0$, and since ${a},{{\rho}}>0$, it
follows that $g^{{\scalebox{0.63}{($\varepsilon$)}}}_{t}\geq 0$, i.e.,
${\mathrm{g}}^{{\scalebox{0.63}{($\varepsilon$)}}}$ is a non-negative impulse
response. Moreover, one can easily see that

$G^{{\scalebox{0.63}{($\varepsilon$)}}}(z)=G(z)+\frac{{a}}{1-{{\rho}}z^{-1}}.$
(109)

Since ${\mathrm{g}}$ belongs to ${\mathscr{P}}$, we know that $G$ is a
rational function. Accordingly, due to (109) and Theorem 2, it follows that
(4) holds for ${\mathrm{g}}^{{\scalebox{0.63}{($\varepsilon$)}}}$. Moreover,
${{\rho}}>r(G)$ implies that $\lim_{t\to\infty}{{\rho}}^{-t}g_{t}=0$.
Subsequently, one has
$\lim_{t\to\infty}{{\rho}}^{-t}g^{{\scalebox{0.63}{($\varepsilon$)}}}_{t}={a}$,
and therefore,
${\mathrm{g}}^{{\scalebox{0.63}{($\varepsilon$)}}}\in{\mathscr{P}}_{\\!{{\rho}}}\subset{\mathscr{P}}_{\\!(0,1)}$.
For the case of $p=\infty$, we have

$\|{\mathrm{g}}-{\mathrm{g}}^{{\scalebox{0.63}{($\varepsilon$)}}}\|_{\infty}=\sup_{t\in{\mathbb{Z}}_{+}}{a}{{\rho}}^{t}={a}<(1-{{\rho}})\varepsilon<\varepsilon.$

On the other hand, for $p\in[1,\infty)$, one can see that

$\|{\mathrm{g}}-{\mathrm{g}}^{{\scalebox{0.63}{($\varepsilon$)}}}\|_{p}={a}\Big{(}\sum_{t=0}^{\infty}{{\rho}}^{pt}\Big{)}^{\frac{1}{p}}\\!\\!=\frac{{a}}{(1-{{\rho}}^{p})^{\frac{1}{p}}}<\frac{(1-{{\rho}})\varepsilon}{(1-{{\rho}}^{p})^{\frac{1}{p}}}\leq\varepsilon,$

where the last inequality is due to ${{\rho}}^{p}+(1-{{\rho}})^{p}\leq 1$
which holds for any ${{\rho}}\in(0,1)$ and $p\in[1,\infty)$. ∎

###  A-C Proof of Theorem 17

Let $c_{00}$ be the space of impulse responses which are finitely non-zero,
i.e., for each ${\mathrm{g}}=(g_{s})_{s=0}^{\infty}$ there exists
$n_{{\mathrm{g}}}\in{\mathbb{Z}}_{+}$ such that $g_{s}=0$, for all $s\geq
n_{{\mathrm{g}}}$. Note that for such ${\mathrm{g}}\in c_{00}$, we have
$\mathrm{Hankel}({\mathrm{g}}){\mathrm{v}}\in{\mathbb{R}}^{n_{{\mathrm{g}}}}\times\\{\mathbf{0}\\}$,
for any ${\mathrm{v}}\in\ell^{\infty}$. This implies that
$\mathrm{rank}\big{(}\mathrm{Hankel}({\mathrm{g}})\big{)}\leq
n_{{\mathrm{g}}}<\infty$, and, therefore, ${\mathrm{g}}$ is a finite Hankel
rank impulse response. Let
$\mathbbm{k}:{\mathbb{Z}}_{+}\times{\mathbb{Z}}_{+}\to{\mathbb{R}}$ be a
finite support kernel, i.e., there exists $n_{\mathbbm{k}}\in{\mathbb{Z}}_{+}$
such that $\mathbbm{k}(s,t)=0$ when $s\geq n_{\mathbbm{k}}$ or $t\geq
n_{\mathbbm{k}}$. One can easily see that $\mathbbm{k}_{t}\in c_{00}$, for any
$t\in{\mathbb{Z}}_{+}$. Therefore, we have
$\mathrm{rank}\big{(}\mathrm{Hankel}(\mathbbm{k}_{t})\big{)}\leq
n_{\mathbbm{k}}<\infty$, and consequently, $\mathbbm{k}$ is a finite Hankel
rank kernel.

Let ${\mathrm{f}}_{\beta}=(f_{t})_{t=0}^{\infty}$ be an impulse response
defined as $f_{t}=\beta^{t}$, for $t\in{\mathbb{Z}}_{+}$. One can easily see
that

$\\{\mathrm{Hankel}({\mathrm{f}}_{\beta}){\mathrm{v}}\ \\!\big{|}\
\\!{\mathrm{v}}\in\ell^{\infty}\Big{\\}}=\\{{\mathrm{f}}_{\beta}v\ \\!\big{|}\
\\!v\in{\mathbb{R}}\\}.$ (110)

Therefore, we have
$\mathrm{rank}\big{(}\mathrm{Hankel}({\mathrm{f}}_{\beta})\big{)}=1$, and
consequently, ${\mathrm{f}}_{\beta}$ is a finite Hankel rank impulse response.
For the TC kernel introduced in (70) and $t\in{\mathbb{Z}}_{+}$, consider the
section of $\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$ at $t$,
i.e.,
$(\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}},t})_{s=0}^{\infty}$.
Note that we have
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}},t}={\mathrm{f}}_{\beta}+{\mathrm{g}}$,
where impulse response ${\mathrm{g}}=(g_{s})_{s=0}^{\infty}$ is defined by
$g_{s}=\beta^{\max(s,t)}-\beta^{s}$, for $s\in{\mathbb{Z}}_{+}$. One can
easily see that $g_{s}=0$, for all $s\geq t$, which implies that
${\mathrm{g}}\in c_{00}$, and subsequently,
$\mathrm{rank}\big{(}\mathrm{Hankel}({\mathrm{g}})\big{)}<\infty$.
Accordingly, since rank of $\mathrm{Hankel}({\mathrm{f}}_{\beta})\big{)}$ is
finite, it follows that
$\mathrm{rank}\big{(}\mathrm{Hankel}(\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}},t})\big{)}<\infty$,
and consequently, $\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$ is a
finite Hankel rank kernel. Based on similar arguments, one can show same
result for $\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}$ and
$\mathbbm{k}_{\SS}$, concluding the proof. ∎

## References

  * [1] W. M. Haddad, V. Chellaboina, and Q. Hui, _Nonnegative and compartmental dynamical systems_. Princeton University Press, 2010.
  * [2] R. F. Brown, “Compartmental system analysis: State of the art,” _IEEE Transactions on Biomedical Engineering_ , no. 1, pp. 1–11, 1980.
  * [3] R. Shorten, F. Wirth, and D. Leith, “A positive systems model of TCP-like congestion control: Asymptotic results,” _IEEE/ACM Transactions on Networking_ , vol. 14, no. 3, pp. 616–629, 2006.
  * [4] U. Krause, _Positive Dynamical Systems in Discrete Time_. de Gruyter, 2015.
  * [5] L. Farina and S. Rinaldi, _Positive Linear Systems: Theory and Applications_. John Wiley & Sons, 2011\.
  * [6] A. Rantzer and M. E. Valcher, “A tutorial on positive systems and large scale control,” in _Conference on Decision and Control_ , 2018, pp. 3686–3697.
  * [7] L. Benvenuti and L. Farina, “A tutorial on the positive realization problem,” _IEEE Transactions on Automatic Control_ , vol. 49, no. 5, pp. 651–664, 2004.
  * [8] D. G. Luenberger, “Introduction to Dynamic Systems; Theory, Models, and Applications,” Tech. Rep., 1979.
  * [9] M. E. Valcher, “Reachability properties of continuous-time positive systems,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 7, pp. 1586–1590, 2009.
  * [10] J. Back and A. Astolfi, “Design of positive linear observers for positive linear systems via coordinate transformations and positive realizations,” _SIAM Journal on Control and Optimization_ , vol. 47, no. 1, pp. 345–373, 2008.
  * [11] D. Hinrichsen and E. Plischke, “Robust stability and transient behaviour of positive linear systems,” _Vietnam Journal of Mathematics_ , vol. 35, no. 4, pp. 429–462, 2007.
  * [12] M. Colombino and R. S. Smith, “A convex characterization of robust stability for positive and positively dominated linear systems,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 7, pp. 1965–1971, 2015.
  * [13] B. Roszak and E. J. Davison, “Necessary and sufficient conditions for stabilizability of positive lti systems,” _Systems and Control Letters_ , vol. 58, no. 7, pp. 474–481, 2009.
  * [14] B. Shafai, J. Chen, and M. Kothandaraman, “Explicit formulas for stability radii of nonnegative and Metzlerian matrices,” _IEEE Transactions on Automatic Control_ , vol. 42, no. 2, pp. 265–270, 1997.
  * [15] A. Oghbaee, B. Shafai, and S. Nazari, “Complete characterisation of disturbance estimation and fault detection for positive systems,” _IET Control Theory and Applications_ , vol. 12, no. 7, pp. 883–891, 2018.
  * [16] N. K. Dhingra, M. Colombino, and M. R. Jovanović, “Structured decentralized control of positive systems with applications to combination drug therapy and leader selection in directed networks,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 1, pp. 352–362, 2018\.
  * [17] Y. Ebihara, D. Peaucelle, and D. Arzelier, “Decentralized control of interconnected positive systems using l 1-induced norm characterization,” in _2012 IEEE 51st IEEE Conference on Decision and Control (CDC)_. IEEE, 2012, pp. 6653–6658.
  * [18] Y. Ebihara, D. Peaucelle, and D. Arzelier, “Stability and persistence analysis of large scale interconnected positive systems,” in _2013 European Control Conference (ECC)_. IEEE, 2013, pp. 3366–3371.
  * [19] A. Rantzer and M. E. Valcher, “Scalable control of positive systems,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 4, pp. 319–341, 2021.
  * [20] A. A. Ahmadi and B. El Khadir, “Learning dynamical systems with side information (short version),” _Proceedings of Machine Learning Research_ , vol. 120, pp. 718–727, 2020.
  * [21] M. Khosravi and R. S. Smith, “Kernel-based identification of positive systems,” in _Conference on Decision and Control_ , 2019, pp. 1740–1745.
  * [22] M. Khosravi and R. S. Smith, “Nonlinear system identification with prior knowledge on the region of attraction,” _IEEE Control Systems Letters_ , vol. 5, no. 3, pp. 1091–1096, 2021.
  * [23] M. Khosravi and R. S. Smith, “Convex nonparametric formulation for identification of gradient flows,” _IEEE Control Systems Letters_ , vol. 5, no. 3, pp. 1097–1102, 2021.
  * [24] J. Umenberger and I. R. Manchester, “Scalable identification of stable positive systems,” in _Conference on Decision and Control_. IEEE, 2016, pp. 4630–4635.
  * [25] L. Benvenuti, A. De Santis, and L. Farina, “On model consistency in compartmental systems identification,” _Automatica_ , vol. 38, no. 11, pp. 1969–1976, 2002.
  * [26] A. De Santis and L. Farina, “Identification of positive linear systems with Poisson output transformation,” _Automatica_ , vol. 38, no. 5, pp. 861–868, 2002.
  * [27] C. Grussler, J. Umenberger, and I. R. Manchester, “Identification of externally positive systems,” in _Conference on Decision and Control_. IEEE, 2017, pp. 6549–6554.
  * [28] M. Zheng and Y. Ohta, “Bayesian positive system identification: Truncated Gaussian prior and hyperparameter estimation,” _Systems and Control Letters_ , vol. 148, p. 104857, 2021.
  * [29] M. Zheng and Y. Ohta, “Positive FIR system identification using maximum entropy prior,” _IFAC-PapersOnLine_ , vol. 51, no. 15, pp. 7–12, 2018.
  * [30] G. Pillonetto and G. De Nicolao, “A new kernel-based approach for linear system identification,” _Automatica_ , vol. 46, no. 1, pp. 81–93, 2010.
  * [31] L. Ljung, T. Chen, and B. Mu, “A shift in paradigm for system identification,” _International Journal of Control_ , vol. 93, no. 2, pp. 173–180, 2020.
  * [32] M. Khosravi and R. S. Smith, “On robustness of kernel-based regularized system identification,” _IFAC-PapersOnLine_ , vol. 54, no. 7, pp. 749–754, 2021, IFAC Symposium on System Identification.
  * [33] G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, and L. Ljung, “Kernel methods in system identification, machine learning and function estimation: A survey,” _Automatica_ , vol. 50, no. 3, pp. 657–682, 2014.
  * [34] M. Khosravi, A. Iannelli, M. Yin, A. Parsi, and R. S. Smith, “Regularized system identification: A hierarchical Bayesian approach,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 406–411, 2020, IFAC World Congress 2020.
  * [35] Y. Fujimoto, I. Maruta, and T. Sugie, “Extension of first-order stable spline kernel to encode relative degree,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 14 016–14 021, 2017.
  * [36] T. Chen, H. Ohlsson, and L. Ljung, “On the estimation of transfer functions, regularizations and Gaussian processe – Revisited,” _Automatica_ , vol. 48, no. 8, pp. 1525–1535, 2012.
  * [37] M. Khosravi, M. Yin, A. Iannelli, A. Parsi, and R. S. Smith, “Low-complexity identification by sparse hyperparameter estimation,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 412–417, 2020, IFAC World Congress 2020.
  * [38] M. A. H. Darwish, G. Pillonetto, and R. Tóth, “The quest for the right kernel in Bayesian impulse response identification: The use of OBFs,” _Automatica_ , vol. 87, pp. 318–329, 2018.
  * [39] T. Chen, “On kernel design for regularized LTI system identification,” _Automatica_ , vol. 90, pp. 109–122, 2018.
  * [40] A. Marconato, M. Schoukens, and J. Schoukens, “Filter-based regularisation for impulse response modelling,” _IET Control Theory and Applications_ , vol. 11, no. 2, pp. 194–204, 2016.
  * [41] M. Khosravi and R. S. Smith, “Kernel-based impulse response identification with side-information on steady-state gain,” _arXiv preprint_ , 2021.
  * [42] R. S. Risuleo, G. Bottegal, and H. Hjalmarsson, “A nonparametric kernel-based approach to Hammerstein system identification,” _Automatica_ , vol. 85, pp. 234–247, 2017.
  * [43] R. S. Risuleo, F. Lindsten, and H. Hjalmarsson, “Bayesian nonparametric identification of Wiener systems,” _Automatica_ , vol. 108, p. 108480, 2019\.
  * [44] N. Everitt, G. Bottegal, and H. Hjalmarsson, “An empirical Bayes approach to identification of modules in dynamic networks,” _Automatica_ , vol. 91, pp. 144–151, 2018.
  * [45] M. Khosravi and R. S. Smith, “Kernel-based identification with frequency domain side-information,” _arXiv preprint_ , 2021.
  * [46] P. A. Fuhrmann, _A polynomial approach to linear algebra_. Springer Science and Business Media, 2011.
  * [47] A. Berlinet and C. Thomas-Agnan, _Reproducing Kernel Hilbert Spaces in Probability and Statistics_. Springer Science and Business Media, 2011.
  * [48] T. Chen and G. Pillonetto, “On the stability of reproducing kernel Hilbert spaces of discrete-time impulse responses,” _Automatica_ , vol. 95, pp. 529–533, 2018.
  * [49] J. Peypouquet, _Convex Optimization in Normed Spaces: Theory, Methods and Examples_. Springer, 2015\.
  * [50] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” in _International Conference on Computational Learning Theory_. Springer, 2001, pp. 416–426.
  * [51] F. Dinuzzo and B. Schölkopf, “The representer theorem for Hilbert spaces: A necessary and sufficient condition,” in _Advances in neural information processing systems_ , 2012, pp. 189–196.
  * [52] F. Dinuzzo, “Kernels for linear time invariant system identification,” _SIAM Journal on Control and Optimization_ , vol. 53, no. 5, pp. 3299–3317, 2015.
  * [53] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” 2014.
  * [54] J. Gondzio, P. González-Brevis, and P. Munari, “Large-scale optimization with the primal-dual column generation method,” _Mathematical Programming Computation_ , vol. 8, no. 1, pp. 47–82, 2016.
  * [55] N. Srinivas, A. Krause, S. M. Kakade, and M. W. Seeger, “Information-theoretic regret bounds for Gaussian process optimization in the bandit setting,” _IEEE Transactions on Information Theory_ , vol. 58, no. 5, pp. 3250–3265, 2012.
  * [56] G. Dullerud and R. Smith, “Sampled-data model validation: An algorithm and experimental application,” _International Journal of Robust and Nonlinear Control_ , vol. 6, no. 9-10, pp. 1065–1078, 1996.
  * [57] B. De Moor, P. De Gersem, B. De Schutter, and W. Favoreel, “Daisy: A database for identification of systems,” _JOURNAL A_ , vol. 38, pp. 4–5, 1997.