good candidates to act as order parameters. Note that, in dimension $d=2$ and
$d=3$, it can even be shown that they are the only ones, in the sense that any
irreducible representation of $SO(d)$ is isomorphic to one of the
$\mathcal{H}_{n}$.

In practice, finding an order parameter that accounts for a spontaneous
symmetry breaking $SO(d)\to G$ can be done as follows: among all the
$f_{n}$’s, with $n\neq 0$, such that $\mathcal{H}_{n}$ contains a (non zero)
$G$–invariant subspace, choose the one with the smallest index $n_{0}$. The
component $f_{n_{0}}$ obtained then satisfies:

  * $\rightarrow$

condition (i) since $n_{0}\neq 0$, and any rotation–invariant function on
$\mathbb{S}^{d-1}$ has all its component $f_{n}=0$, for all $n>0$;

  * $\rightarrow$

condition (ii) since, if $f$ is $G$-invariant, so are all the $f_{n}$, and
$f_{n_{0}}$ in particular. If it is the case, then $f_{n_{0}}$ belongs to the
subspace of $\mathcal{H}_{n}$ that is $G$–invariant. Note that, in general,
this component can vanish for very peculiar $f$, but is expected to be non-
zero in a generic case. For instance, if $d=2$ and $G$ contains the identity
and the $180^{\circ}$–rotation only, then $n_{0}=2$. If $f=\cos(4\theta)$,
then $f$ is indeed $G-$invariant while its component on $\mathcal{H}_{2}$
vanishes, but its full symmetry group is actually bigger since it contains
also $90$ and $270^{\circ}$ rotations.

  * $\rightarrow$

the linear–transformation criterion of (iii). In addition, by choosing the
smallest $n_{0}$ satisfying the aforementioned properties, we chose the
simplest possible parameter, in a sense that will become clear in the next
section.

In practice, an order parameter obtained by the procedure described above is
not very convenient yet. Indeed, to describe the dynamics of an
inhomogeneous—e.g nematic or polar—system, it is not very convenient to model
its evolution by means of a field that assigns to each point in space an
element of $\mathcal{H}_{n_{0}}$, i.e. a $\mathbb{R}$–valued function on the
unit sphere. We would usually rather like to be able to locally describe the
system by arrays of numbers. In addition, these arrays of numbers should
transform nicely under change of coordinates. To do so, one could try to pick
an arbitrary basis of $\mathcal{H}_{n_{0}}$, and stack the coordinates of the
local $\mathcal{H}_{n_{0}}$–value of the field into a column vector.
Unfortunately, the resulting array would generically transform in an awful way
under change of basis. Luckily, there exists a slightly more subtle,
basis–invariant, way of associating to any element of $\mathcal{H}_{n_{0}}$ a
“nice array of numbers”, i.e. a tensor, as described in the following
paragraph.

##### From harmonic scalar functions to harmonic tensors.

At this stage, the properties of decomposition (A.1) haven’t yet been
exploited to the full extent of their possibility. Indeed, any function
belonging to a given $\mathcal{H}_{n}$ turns out to be the restriction of a
homogeneous, harmonic polynomial444We recall that a polynomial on
$\mathbb{R}^{d}$ is said to be harmonic whenever it lies in the kernel of the
Laplacian $\Delta$ of $\mathbb{R}^{d}$. of order $n$ on $\mathbb{R}^{d}$. In
particular, for each element $Y_{n\ell}(\mathbb{u})$ of the arbitrarily–chosen
basis appearing in Eq. (A.2), we denote this polynomial on $\mathbb{R}^{d}$ by
$Y_{n\ell}(\mathbb{r})$. Then, we recall that for, any $n^{\rm th}$–order
homogeneous polynomial on $\mathbb{R}^{d}$, there exists a unique $n^{\rm
th}$–order symmetric tensor $\mathbb{Y}_{n\ell}$—determined by
polarization—such that

$Y_{n\ell}(\mathbb{u})=\mathbb{Y}_{n\ell}\cdot\mathbb{u}^{\otimes
n}=\sum_{i_{1},\dots,i_{n}=1}^{d}Y^{i_{1}\dots i_{n}}u_{i_{1}}\dots u_{i_{n}}\
.$ (A.3)

Furthermore, the harmonicity of $Y_{n\ell}(\mathbb{r})$ is equivalent to the
traceless–ness of the associated tensor:
$\mathrm{tr}(\mathbb{Y}_{n\ell})\equiv\sum_{j}Y^{jji_{3}\dots i_{n}}=0$—hence
the name “harmonic tensors” assigned to symmetric traceless tensors.

Each $f_{n}$ in decomposition (A.2) being itself a homogeneous polynomial of
order $n$, there exist a unique corresponding harmonic tensor $\mathbb{A}_{n}$
such that

$f(\mathbb{u})=\sum_{n\in\mathbb{N}}\mathbb{A}_{n}\cdot\mathbb{u}^{\otimes
n}=\sum_{n\in\mathbb{N}}\left(\sum_{\ell=1}^{\dim\mathcal{H}_{n}}c_{n\ell}\mathbb{Y}_{n\ell}\right)\cdot\mathbb{u}^{\otimes
n}\ .$ (A.4)

It can then be shown that

$\mathbb{A}_{n}=\frac{1}{\Omega}\frac{(2n+d-2)!!}{n!!(d-2)!!}\int_{\mathbb{S}_{d-1}}f(\mathbb{u})\widehat{\mathbb{u}^{\otimes
n}}\;\mathrm{d}\mathbb{u}\ ,$ (A.5)

which is, together with (A.4) and up to a rescaling of the $\mathbb{A}_{n}$,
the decomposition described in the section 2.2.1 of the main text.

In addition to be very convenient for computing the diffusive limits ABP and
RTP in a coordinate–free manner as done in this article, the tensors
$\mathbb{A}_{n}$ are, up to a rescaling (and potential rewriting in complex
form for the case $d=2$), the usual order parameters chosen to describe phase
transitions accompanied by a spontaneous symmetry breaking as
$SO(d)\rightarrow G\subset SO(d)$ (i.e. isotropic to e.g. polar, nematic,
hexatic,$\dots$ phases).

##  Appendix B Numerical details

###  B.1 Microscopic simulations

Microscopic simulations of RTPs, ABPs and AOUPs are carried out in $d=1$ and
$d=2$ dimensions in continuous space with periodic boundary conditions. For
non-interacting simulations, at each time-step $dt$ we first compute the
space-dependent motility parameters $\\{\gamma({\bf r}_{i})\\}$ for each
particle $i$. We then update the position ${\bf r}_{i}$ by an amount
$v_{i}{\bf u}_{i}dt+\sqrt{2D_{t}dt}\boldsymbol{\Delta\eta}_{i}$, where
$\boldsymbol{\Delta\eta}_{i}$ is a vector of $d$ independent Gaussian random
variables of unit variance and zero mean. After displacing the particle we
update its orientation vector ${\bf u}_{i}$ depending on the specific
dynamics:

  * •

For RTPs, we draw a random number $\delta$ uniformly distributed on $[0,1]$
and compare it with $\varepsilon=\alpha({\bf r}_{i})dt$. If
$\delta<\varepsilon$, we draw a new orientation ${\bf u}_{i}^{\prime}$ from
the unit sphere $\mathbb{S}^{d-1}$, otherwise we keep the same orientation
${\bf u}_{i}(t+dt)={\bf u}_{i}(t)$. In $1d$, this operation simply corresponds
to flipping the orientation $\sigma\in\\{\pm 1\\}$ with probability
$\varepsilon=\frac{1}{2}\alpha(x_{i})dt$.

  * •

For ABPs in $2d$, we draw a random number from a Gaussian distribution
$\mathcal{N}(0,1)$ and update the orientation angle as
$\theta_{i}(t+dt)=\theta_{i}(t)+\sqrt{2\Gamma({\bf
r}_{i})dt}\mathcal{N}(0,1)$, with ${\bf
u}_{i}=(\cos\theta_{i},\sin\theta_{i})$.

  * •

For AOUPs, we update the orientation by ${\bf e}_{i}(t+dt)={\bf
e}_{i}+\sqrt{2\omega({\bf r}_{i})dt}\boldsymbol{\Delta\xi}_{i}$, where
$\boldsymbol{\Delta\xi}_{i}$ is a $d$-dimensional vector of independent,
centered Gaussian random variables with unit variance.

For simulations of QS-RTPs in $1d$ we resort to spatial hashing: we divide the
simulation domain in squared boxes of linear size $r_{QS}=1$, so we can
compute the local density around each particle by looking at the box where it
belongs and the neighboring boxes using Eq. (3.40). Spatial hashing reduces
the computational complexity to $\mathcal{O}(N)$. Based on the local density,
we update the self-propulsion speed $v(\rho)$ every $dt$ according to Eq.
(3.39). The tumbling events are instead simulated as follows: at the $n$-th
tumbling-event occurring at time $t_{n}$, we draw a waiting time $\delta t$
from an exponential distribution $p(\delta
t)\propto\exp(\frac{\alpha_{0}}{2}\delta t)$, so that $t_{n+1}=t_{n}+\delta
t$. Between $t_{n}$ and $t_{n+1}$ the particle maintains the same orientation
$\sigma(t_{n})$, while at $t_{n+1}$ it flips, so that
$\sigma(t_{n+1})=-\sigma(t_{n})$.

###  B.2 Structure factor

To measure the static structure factor $S(q)$ in $1d$ simulations, we first
define the interval of Fourier modes $q_{n}=\frac{2\pi}{L_{x}}n$, with
$n\in\left\\{1,\dots,N_{q}-1\right\\}$, where $L_{x}$ is the domain size.
After waiting a relaxation time $\tau_{\rm relax}$, we start measuring the
structure factor at intervals of $d\tau_{S_{Q}}$. For each wavevector $q_{n}$
we compute the associated Fourier component of the density $\hat{\rho}_{n}(t)$
as:

$\hat{\rho}_{n}(t)=\sum_{i=1}^{N}e^{-iq_{n}x_{i}(t)}\;,$ (B.1)

where $i$ is the particle index and $x_{i}$ the particle’s position. The
structure factor sampled at time $t$ is then given by:

$S(q_{n},t)=\frac{1}{N}|\hat{\rho}_{n}(t)|^{2}\delta_{n,0}\;.$ (B.2)

Note that $\delta_{n,0}$ ensures mass conservation, hence Eq. (B.2) coincides
with the definition (3.27) of $S(q)$ given in the main text. The final curve
for $S(q_{n})$ is obtained by averaging over all samples at different times
$t$: $S(q_{n})=\langle S(q_{n},t)\rangle_{t}\;$.

###  B.3 Spatial correlation

To compute the theoretical prediction for the spatial correlation function
$G(r)$ we need to (inverse) Fourier transform $S(q_{n})$ according to Eq.
(3.34). We remind that $S(q)\to 1$ for $q\to\infty$, hence the need of
shifting it by a constant $-1$ before applying the Fourier transform.
To measure the density-density correlation $G(r)$ in $1d$ simulations, we
first measure the pair distribution function $g(r)$, which is related to the
$G(r)$ in a homogeneous and isotropic system. Let $\rho_{0}$ be the average
particle density of the system. Following [77]:

$\displaystyle G(r)$ $\displaystyle=$
$\displaystyle\langle\delta\rho(r)\delta\rho(0)\rangle=\langle\rho(r)\rho(0)\rangle-\rho_{0}^{2}$
(B.3) $\displaystyle=$ $\displaystyle\left\langle\sum_{i=1}^{N}\sum_{j\neq
i}^{N}\delta(r-x_{i})\delta(x_{j})\right\rangle+\left\langle\sum_{i=1}^{N}\delta(r-x_{i})\delta(x_{i})\right\rangle-\rho_{0}^{2}$
$\displaystyle=$
$\displaystyle\rho_{0}^{2}g(r)+\delta(r)\left\langle\sum_{i=1}^{N}\delta(r-x_{i})\right\rangle-\rho_{0}^{2}=\rho_{0}^{2}\left[g(r)-1\right]+\rho_{0}\delta(r)\;.$

The function $g(r)$ gives the probability of finding a particle at distance
$[r,r+dr]$ from a particle located at the origin. In simulations, we proceed
as follows:

  * •

Choose an interval $[r_{\rm min},r_{\rm max}]$ over which to sample $g(r)$.
This interval is divided into bins with spacing $dr$.

  * •

Choose a subset of $K=N/10$ reference particles.

  * •

Every $dt_{\rm meas}$:

    * –

For each reference particle $k$, we look at its neighbors and measure
$n_{i}^{(k)}$, the number of particles at distance $\Delta$, with $r_{\rm
min}+i\>dr\leq\Delta\leq r_{\rm min}+(i+1)dr$. The $i^{th}$ bin considers both
particles to the left and right of the reference one.

    * –

We average $n_{i}^{(k)}$ over the $K$ reference particles:
$n_{i}=\frac{1}{K}\sum_{k=1}^{K}n_{i}^{(k)}$.

    * –

In $1d$, the pair distribution function $g_{i}$ for the $i^{th}$ spatial bin
can be computed as [77]:

$g_{i}=\frac{1}{\rho_{0}}\frac{n_{i}}{2dr}\;,$ (B.4)

corresponding to the ratio between the density at distance $r_{i}$ from a
reference particle and the homogeneous density $\rho_{0}$.

  * •

We average $g_{i}$ over the different time measurements to obtain the final
curve for $g$. Using Eq. (B.3) we compute the associated correlation function
$G(r)$.

###  B.4 Intermediate scattering function

To measure the intermediate scattering function $F(q,t)$ in our simulations,
we first choose a subset of $n$ Fourier modes $q_{n}\in\\{q_{\rm min},q_{\rm
min}+dq,\dots,q_{\rm min}+(n-1)dq\\}$. After waiting a relaxation time
$\tau_{\rm relax}$ we start measuring $F(q,t)$ for each mode $q_{n}$ over a
time window $[\tau_{i},\tau_{i}+t_{\rm max}]$, with a time separation $dt_{\rm
meas}$ between each measurement inside this interval. The quantity we measure
is:

$F_{i}(q_{n},t)=\hat{\rho}_{n}(\tau_{i})\hat{\rho}_{n}(\tau_{i}+t)\;\;,\quad
t\in\\{0,dt_{\rm meas},2dt_{\rm meas},\dots,t_{\rm max}\\}\;,$ (B.5)

where $\hat{\rho}_{n}$ is defined by Eq. (B.1). After storing a matrix of
$F_{i}(q_{n},t)$ for the $i^{th}$ time window, we start a new measurement at
$\tau_{i+1}=\tau_{i}+\mathrm{d}\tau$. In this way, we eventually collect $K$
samples of $F_{i}(q,t)$ over different time windows. The final intermediate
scattering function is given by averaging $F_{i}(q_{n},t)$ over the different
samples $i$:

$F(q_{n},t)=\frac{1}{K}\sum_{i=1}^{K}F_{i}(q_{n},t)\;.$ (B.6)

##  Appendix C Stationary distributions obtained from coarse-grained equations

In this Appendix we show how the coarse-grained equations (2.59), (2.88) can
be used to compute the steady-state distribution of active particles with
space-dependent motility parameters. As explained in Sec. 2.5 of the main
text, the scope of validity of these solutions is limited to spatial
modulations that occur over large spatial scales compared to the persistence
length $\ell_{p}$.

###  C.1 ABP-RTPs with position-dependent speed and translational noise

We consider the case of ABP-RTP with:

$\displaystyle v=v_{0}({\bf
r})\;,\quad\Gamma=\Gamma_{0}\;,\quad\alpha=\alpha_{0}$

and translational noise $D_{t}$. In the limit case where $D_{t}=0$ the
associated master equation (2.44) is exactly solvable at steady-state, and the
corresponding solution reads [32]:

$\mathcal{P}_{s}({\bf r},{\bf u})\propto\frac{1}{v_{0}({\bf r})}\;.$ (C.1)

Whenever $D_{t}>0$, instead, no exact solution for the microscopic dynamics
(4.5) is known for generic $v_{0}({\bf r})$. However, we can look for an
approximate solution $\tilde{p}_{s}({\bf r})$ for the coarse-grained dynamics:

$\partial_{t}\tilde{p}_{s}({\bf r})=-\nabla_{{\bf
r}}\cdot\left[\mathbf{V}\tilde{p}_{s}({\bf r})-\mathcal{D}\nabla_{{\bf
r}}\tilde{p}_{s}({\bf r})\right]=0\;,$ (C.2)

with:

$\mathbf{V}=-\frac{v_{0}}{d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)}\nabla_{\bf
r}v_{0}\quad\text{and}\quad\mathcal{D}=\frac{v_{0}^{2}}{d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)}+D_{t}\
.$ (C.3)

If we look for a flux-less stationary solution of Eq. (C.2):

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\mathbf{V}\tilde{p}_{s}-\mathcal{D}\nabla_{\bf r}\tilde{p}_{s}$
(C.4) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{v_{0}}{d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)}\nabla_{\bf
r}v_{0}+\left[\frac{v_{0}^{2}}{d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)}+D_{t}\right]\nabla_{\bf
r}\log\tilde{p}_{s}$ (C.5) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\nabla_{\bf
r}\frac{v_{0}^{2}}{2}+\left[v_{0}^{2}+d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)D_{t}\right]\nabla_{\bf
r}\log\tilde{p}_{s}\;.$ (C.6)

Setting $\tilde{p}_{s}({\bf r})=\exp[-\Phi({\bf r})]$ in Eq. (C.6) we obtain:

$\nabla_{\bf r}\Phi({\bf r})=\frac{1}{2}\>\dfrac{\nabla_{\bf r}v_{0}^{2}({\bf
r})}{v_{0}^{2}({\bf r})+d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)D_{t}}\;,$
(C.7)

which is solved by:

$\Phi({\bf r})=\frac{1}{2}\log\left[v_{0}^{2}({\bf
r})+d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)D_{t}\right]\;.$ (C.8)

In conclusion, the coarse-grained stationary distribution reads:

$\tilde{p}_{s}({\bf r})\propto\frac{1}{v_{0}({\bf
r})}\dfrac{1}{\sqrt{1+\dfrac{d\left(\alpha_{0}+(d-1)\Gamma_{0}\right)D_{t}}{v_{0}^{2}({\bf
r})}}}\;.$ (C.9)

Note that for $D_{t}=0$ the coarse-grained solution recovers exactly the
microscopic one (C.1).

###  C.2 AOUPs with position-dependent speed and translational noise

We consider the case of an AOUP with $v=v_{0}({\bf r})$, constant
$\omega=\omega_{0}$ and constant translational noise $D_{t}$.
For finite $D_{t}>0$ the associated master equation (2.65) is not exactly
solvable, so we look for an approximate solution $\tilde{p}_{s}({\bf r})$ of
the coarse-grained dynamics:

$\partial_{t}\tilde{p}_{s}({\bf r})=-\nabla_{{\bf
r}}\cdot\left[\mathbf{V}\tilde{p}_{s}({\bf r})-\mathcal{D}\nabla_{{\bf
r}}\tilde{p}_{s}({\bf r})\right]=0$ (C.10)

with:

$\mathbf{V}=-\frac{v_{0}}{d\omega_{0}}\nabla_{\bf
r}v_{0}\quad\text{and}\quad\mathcal{D}=\frac{v_{0}^{2}}{d\omega_{0}}+D_{t}\ .$
(C.11)

The stationary solution of Eq. (C.10) can be found following the same steps as
in Sec. C.1. Eventually one finds:

$\tilde{p}_{s}({\bf r})\propto\frac{1}{v_{0}({\bf
r})}\dfrac{1}{\sqrt{1+\dfrac{d\omega_{0}D_{t}}{v_{0}^{2}({\bf r})}}}\;.$
(C.12)

##  Appendix D Diffusion-drift approximation for active mixtures

###  D.1 ABP-RTP mixture: Single-particle coarse-graining

In this Appendix we derive the mesoscopic Langevin equation for mixtures of
ABP-RTPs starting from the associated master equation (4.5). The derivation is
carried out under the assumption of space-dependent motility parameters
(4.6)-(4.8), but once again we remark its relevance for the interacting case,
thanks to the frozen-field approximation.
To begin with, we expand $\mathcal{P}_{\mu}({\bf r},{\bf u})$ on the spherical
harmonics basis:

$\mathcal{P}_{\mu}({\bf r},{\bf
u})=\sum_{p=0}^{\infty}\frac{1}{\Omega}\frac{(d-2+2p)!!}{p!(d-2)!!}\>{\bf{a}}_{\mu}^{p}(\mathbb{r})\cdot\widehat{{\bf
u}^{\otimes
p}}\;,\qquad{\bf{a}}_{\mu}^{p}=\int_{\mathbb{S}^{d-1}}\mathrm{d}{\bf
u}\mathcal{P}_{\mu}\widehat{{\bf u}^{\otimes p}}=\langle\widehat{{\bf
u}^{\otimes p}}\;\mathcal{P}_{\mu}\rangle\;.$ (D.1)

Integrating Eq. (4.5) over all possible orientations ${\bf u}$ in the unit
sphere $\mathbb{S}^{d-1}$, we obtain:

$\partial_{t}\langle\mathcal{P}_{\mu}\rangle=-\nabla_{\bf
r}\cdot\left[v_{0\mu}\langle{\bf
u}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle{\bf u}^{\otimes 2}\;\mathcal{P}_{\mu}\rangle-
D_{t\mu}\nabla_{\bf r}\langle\mathcal{P}_{\mu}\rangle\right]\;,$ (D.2)

Using the definition (D.1) of harmonic components ${\bf{a}}_{\mu}^{p}$ and the
decomposition ${\bf u}^{\otimes 2}=\frac{{\bf{I}}}{d}+\widehat{{\bf
u}^{\otimes 2}}$, we obtain:

$\partial_{t}{\bf{a}}^{0}_{\mu}=-\nabla_{\bf
r}\cdot\left[v_{0\mu}{\bf{a}}^{1}_{\mu}-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\left({\bf{a}}^{2}_{\mu}+\frac{{\bf{I}}}{d}{\bf{a}}^{0}_{\mu}\right)-D_{t\mu}\nabla_{\bf
r}{\bf{a}}_{\mu}^{0}\right]\;,$ (D.3)

Next, we multiply our master Eq. (4.9) by $\widehat{{\bf u}^{\otimes 1}}={\bf
u}$ and integrate over ${\bf u}$ to obtain the dynamics of the first harmonic
component:

$\displaystyle\partial_{t}\langle\widehat{{\bf u}^{\otimes
1}}\;\mathcal{P}_{\mu}\rangle=$ $\displaystyle-\nabla_{\bf
r}\cdot\left[v_{0\mu}\langle{\bf u}^{\otimes
2}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle{\bf u}^{\otimes 3}\;\mathcal{P}_{\mu}\rangle-
D_{t\mu}\nabla_{\bf r}\langle\widehat{{\bf u}^{\otimes
1}}\;\mathcal{P}_{\mu}\rangle\right]$ (D.4)
$\displaystyle-\alpha_{0\mu}\langle\widehat{{\bf u}^{\otimes
1}}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}\alpha_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle{\bf u}^{\otimes
2}\;\mathcal{P}_{\mu}\rangle+\Gamma_{0}\langle\widehat{{\bf u}^{\otimes
1}}\Delta_{\bf u}\mathcal{P}_{\mu}\rangle$
$\displaystyle+\left(\sum_{h=1}^{n}\Gamma_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle\widehat{{\bf u}^{\otimes 1}}\otimes\Delta_{\bf
u}({\bf u}\mathcal{P}_{\mu})\rangle$

To simplify this equation, we first write ${\bf u}^{\otimes 3}$ and ${\bf
u}^{\otimes 2}$ as a combination of harmonic tensors, using Eqs.
(2.26)-(2.30). Then, using the fact that $\widehat{{\bf u}^{\otimes p}}$ is an
eigenvector of the Laplacian $\Delta_{\bf u}$ with eigenvalue $-p(p+d-2)$, we
obtain:

$\displaystyle\partial_{t}{\bf{a}}^{1}_{\mu}=$ $\displaystyle-\nabla_{\bf
r}\cdot\left[v_{0\mu}\left({\bf{a}}^{2}_{\mu}+\frac{{\bf{I}}}{d}{\bf{a}}^{0}_{\mu}\right)-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\left({\bf{a}}^{3}_{\mu}+\frac{3}{d+2}{\bf{a}}^{1}_{\mu}\odot{\bf{I}}\right)-D_{t\mu}\nabla_{\bf
r}{\bf{a}}_{\mu}^{1}\right]$ (D.5)
$\displaystyle-\left[\alpha_{0\mu}+(d-1)\Gamma_{0\mu}\right]{\bf{a}}_{\mu}^{1}-\left[\sum_{h=1}^{n}\left(\alpha_{1\mu}^{h}+(d-1)\Gamma_{1\mu}^{h}\right)\nabla_{\bf
r}c_{h}\right]\cdot\left({\bf{a}}^{2}_{\mu}+\frac{{\bf{I}}}{d}{\bf{a}}^{0}_{\mu}\right)\;.$

The dynamics of the second harmonic component ${\bf{a}}_{\mu}^{2}$ can be
obtained in a similar way, by multiplying the master Eq. (4.9) by
$\widehat{{\bf u}^{\otimes 2}}$ and integrating over ${\bf u}$:

$\displaystyle\partial_{t}\langle\widehat{{\bf u}^{\otimes
2}}\;\mathcal{P}_{\mu}\rangle=$ $\displaystyle-\nabla_{\bf
r}\cdot\left[v_{0\mu}\langle{\bf u}\otimes\widehat{{\bf u}^{\otimes
2}}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle{\bf u}^{\otimes 2}\otimes\widehat{{\bf u}^{\otimes
2}}\;\mathcal{P}_{\mu}\rangle-D_{t\mu}\nabla_{\bf r}\langle\widehat{{\bf
u}^{\otimes 2}}\;\mathcal{P}_{\mu}\rangle\right]$ (D.6)
$\displaystyle-\alpha_{0\mu}\langle\widehat{{\bf u}^{\otimes
2}}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}\alpha_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle{\bf u}\otimes\widehat{{\bf u}^{\otimes
2}}\;\mathcal{P}_{\mu}\rangle+\Gamma_{0}\langle\widehat{{\bf u}^{\otimes
2}}\Delta_{\bf u}\mathcal{P}_{\mu}\rangle$
$\displaystyle+\left(\sum_{h=1}^{n}\Gamma_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle\widehat{{\bf u}^{\otimes 2}}\otimes\Delta_{\bf
u}({\bf u}\mathcal{P}_{\mu})\rangle$

Using the expressions (2.26)-(2.30), one gets:

$\displaystyle{\bf u}\otimes\widehat{{\bf u}^{\otimes 2}}$ $\displaystyle=$
$\displaystyle{\bf u}^{\otimes 3}-\frac{1}{d}{\bf
u}\otimes{\bf{I}}\quad\>=\>\widehat{{\bf u}^{\otimes 3}}+\frac{3}{d+2}{\bf
u}\odot{\bf{I}}-\frac{1}{d}{\bf u}\otimes{\bf{I}}\;,$ (D.7)
$\displaystyle\mathbb{u}^{\otimes 2}\otimes\widehat{\mathbb{u}^{\otimes 2}}$
$\displaystyle=$ $\displaystyle\mathbb{u}^{\otimes
4}-\frac{1}{d}\mathbb{u}^{\otimes
2}\otimes{\bf{I}}\>=\>\widehat{\mathbb{u}^{\otimes
4}}+\frac{6}{d+4}\widehat{\mathbb{u}^{\otimes
2}}\odot{\bf{I}}+\frac{3}{d(d+2)}{\bf{I}}^{\odot
2}-\frac{1}{d}\widehat{\mathbb{u}^{\otimes
2}}\otimes{\bf{I}}-\frac{1}{d^{2}}{\bf{I}}^{\otimes 2}\;,$ (D.8)

which can be inserted in Eq. (D.6) to obtain:

$\displaystyle\partial_{t}{\bf{a}}_{\mu}^{2}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}\left({\bf{a}}_{\mu}^{3}+\frac{3}{d+2}{\bf{a}}_{\mu}^{1}\odot{\bf{I}}-\frac{1}{d}{\bf{a}}_{\mu}^{1}\otimes{\bf{I}}\right)-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\Big{(}{\bf{a}}_{\mu}^{4}+\frac{6}{d+4}{\bf{a}}_{\mu}^{2}\odot{\bf{I}}\right.$
(D.9) $\displaystyle\left.\hskip
34.14322pt-\frac{1}{d}{\bf{a}}_{\mu}^{2}\otimes{\bf{I}}+\frac{3\Omega}{d(d+2)}{\bf{I}}^{\odot
2}-\frac{\Omega}{d^{2}}{\bf{I}}^{\otimes
2}\Big{)}-D_{t\mu}\nabla_{\mathbb{r}}\cdot{\bf{a}}_{\mu}^{2}\right]$
$\displaystyle-(\alpha_{0\mu}+2d\Gamma_{0\mu}){\bf{a}}_{\mu}^{2}-\left[\sum_{h=1}^{n}\left(\alpha_{1\mu}^{h}+2d\Gamma_{1\mu}^{h}\right)\nabla_{\bf
r}c_{h}\right]\cdot\left({\bf{a}}_{\mu}^{3}+\frac{3}{d+2}{\bf{a}}_{\mu}^{1}\odot{\bf{I}}-\frac{1}{d}{\bf{a}}_{\mu}^{1}\otimes{\bf{I}}\right)\
.$

In general, projecting the master Eq. (4.5) onto $\widehat{{\bf u}^{\otimes
p}}$ leads to the following dynamics for the harmonic component
${\bf{a}}^{p}$:

$\displaystyle\partial_{t}{\bf{a}}^{p}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}\langle\widehat{{\bf
u}^{\otimes p}}\otimes{\bf
u}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle\widehat{{\bf u}^{\otimes p}}\otimes{\bf u}^{\otimes
2}\;\mathcal{P}_{\mu}\rangle-
D_{t\mu}\nabla_{\mathbb{r}}\cdot{\bf{a}}_{\mu}^{p}\right]$ (D.10)
$\displaystyle-\left[\alpha_{0\mu}+p(p+d-2)\Gamma_{0\mu}\right]{\bf{a}}_{\mu}^{p}-\left[\sum_{h=1}^{n}(\alpha_{1\mu}^{h}+p(p+d-2)\Gamma_{1\mu}^{h})\nabla_{\bf
r}c_{h}\right]\cdot\langle\widehat{{\bf u}^{\otimes p}}\otimes{\bf
u}\mathcal{P}\rangle\;.$

As in the single-species case, the zeroth-order harmonics is a conserved mode
and thus evolves over a slow, diffusive timescale. On the contrary, all modes
$p\geq 1$ undergo a fast exponential relaxation with finite relaxation times
$[\alpha_{0\mu}+p(p+d-2)\Gamma_{0\mu}]^{-1}$. When studying the diffusive
dynamics of ${\bf{a}}^{0}_{\bf u}$ we can therefore assume
$\partial_{t}{\bf{a}}_{\mu}^{p}=0$ for $p\geq 1$, which leads to:

$\displaystyle{\bf{a}}_{\mu}^{p}$ $\displaystyle=$
$\displaystyle\mathcal{O}(\nabla_{\bf r})\;,\hskip 28.45274ptp>2$ (D.11)
$\displaystyle{\bf{a}}_{\mu}^{2}$ $\displaystyle=$
$\displaystyle\mathcal{O}(\nabla_{\bf r}^{2})\;,$ (D.12)
$\displaystyle{\bf{a}}_{\mu}^{1}$ $\displaystyle=$
$\displaystyle-\frac{\nabla_{\mathbb{r}}(\>v_{0\mu}{\bf{a}}_{\mu}^{0}\>)}{d[\alpha_{0\mu}+(d-1)\Gamma_{0\mu}]}-{\bf{a}}_{\mu}^{0}\sum_{h=1}^{n}\frac{\alpha_{1\mu}^{h}+(d-1)\Gamma_{1\mu}^{h}}{d[\alpha_{0\mu}+(d-1)\Gamma_{0\mu}]}\nabla_{\mathbb{r}}c_{h}+\mathcal{O}(\nabla_{\bf
r}^{2})\;.$ (D.13)

To conclude, we insert Eq. (D.11), (D.13) inside the dynamics (D.3) of the
zeroth-order harmonics ${\bf{a}}_{\mu}^{0}$. In the diffusion-drift
approximation, we can truncate the resulting equation including all terms up
to $\mathcal{O}(\nabla_{\bf r}^{2})$. Finally, we obtain the dynamics of
${\bf{a}}^{0}_{\mu}({\bf r},t)$, i.e. the probability of finding a
$\mu$-particle at position ${\bf r}$:

$\partial_{t}{\bf{a}}^{0}_{\mu}=-\nabla_{{\bf
r}}\cdot\left[\mathbf{V}_{\mu}{\bf{a}}_{\mu}^{0}-\mathcal{D}_{\mu}\nabla_{{\bf
r}}{\bf{a}}^{0}_{\mu}\right]\;,$ (D.14)

where the drift and diffusion coefficients read, respectively:

$\begin{split}&\mathbf{V}_{\mu}=-\frac{v_{0\mu}\nabla
v_{0\mu}}{d\left[\alpha_{0\mu}+(d-1)\Gamma_{0\mu}\right]}-\frac{1}{d}\sum_{h=1}^{n}\left[v_{1\mu}^{h}+v_{0\mu}\frac{\alpha_{1\mu}^{h}+(d-1)\Gamma_{1\mu}^{h}}{\alpha_{0\mu}+(d-1)\Gamma_{0\mu}}\right]\nabla_{{\bf
r}}c_{h}\\\
&\mathcal{D}_{\mu}=\frac{v_{0\mu}^{2}}{d\left[\alpha_{0\mu}+(d-1)\Gamma_{0\mu}\right]}+D_{t\mu}\
.\end{split}$ (D.15)

One can then associate to Eq. (D.14) the following Itô-Langevin dynamics for
particle $i$ of species $\mu$:

$\dot{{\bf r}}_{i,\mu}=\mathbf{V}_{\mu}({\bf
r}_{i,\mu},[\\{\rho_{\nu}\\}])+\nabla_{{\bf r}_{i,\mu}}\mathcal{D}_{\mu}({\bf
r}_{i,\mu},[\\{\rho_{\nu}\\}])+\sqrt{2\mathcal{D}_{\mu}({\bf
r}_{i,\mu},[\\{\rho_{\nu}\\}])}\boldsymbol{\xi}_{i,\mu}(t)\;,$ (D.16)

where the $\\{\boldsymbol{\xi}_{i,\mu}(t)\\}$ are delta-correlated, centred
Gaussian white noise terms.

###  D.2 AOUP mixture: Single-particle coarse-graining

In this Appendix we carry out the coarse-graining for AOUPs with spatial
dependent motility parameters according to (4.10), (4.11). We start from the
master equation for $\mathcal{P}_{\mu}({\bf r},\mathbb{e})$:

$\displaystyle\partial_{t}\mathcal{P}_{\mu}(\mathbb{r},\mathbb{e})$
$\displaystyle=$ $\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}({\bf
r})\mathbb{e}\mathcal{P}_{\mu}-\mathbb{e}^{\otimes
2}\cdot\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf r}c_{h}({\bf
r})\right)\mathcal{P}_{\mu}-D_{t\mu}\nabla_{\mathbb{r}}\mathcal{P}_{\mu}\right]$
(D.17) $\displaystyle+\nabla_{\mathbb{e}}\cdot\left[\omega_{0\mu}({\bf
r})\mathbb{e}\mathcal{P}_{\mu}+\mathbb{e}^{\otimes
2}\cdot\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf r}c_{h}({\bf
r})\right)\mathcal{P}_{\mu}\right]$ $\displaystyle+\frac{1}{d}\Delta_{\bf
e}\left[\omega_{0\mu}({\bf
r})\mathcal{P}_{\mu}+\mathbb{e}\cdot\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}({\bf r})\right)\mathcal{P}_{\mu}\right]\;.$

As we did in Sec. 2.3, we define the $p$-th order moment of
$\mathcal{P}_{\mu}$ as:

${\bf m}^{p}_{\mu}({\bf
r})\equiv\int_{\mathbb{R}^{d}}\mathrm{d}\mathbb{e}\>\mathbb{e}^{\otimes
p}\mathcal{P}_{\mu}({\bf r},\mathbb{e})=\langle{\bf e}^{\otimes
p}\;\mathcal{P}_{\mu}\rangle\;.$ (D.18)

In particular, we note that ${\bf m}_{\mu}^{0}({\bf r},t)$ corresponds to the
probability of finding one particle of type $\mu$ at position ${\bf r}$ at
time $t$, marginalized over all possible orientations ${\bf e}$.
Multiplying the master Eq. (D.17) by $\mathbb{e}^{\otimes p}$ and integrating
over ${\bf e}$ we obtain:

$\displaystyle\partial_{t}\langle{\bf e}^{\otimes
p}\;\mathcal{P}_{\mu}\rangle$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}\langle{\bf e}^{\otimes
p+1}\;\mathcal{P}_{\mu}\rangle-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot\langle{\bf e}^{\otimes p+2}\;\mathcal{P}_{\mu}\rangle-
D_{t\mu}\nabla_{\mathbb{r}}\langle{\bf e}^{\otimes
p}\;\mathcal{P}_{\mu}\rangle\right]$ (D.19) $\displaystyle+\Bigg{\langle}{\bf
e}^{\otimes
p}\nabla_{\mathbb{e}}\cdot\left[\omega_{0\mu}\mathbb{e}\mathcal{P}_{\mu}+\mathbb{e}^{\otimes
2}\cdot\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right)\mathcal{P}_{\mu}\right]\Bigg{\rangle}$
$\displaystyle+\frac{1}{d}\Bigg{\langle}\>{\bf e}^{\otimes p}\Delta_{\bf
e}\left[\omega_{0\mu}\mathcal{P}_{\mu}+\mathbb{e}\cdot\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right)\mathcal{P}_{\mu}\right]\Bigg{\rangle}\;.$

Using the results obtained in Sec. 2.3:

$\displaystyle\omega_{0\mu}\Bigg{\langle}{\bf e}^{\otimes
p}\nabla_{\mathbb{e}}\cdot(\mathbb{e}\mathcal{P}_{\mu})\Bigg{\rangle}$
$\displaystyle=$ $\displaystyle-p\,\omega_{0\mu}{\bf m}^{p}_{\mu}$ (D.20)
$\displaystyle\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right)\cdot\Bigg{\langle}{\bf e}^{\otimes
p}\nabla_{\mathbb{e}}\cdot(\mathbb{e}^{\otimes
2}\mathcal{P}_{\mu})\Bigg{\rangle}$ $\displaystyle=$
$\displaystyle-p\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right){\bf m}^{p+1}_{\mu}$ (D.21)
$\displaystyle\omega_{0\mu}\Bigg{\langle}{\bf e}^{\otimes p}\Delta_{\bf
e}\mathcal{P}_{\mu}\Bigg{\rangle}$ $\displaystyle=$ $\displaystyle
p(p-1)\omega_{0\mu}{\bf m}_{\mu}^{p-2}\odot{\bf{I}}$ (D.22)
$\displaystyle\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right)\cdot\Bigg{\langle}{\bf e}^{\otimes p}\Delta_{\bf
e}(\mathbb{e}\mathcal{P}_{\mu})\Bigg{\rangle}$ $\displaystyle=$ $\displaystyle
p(p-1)\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf r}c_{h}\right){\bf
m}^{p-1}_{\mu}\odot{\bf{I}}$ (D.23)

we can eventually rewrite Eq. (D.19) as:

$\displaystyle\partial_{t}{\bf m}^{p}_{\mu}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}{\bf
m}^{p+1}_{\mu}-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot{\bf m}^{p+2}_{\mu}-D_{t\mu}\nabla_{\mathbb{r}}{\bf
m}^{p}_{\mu}\right]-p\,\omega_{0\mu}{\bf
m}^{p}_{\mu}-p\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right){\bf m}^{p+1}_{\mu}$ (D.24)
$\displaystyle+p(p-1)\frac{\omega_{0\mu}}{d}{\bf
m}_{\mu}^{p-2}\odot{\bf{I}}+p(p-1)\frac{1}{d}\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right){\bf m}^{p-1}_{\mu}\odot{\bf{I}}$

From Eq. (D.24) we find that the zeroth-order mode is a slow conserved field,
whose dynamics is given by:

$\displaystyle\partial_{t}{\bf m}^{0}_{\mu}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}{\bf
m}^{1}_{\mu}-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot{\bf m}^{2}_{\mu}-D_{t\mu}\nabla_{\mathbb{r}}{\bf
m}^{0}_{\mu}\right]$ (D.25)

Conversely, all modes $\mathbb{m}_{\mu}^{p}$ with $p\geq 1$ relax
exponentially fast with a characteristic time $(p\,\omega_{0\mu})^{-1}$. We
can therefore use a fast–variable approximation for all
$\mathbb{m}_{\mu}^{p}$, $p\geq 1$, setting and set
$\partial_{t}\mathbb{m}_{\mu}^{p}$ for all $p\geq 1$ to zero in Eq. (D.24):

$\displaystyle p\,\omega_{0\mu}\mathbb{m}^{p}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbb{r}}\cdot\left[v_{0\mu}{\bf
m}^{p+1}_{\mu}-\left(\sum_{h=1}^{n}v_{1\mu}^{h}\nabla_{\bf
r}c_{h}\right)\cdot{\bf m}^{p+2}_{\mu}-D_{t\mu}\nabla_{\mathbb{r}}{\bf
m}^{p}_{\mu}\right]-p\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right){\bf m}^{p+1}_{\mu}$ (D.26)
$\displaystyle+p(p-1)\frac{\omega_{0\mu}}{d}{\bf
m}_{\mu}^{p-2}\odot{\bf{I}}+p(p-1)\frac{1}{d}\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right){\bf m}^{p-1}_{\mu}\odot{\bf{I}}$

Eq. (D.26) provides a bound on the scaling of the moments
$\mathbb{m}_{\mu}^{p}$ in gradients:

$\forall l\in\mathbb{N}\
,\quad\mathbb{m}_{\mu}^{2l}=\mathcal{O}(1)\quad\text{while}\quad\mathbb{m}_{\mu}^{2l+1}=\mathcal{O}(\nabla_{\mathbb{r}})\
,$ (D.27)

as well as the more precise scalings of $\mathbb{m}_{\mu}^{2}$:

$\mathbb{m}_{\mu}^{2}=\frac{1}{d}\mathbb{m}_{\mu}^{0}{\bf{I}}+\mathcal{O}(\nabla_{\mathbb{r}})\
,$ (D.28)

and $\mathbb{m}_{\mu}^{1}$:

$\mathbb{m}_{\mu}^{1}=-\frac{1}{d\omega_{0\mu}}\nabla_{\mathbb{r}}(v_{0\mu}\mathbb{m}_{\mu}^{0})-\left(\sum_{h=1}^{n}\omega_{1\mu}^{h}\>\nabla_{\bf
r}c_{h}\right)\mathbb{m}_{\mu}^{0}+\mathcal{O}(\nabla_{\mathbb{r}}^{2})\ .$
(D.29)

Finally, we can insert Eqs. (D.28)– (D.29) into Eq. (D.25) to close the
dynamics of $\mathbb{m}_{\mu}^{0}$. We truncate the resulting equation at the
second order in gradient, thus obtaining:

$\partial_{t}\mathbb{m}_{\mu}^{0}=-\nabla_{{\bf
r}}\cdot\left[\mathbf{V}_{\mu}\mathbb{m}_{\mu}^{0}-\mathcal{D}_{\mu}\nabla_{{\bf
r}}\mathbb{m}_{\mu}^{0}\right]\;,$ (D.30)

where the mesoscopic drifts and diffusivities are respectively given by:

$\mathbf{V}_{\mu}=-\frac{v_{0\mu}\nabla
v_{0\mu}}{d\omega_{0\mu}}-\frac{1}{d}\sum_{h=1}^{n}\left[v_{1\mu}^{h}+v_{0\mu}\frac{\omega_{1\mu}^{h}}{\omega_{0\mu}}\right]\nabla_{{\bf
r}}c_{h}\qquad\text{and}\qquad\mathcal{D}_{\mu}=\frac{v_{0\mu}^{2}}{d\omega_{0\mu}}+D_{t\mu}\
.$ (D.31)

One can then associate to Eq. (D.30) the following Itô-Langevin dynamics for
particle $i$ of species $\mu$:

$\dot{{\bf r}}_{i,\mu}=\mathbf{V}_{\mu}({\bf
r}_{i,\mu},[\\{\rho_{\nu}\\}])+\nabla_{{\bf r}_{i,\mu}}\mathcal{D}_{\mu}({\bf
r}_{i,\mu},[\\{\rho_{\nu}\\}])+\sqrt{2\mathcal{D}_{\mu}({\bf
r}_{i,\mu},[\\{\rho_{\nu}\\}])}\boldsymbol{\xi}_{i,\mu}(t)\;,$ (D.32)

where the $\\{\boldsymbol{\xi}_{i,\mu}(t)\\}$ are delta-correlated, centred
Gaussian white noise terms.