# Nonreciprocal collective dynamics in a mixture of phoretic Janus colloids

Gennaro Tucci  Max Planck Institute for Dynamics and Self-Organization
(MPIDS), D-37077 Göttingen, Germany  Ramin Golestanian
<EMAIL_ADDRESS>Max Planck Institute for Dynamics and Self-
Organization (MPIDS), D-37077 Göttingen, Germany  Rudolf Peierls Centre for
Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom
Suropriya Saha<EMAIL_ADDRESS>Max Planck Institute for Dynamics and
Self-Organization (MPIDS), D-37077 Göttingen, Germany

###### Abstract

A multicomponent mixture of Janus colloids with distinct catalytic coats and
phoretic mobilities is a promising theoretical system to explore the
collective behavior arising from nonreciprocal interactions. An active colloid
produces (or consumes) chemicals, self-propels, drifts along chemical
gradients, and rotates its intrinsic polarity to align with a gradient. As a
result the connection from microscopics to continuum theories through coarse-
graining couples densities and polarization fields in unique ways. Focusing on
a binary mixture, we show that these couplings render the unpatterned
reference state unstable to small perturbations through a variety of
instabilities including oscillatory ones which arise on crossing an
exceptional point or through a Hopf bifurcation. For fast relaxation of the
polar fields, they can be eliminated in favor of the density fields to obtain
a microscopic realization of the Nonreciprocal Cahn-Hilliard model for two
conserved species with two distinct sources of non-reciprocity, one in the
interaction coefficient and the other in the interfacial tension. Our work
establishes Janus colloids as a versatile model for a bottom-up approach to
both scalar and polar active mixtures.

##  I Introduction

As a route to manifesting active matter systems [1], the breaking of action-
reaction symmetry in effective interactions, or non-reciprocity, has recently
garnered increased attention [2]. Reciprocity in interactions is synonymous
with the existence of an interaction potential. If the concept of an effective
free energy cannot be applied, as is very likely to be the case for
interactions driven by chemicals [3, 4], social interactions [5], velocity
fields [6, 7, 8], or information transfer [9, 10, 11] non-reciprocity will
inevitably emerge, whether or not it is significant at long timescales [12].

In a system of particles without an intrinsic polarity, non-reciprocity is
apparent only in active mixtures which allow breaking Newton’s third law in
pairwise interactions. For example, activity manifests in the formation of
novel bound states in collections of uniformly coated active colloids [13, 14,
15, 16]. Striking collective behavior emerges in large collections of
chemically active colloids [17, 18, 19, 20] involving chasing dynamics.
Continuum theories that minimally capture the essence of nonreciprocal
interactions in scalar mixtures have been proposed [21, 22] and are being
explored intensely [23, 24, 25]. In polar active matter, with orientation as a
relevant degree of freedom, non-reciprocity can be incorporated in a multitude
of ways - directly in the alignment rules [26, 27], through a dependence of
the spin-spin interactions on spatial anisotropy [28, 29, 30], or through
quorum sensing [31]. Explorations of the collective behavior of systems with
polar, nematic, or chiral order constitute an active sub-field of research
[26, 10, 32]. Novel steady states arise in all the examples mentioned here due
to the simultaneous breaking of parity and time reversal symmetry leading to
chiral motion in polar mixtures [26, 6], traveling waves in scalar mixtures
[21], and stress-strain cycles in odd solids [33, 34] relating this class of
phenomenon to odd response [35].

A few examples of experimental systems that exhibit nonreciprocal interactions
are - active Janus colloids [36, 37, 38], light actuated colloids with a
vision cone [39, 40], and dusty plasma [41]. In this work, we will focus on
Janus colloids, where nonreciprocal interactions between densities and
polarities are realized through chemical field-mediated interactions. The
nonequilibrium active dynamics of self-propelled Janus colloids are due to
self-phoresis [36], which harnesses the force-free mechanism of
diffusiophoresis at microscopic lengthscales [42, 4]. These particles are able
to catalyze a chemical reaction on their surface which modifies the density
profile of the involved reactants and products [43]. Depending on the
geometrical properties of the colloid, changes in these chemical substrates’
concentration may lead to the particle net motion [44]. Moreover, the coupling
of different colloids to the same substrates induces an effective long-ranged
interaction among different particles.

The versatility of a system of Janus colloids arises from the variety in their
dynamical response [45, 46]. In a collection of identical colloids, possible
collective dynamics include system-wide phase separation, pattern formation
with a selected lengthscale, and oscillations similar to the Jeans instability
in gravitational systems as explored in [45]. Numerical solutions of the
dynamics for chemorepulsive colloids in [47, 48] show that the linear
instabilities indeed pave the route to dynamic aggregates, and spatiotemporal
patterns including traveling waves. For two interacting Janus colloids, the
interplay of orientational dynamics, self-propulsion, and drift produces a
complex effective potential landscape leading to bound orbits with internal
jiggling and chiral bound state [49].

In this paper, we have studied the collective behavior of two species of
phoretic colloids starting from their microscopic dynamics and building the
continuum field theories. In section II, we introduce the model of multi-
species Janus colloids. In section III, we discuss the linear stability
analysis of the two-species case, while in appendix C the one-species one.
Finally in section IV, we study what happens to the two-species system
whenever the polarization fields relax fast enough so that they can be
adiabatically eliminated.

Figure 1: Schematic showcasing the complexity in the dynamics that can be
harnessed to tune the system to the collective behavior of choice. (a) shows a
Janus colloid that produces a chemical field (profile shown by red heatmap)
with rate $\alpha>0$. $\beta$ determines the asymmetry in production along the
polarity $\bm{n}$. (b) is similar to (a) but with the sign of $\beta$ reversed
thus flipping the profile of the chemical field in (a). (c) lists all possible
types of center-to-center interaction between two different species of Janus
colloids (blue and green) due to chemotactic drift. To fix ideas, we
illustrate pairwise interactions only, and the two involved particles are
joined by a dotted black line. For both $\mu_{a}$ positive (negative), all
interactions are repulsive (attractive). For
$\operatorname{sign}(\mu_{1})\neq\operatorname{sign}({\mu}_{2})$,
nonreciprocal interactions emerge where one species chases the other. (d)
shows the effect of chemotactic alignment: whenever $\Omega_{a}$ is positive
(negative), particles re-orient towards (away from) high-density regions
resulting in a novel form of orientational order. Orientational dynamics leads
to effective pairwise attraction-repulsion or chasing depending on the sign of
$\Omega_{a}$.

##  II Multi-species Janus particles

We study the dynamics of $n$ different species of Janus colloids [50, 45]
interacting with a chemical substrate. These particles are sensitive to the
spatial gradient of the substrate and respond accordingly in two different
ways: (i) their velocity varies proportionally to the gradient of the chemical
concentration (chemotactic drift), (ii) they re-orient along it (chemotactic
alignment). Moreover, each Janus particle contributes to the production (or
consumption) of the substrate particles. This inhomogeneity in the chemical in
the proximity of the Janus particle’s surface produces a slip velocity, which
induces self-propulsion in the direction of the particle’s axis [44]. An
effective nonreciprocal interaction between the two species follows from the
direct interplay with the chemical. We describe the system by looking at the
3-dimensional dynamics of Janus colloids of species $a\in\\{1,2,\dots,n\\}$
with position $\bm{r}_{a}$ and orientation $\bm{n}_{a}$ which specifies the
direction of self-propulsion. They evolve in time according to the Langevin
equations

$\displaystyle\frac{\mathrm{d}\bm{r}_{a}}{\mathrm{d}t}$
$\displaystyle=-\mu_{a}\nabla c+v_{a}\bm{n}_{a}+\bm{\xi}_{a},$ (1)
$\displaystyle\frac{\mathrm{d}\bm{n}_{a}}{\mathrm{d}t}$
$\displaystyle=\bm{\omega}_{a}\times\bm{n}_{a}$
$\displaystyle=(\Omega_{a}\bm{n}_{a}\times\nabla
c+\bm{\zeta}_{a})\times\bm{n}_{a}$
$\displaystyle=\Omega_{a}\left(\mathbb{I}-\bm{n}_{a}\bm{n}_{a}\right)\cdot\nabla
c+\bm{\zeta}_{a}\times\bm{n}_{a},$

here $\mathbb{I}$ denotes the 3-dimensional identity matrix, and $v_{a}$ the
self-propulsion speed. The sign of the chemotactic mobility $\mu_{a}$
prescribes whether the particle moves following increasing ($\mu_{a}<0$) or
decreasing ($\mu_{a}>0$) gradients of the substrate density field
$c(\bm{r},t)$. Similarly, the coefficient $\Omega_{a}$ leads to alignment
($\Omega_{a}>0$) or anti-alignment ($\Omega_{a}<0$) of the particle velocity
along the gradient $\nabla c(\bm{r},t)$. In writing equation (1) we ignore
terms that are quadratic in $\bm{n}_{a}$, keeping only those that are leading
order in $\bm{n}_{a}$. This simplification is consistent with ignoring nematic
order during coarse-graining as we will discuss later; we refer to appendix A
for more details. The fluctuating nature of equation (1) is encoded in the
Gaussian white noises $\bm{\xi}_{a}(t)$ and $\bm{\zeta}_{a}(t)$, characterized
by zero mean and variance

$\langle\zeta_{a,i}(t)\zeta_{b,j}(t)\rangle=2D_{r,a}\delta_{ab}\,\delta_{ij}\,\delta(t-t^{\prime}),\;\;\text{and}\;\;\langle\xi_{a,i}(t)\xi_{b,j}(t)\rangle=2D_{a}\delta_{ab}\,\delta_{ij}\,\delta(t-t^{\prime}),$
(2)

where the labels $i,\,j$ identify the 3-dimensional space components of the
vectors. Note that the conservation of the modulus $|\bm{n}_{a}|^{2}=1$ from
equation (1) implicitly assumes the Stratonovich representation of the
stochastic differential equation [51, 52]. The set of dynamical equations (1)
is completed by including the evolution of the substrate density field $c(t)$,
which is given by

$\partial_{t}c-D_{c}(\nabla^{2}-\kappa^{2})c=\sum\left[\alpha_{a}\rho_{a}-\beta_{a}\nabla\cdot\bm{p}_{a}\right].$
(3)

The fields $\rho_{a}(\bm{r},t)$ and $\bm{p}_{a}(\bm{r},t)$ denote respectively
the particle density and the polarization field associated with the species
$a$, and they are defined as follows

$\rho_{a}=\Big{\langle}\sum_{i}\delta(\bm{r}-\bm{r}_{a,i})\Big{\rangle},\;\;\;\;\bm{p}_{a}=\Big{\langle}\sum_{i}\bm{n}_{a,i}\delta(\bm{r}-\bm{r}_{a,i})\Big{\rangle},$
(4)

where $i$ runs over all particles of the $a-$th species, and the average is
evaluated with respect to the noise realizations. As clear from equation (3),
the substrate diffuses with diffusion constant $D_{c}$ and degrades
exponentially over the timescale $1/D_{c}\kappa^{2}$, associated with the
screening lengthscale $\kappa^{-1}$. We also assume that Janus particles
constitute point-like sinks or sources for the substrate. If the sign of the
parameter $\alpha_{a}$ in equation (3) is positive (negative) the particles of
species $a$ produce (consume) the substrate. Similarly, the dipole term
$\beta_{a}$ accounts for the head-tail asymmetry in the production (or
consumption) along $\bm{n}_{a}$. If $\beta_{a}>0$ the Janus particles produce
more (or consume less) substrate around the catalytic cap of $\bm{n}_{a}$,
while they produce less (or consume more) otherwise. We represent this in
figure 1 (a)-(b), where, the net production of the chemical ($\alpha>0$) is
modulated along its symmetry axis by the sign of $\beta$. In figure 1 (c), we
consider the joint effect of the production of the chemical and chemotaxis for
two different species of Janus colloids, marked by the two different blue and
green colors. In particular, we assume the scenario represented in figure 1
(a), where the chemical gradient grows in the self-propelling direction. If we
consider the interactions to be pairwise, given a Janus particle in the
picture, the bold blue and green arrows identify the net qualitative force due
to the colloid of the same color along the dotted line. In the case where
$\mu_{1,2}>0$, the particles, being a source of chemical, are effectively
repelled by each other. If $\mu_{2}<0$ changes sign, green colloids are
attracted by other particles, while blue continue to be repelled: this is a
prototypical manifestation of effective nonreciprocal interactions. In the
last case, being $\mu_{1,2}<0$, the particles are all effectively attracted by
each other. Similarly, in figure 1 (d), we represent the effect of pairwise
chemotactic alignment for the two species of particles. Whenever
$\Omega_{1,2}>0$ the particles tend to align following higher gradients of the
chemical concentration $\nabla c$, which in this example coincides with the
self-propelling direction (tiny arrow), while they anti-align along it
otherwise.

Our goal is to characterize the various dynamical steady states by looking at
the effect of chemically mediated interactions on the particle distribution
$\rho_{a}$ and orientation $\bm{p}_{a}$. It follows from equation (4) that
$\bm{p}_{a}$ and $\rho_{a}$ are coupled fields, for instance, $\bm{p}_{a}$
vanishes at a point where $\rho_{a}$ is zero. As we discuss in A, $\rho_{a}$
and $\bm{p}_{a}$ are respectively the zeroth and first moment in $\bm{n}_{a}$
of the joint probability density
$\mathcal{P}_{a}=\langle\sum_{i}\delta(\bm{r}-\bm{n}_{a,i})\delta(\bm{r}-\bm{r}_{a,i})\rangle$.
In general, one can construct an infinite hierarchy of equations where the
time evolution of the $m-$th moment of $\bm{n}$ depends on the $m^{\prime}-$th
moment, with $m^{\prime}\geq m$. In the spirit of describing the collective
behavior and keeping only the most relevant fields, i.e., those reflecting
conservation laws or broken symmetry in the system, we restrict our analysis
to the first and second moments of the orientation $\bm{n}$ by truncating the
corresponding hierarchy of infinitely many equations in the $\bm{n}$ moment
expansion [3]; the details of the derivation are reported in appendix A. Thus,
we get equations for the coarse-grained fields $\rho_{a}$ and $\bm{p}_{a}$

$\partial_{t}\rho_{a}=-\nabla\cdot\left(v_{a}\bm{p}_{a}-\mu_{a}\rho_{a}\nabla
c\right)+D_{a}\nabla^{2}\rho_{a},$ (5)
$\partial_{t}\bm{p}_{a}=(-2D_{r,a}+D_{a}\nabla^{2})\bm{p}_{a}+\mu_{a}\nabla\cdot\left(\bm{p}_{a}\nabla
c\right)+\frac{2}{3}\Omega_{a}\rho_{a}\nabla c-\frac{v_{a}}{3}\nabla\rho_{a},$
(6)

that allow us, together with equation (3), to describe the macroscopic
behavior of the system under the assumption of negligible nematic order
parameter
$\mathbb{Q}_{a}=\langle\sum_{i}(\bm{n}_{a,i}\bm{n}_{a,i}-\mathbb{I}/3)\delta(\bm{r}-\bm{r}_{a,i})\rangle$.
Equation (5) describes the conserved evolution of the particle density
$\rho_{a}$ via three different contributions: the first accounts for advection
of particles because of self-propulsion, the second is the effect of
chemotaxis, and the third is translational diffusion. Similarly, in equation
(6) for the polar field $\bm{p}_{a}$, the first term represents the effect of
orientational and translational diffusion. By its very definition in equation
(4), $\bm{p}_{a}$ is coupled with $\rho_{a}$ – consequently, both diffusion
and phoretic drift of $\rho_{a}$ affect $\bm{p}_{a}$ through the second and
the third terms respectively. To the lowest order in spatial gradients,
$\bm{p}_{a}$ rotates to align with the local substrate gradient, an effect
that is encoded in the third term proportional to the coefficient of alignment
$\Omega_{a}$. Finally, the last term can be interpreted as a pressure term
that measures how self-propulsion influences local order.

###  II.1 Linearized dynamics

In general, equations (5), (6), complemented with equation (3) constitute a
set of $2n+1$ nonlinear partial differential coupled equations for an equal
number of scalar and vector fields. A great reduction of their complexity is
achieved by studying their linearized form, by looking at perturbations with
respect to a simple and physically relevant solution. We perform linear
stability analysis of equations (3), (5), and (6) by considering small
perturbations around the spatially homogeneous solution

$\rho_{a}=\bar{\rho}_{a},\quad\bm{p}_{a}=0,\quad
c=\bar{c}=\frac{\sum_{a}\alpha_{a}\bar{\rho}_{a}}{D_{c}\kappa^{2}},$ (7)

corresponding to a state where the $n$ species are well mixed, i.e., there are
modulations in the density, and no orientational order is present. The
stationary substrate density $\bar{c}>0$ is constant, as its net production by
Janus colloids is balanced by its degradation rate. We are interested in
solutions that are perturbations of equation (7) in the form of

$\displaystyle\rho_{a}(\bm{r},t)=\bar{\rho}_{a}+\delta\rho_{a}(\bm{r},t),\quad\bm{p}_{a}(\bm{r},t)=\delta\bm{p}_{a}(\bm{r},t),\quad
c(\bm{r},t)=\bar{c}+\delta c(\bm{r},t).$ (8)

We can now build the dynamical equations for $\delta\rho_{a}$,
$\delta\bm{p}_{a}$ and $\delta c$ by discarding contributions that are higher
than the linear order in the perturbations. This procedure leads to the
following set of linearized equations

$\displaystyle\partial_{t}\delta\rho_{a}(\bm{r},t)=\mu_{a}\bar{\rho}_{a}\nabla^{2}\delta
c(\bm{r},t)+D_{a}\nabla^{2}\delta\rho_{a}(\bm{r},t)-v_{a}\nabla\cdot\delta\bm{p}_{a}(\bm{r},t),$
(9)
$\displaystyle\partial_{t}\delta\bm{p}_{a}(\bm{r},t)=-\frac{v_{a}}{3}\nabla\delta\rho_{a}(\bm{r},t)+\frac{2}{3}\Omega_{a}\bar{\rho}_{a}\nabla\delta
c(\bm{r},t)+(D_{a}\nabla^{2}-2D_{r,a})\delta\bm{p}_{a}(\bm{r},t),$
$\displaystyle\partial_{t}\delta
c(\bm{r},t)=-D_{c}\left(\kappa^{2}-\nabla^{2}\right)\delta
c(\bm{r},t)+\sum_{a}\left[\alpha_{a}\delta\rho_{a}(\bm{r},t)-\beta_{a}\nabla\cdot\delta\bm{p}_{a}(\bm{r},t)\right].$

The system of equations (9) represents the starting point to discuss the
linear stability of the disordered state and the onset of order in the system.
Although linear, the current form of equations (9) is still very complex due
to the large number of fields involved and the related parameters. A further
simplification follows from the physical assumption that the deviations of the
substrate density $\delta c$ from the space homogeneous solution relax much
faster than those of the density and polarization field of the Janus
particles, that is $\partial_{t}\delta c(\bm{r},t)\simeq 0$, or equivalently

$D_{c}\left(\kappa^{2}-\nabla^{2}\right)\delta
c(\bm{r},t)=\sum_{a}\left[\alpha_{a}\delta\rho_{a}(\bm{r},t)-\beta_{a}\nabla\cdot\delta\bm{p}_{a}(\bm{r},t)\right].$
(10)

Furthermore, it is convenient to express the time evolution of
$\delta\rho_{a}(\bm{k},t)$, $\delta\bm{p}_{a}(\bm{k},t)$, and $\delta
c(\bm{k},t)$ in terms of their Fourier modes, whose dynamics directly follow
from equations (9) and (10) as

$\displaystyle\partial_{t}\delta\rho_{a}(\bm{k},t)=-k^{2}\left[\mu_{a}\bar{\rho}_{a}\delta
c(\bm{k},t)+D_{a}\delta\rho_{a}(\bm{k},t)\right]-v_{a}(i\bm{k})\cdot\delta\bm{p}_{a}(\bm{k},t),$
(11)
$\displaystyle\partial_{t}\delta\bm{p}_{a}(\bm{k},t)=-\frac{i\bm{k}}{3}\left[v_{a}\delta\rho_{a}(\bm{k},t)-2\Omega_{a}\bar{\rho}_{a}\delta
c(\bm{k},t)\right]-(D_{a}k^{2}+2D_{r,a})\delta\bm{p}_{a}(\bm{k},t),$
$\displaystyle D_{c}\left(\kappa^{2}+k^{2}\right)\delta
c(\bm{k},t)=\sum_{a}\left[\alpha_{a}\delta\rho_{a}(\bm{k},t)-\beta_{a}(i\bm{k})\cdot\delta\bm{p}_{a}(\bm{k},t)\right];$

as a convention the Fourier transform of a given function $f(\bm{r})$ reads
$f(\bm{k})\equiv\int_{-\infty}^{+\infty}\mathrm{d}\bm{r}\,e^{-i\bm{k}\cdot\bm{r}}\,f(\bm{r})$.
For the sake of completeness, we mention that the closure of the set
linearized equations at the nematic field $\mathbb{Q}_{a}$ order would simply
lead to renormalization of $D_{r,a}$ to $D_{r,a}+k^{2}v_{a}^{2}/(45D_{r,a})$,
thus leaving the qualitative behavior of the system unaltered.

We can reduce the degrees of freedom of the problem by looking at the
transverse and longitudinal components of $\delta\bm{p}_{a}(\bm{k},t)$. The
initial value of the perturbation $\delta\bm{p}_{a}(\bm{k},t)$ in the
direction transverse to the wave vector $\bm{k}$ decays with relaxation time
$D_{r,a}+D_{a}k^{2}$, which remains positive and finite at all wave numbers.
Moreover, the dynamics of $\delta\rho_{a}(\bm{k},t)$ and $\delta c(\bm{k},t)$
in equations (11) depends on those of $\delta\bm{p}_{a}(\bm{k},t)$ only via
$i\bm{k}\cdot\delta\bm{p}_{a}(\bm{k},t)$, i.e., the divergence of the polar
field in real-space coordinates. Accordingly, we decompose the polar field
along $\hat{\bm{k}}=\bm{k}/k$ and in the transverse direction as

$\displaystyle\delta\bm{p}_{a}(\bm{k},t)$
$\displaystyle=(\hat{\bm{k}}\cdot\delta\bm{p}_{a}(\bm{k},t))\,\hat{\bm{k}}+(\mathbb{I}-\hat{\bm{k}}\hat{\bm{k}})\cdot\delta\bm{p}_{a}(\bm{k},t)$
(12) $\displaystyle=\delta
p_{\parallel,a}(\bm{k},t)\,\hat{\bm{k}}+\delta\bm{p}_{\perp,a}(\bm{k},t),$

where $\delta p_{\parallel,a}(\bm{k},t)$ is the parallel (longitudinal)
component of $\delta\bm{p}_{a}(\bm{k},t)$ to $\hat{\bm{k}}$ and
$\delta\bm{p}_{\perp,a}(\bm{k},t)$ is a vector which belongs to the plane
perpendicular to $\hat{\bm{k}}$ (transverse component). As anticipated, one
can check that the perpendicular component $\delta\bm{p}_{\perp,a}(\bm{k},t)$
of the polar field is decoupled from the other fields and it relaxes
exponentially in time according to
$\partial_{t}\delta\bm{p}_{\perp,a}(\bm{k},t)=-(D_{a}k^{2}+2D_{r,a})\delta\bm{p}_{\perp,a}(\bm{k},t)$.
Finally, by substituting the explicit expression of $\delta c(\bm{k},t)$ in
equations (11) for $\delta\rho_{a}(\bm{k},t)$ and
$\delta\bm{p}_{a}(\bm{k},t)$, the dynamics of the system are encoded in those
of the two slow (scalar) fields

$\displaystyle\partial_{t}\delta\rho_{a}(\bm{k},t)=\sum_{b}\left[k^{2}\mathcal{G}_{ab}^{tt}(k)\,\delta\rho_{b}(\bm{k},t)+ik\,\mathcal{G}^{tr}_{ab}(k)\,\delta
p_{\parallel,b}(\bm{k},t)\right],$ (13) $\displaystyle\partial_{t}\delta
p_{\parallel,a}(\bm{k},t)=\sum_{b}\left[ik\,\mathcal{G}^{rt}_{ab}(k)\,\delta\rho_{b}(\bm{k},t)+\mathcal{G}_{ab}^{rr}(k)\,\delta
p_{\parallel,b}(\bm{k},t)\right],$

where we have defined the $k-$dependent matrices of the inter-species
couplings

$\begin{split}&\mathcal{G}_{ab}^{tt}(k)\equiv-
D_{a}\delta_{ab}-\mu_{a}\bar{\rho}_{a}\alpha_{b}/[D_{c}(\kappa^{2}+k^{2})],\\\
&\mathcal{G}^{rt}_{ab}(k)\equiv-\delta_{ab}v_{a}/3+2\Omega_{a}\bar{\rho}_{a}\alpha_{b}/[3D_{c}(\kappa^{2}+k^{2})],\end{split}\quad\quad\begin{split}&\mathcal{G}^{tr}_{ab}(k)\equiv-
v_{a}\delta_{ab}+\mu_{a}\bar{\rho}_{a}\beta_{b}\,k^{2}/[D_{c}(\kappa^{2}+k^{2})],\\\
&\mathcal{G}_{ab}^{rr}(k)\equiv-(D_{a}k^{2}+2D_{r,a})\delta_{ab}+2\Omega_{a}\bar{\rho}_{a}\beta_{b}\,k^{2}/[3D_{c}(\kappa^{2}+k^{2})].\end{split}$
(14)

All the elements of these interaction matrices have a similar structure

$\mathcal{G}_{ab}=\mbox{single particle dynamics}+\mbox{mobility coefficient
for species }a\times\mbox{activity coefficient for species }b.$ (15)

The dynamical matrices in equations (14) are not symmetric, which means that
the effective interactions between the various species, following from the
different coupling with the chemical substrate, are nonreciprocal. Note that
while the homogeneous production/consumption rate $\alpha_{a}$ contributes at
all scales, the dipole contribution proportional to $\beta_{a}$ is expected to
be sub-leading in the macroscopic limit $k\rightarrow 0$, as it scales as
$k^{2}$.

In the next two sections, we present a stability analysis of the two-species
system in two particularly significant cases. In the first we consider the
fully coupled dynamics focusing mainly on the simpler case where the bare
translational and rotational diffusivities of the species are identical. In
this case, the dynamics reduce to the coupled dynamics of a conserved and a
non-conserved field which can undergo an instability to a spontaneously
oscillating state. In the second case, we assume that the timescales
$D_{r,a}^{-1}$ are smaller than all other timescales in the system, such that
the polarization fields can be enslaved to the density ones to obtain the
equations for two active densities, representing a microscopic realization of
the Nonreciprocal Cahn-Hilliard model [21].

##  III Coupled dynamics of two species

Here we simplify our analysis by restricting to the case of two species of
Janus colloids. This represents the simplest, yet physically relevant, system
of Janus particles that display nonreciprocal interactions due to phoretic
coupling between similar types of fields - the density fields or the polar
fields. To reduce the large parameter space, we first assume equal diffusivity
$D_{1}=D_{2}=D$, orientational noise timescale $D_{r,1}=D_{r,2}=D_{r}$ for the
two species, and equal self-propelling velocity $v_{1}=v_{2}=v$, while
allowing the phoretic mobility coefficients to be different. The dynamics of
the system are described by the spectral properties of the full dynamical
matrix $\mathcal{G}$, which is here defined via equation (14) as

$\mathcal{G}(k)=\left(\begin{array}[]{cc}k^{2}\,\mathcal{G}_{ab}^{tt}(k)&ik\,\mathcal{G}_{ab}^{tr}(k)\\\\[5.0pt]
ik\,\mathcal{G}_{ab}^{rt}(k)&\mathcal{G}_{ab}^{rr}(k)\\\ \end{array}\right).$
(16)

In particular, the eigenvalues of $\mathcal{G}$ quantify the rate with which
the different modes in the linearized system of equations grow or decay
exponentially with time. In addition, complex eigenvalues imply an oscillatory
response to perturbations. Here we list the four eigenvalues $\Lambda_{i}(k)$
with $i=1-4$ of the $\mathcal{G}$ matrix. The first pair of eigenvalues
$\Lambda_{1,2}$ are

$\Lambda_{1,2}(k)=-\left(D_{r}+k^{2}D\pm\sqrt{\Delta_{1}(k)}\right),\quad\quad\Delta_{1}(k)=\left(D_{r}\right)^{2}-\frac{k^{2}v^{2}}{3},$
(17)

and they are associated with stable modes that are independent of the phoretic
effects. The second pair of eigenvalues $\Lambda_{3,4}$ are

$\displaystyle\Lambda_{3,4}(k)=-\left(D_{r}+Dk^{2}+\frac{K(k)\bar{\rho}}{2}\frac{\bm{\alpha}\cdot\bm{\mu}-\bm{\beta}\cdot\bm{\Omega}}{D_{c}}\pm\sqrt{\Delta_{2}(k)}\right),$
(18)
$\displaystyle\Delta_{2}(k)=\left(\bar{\rho}K(k)\frac{\bm{\alpha}\cdot\bm{\mu}-\bm{\beta}\cdot\bm{\Omega}}{2D_{c}}\right)^{2}+\frac{\bar{\rho}K(k)}{D_{c}}\left[v\left(\bm{\alpha}\cdot\bm{\Omega}+\frac{k^{2}}{3}\bm{\beta}\cdot\bm{\mu}\right)-D_{r}\left(\bm{\beta}\cdot\bm{\Omega}+\bm{\alpha}\cdot\bm{\mu}\right)\right]+\Delta_{1}(k),$

where we define the screening parameter $K(k)=k^{2}/(k^{2}+\kappa^{2})$, total
average density $\bar{\rho}=\bar{\rho}_{1}+\bar{\rho}_{2}$, and we group the
phoretic coefficients for the two species as pairs of numbers which we write
in a compact manner adapting the vector notation
$\bm{\alpha}=(\alpha_{1},\alpha_{2})$, $\bm{\beta}=(\beta_{1},\beta_{2})$,
$\bm{\mu}=(\bar{\rho}_{1}\mu_{1}\bar{\rho}^{-1},\bar{\rho}_{2}\mu_{2}\bar{\rho}^{-1})$,
$\bm{\Omega}=(\bar{\rho}_{1}\Omega_{1}\bar{\rho}^{-1},\bar{\rho}_{2}\Omega_{2}\bar{\rho}^{-1})/3$.
The factor $K(k)$ is determined by the relative magnitude of $k/\kappa$ \- at
fixed $\kappa$ it vanishes when $k\ll\kappa$ and approaches unity when
$k\gg\kappa$. We recall that the linear regime relies on the condition of
steady chemical density
$\bar{c}=(\bar{\rho}_{1}\alpha_{1}+\bar{\rho}_{2}\alpha_{2})/(D_{c}\kappa^{2})\geq
0$, which is obtained by the balance of net positive production of the
substrate via the activity of the Janus colloids and its spontaneous
degradation, expressed by $\kappa>0$. Contrarily to the eigenvalues
$\Lambda_{1,2}$, the pair $\Lambda_{3,4}$ depends on the phoretic
coefficients. This implies that the mode structure of the linearized dynamics
splits into two parts. The first set of modes $\Lambda_{1,2}$ describes the
linear response of a self-propelled polar field coupled to a density field and
no other source of activity or interaction. On the other hand, $\Lambda_{3,4}$
receives contributions from groups of phoretic parameters whose relative
magnitudes can be tuned to give rise to linear instabilities. The reason
behind this partition at the linear level is the consideration of a single
chemical species $c$ which leads to the factorized structure in the dynamical
matrix discussed in equations (14) and (15). We will show later in this
section that this decomposition does not hold if the diffusivities and the
self-propelling velocities of the two species are unequal. For $k\ll\kappa$,
as $K(k)$ approaches the ratio $k^{2}/\kappa^{2}$, and the pair
$\Lambda_{3,4}$ approaches $\Lambda_{1,2}$, i.e. phoretic effects are
suppressed in dynamics occurring at lengthscales much larger than
$\kappa^{-1}$.

Figure 2:  Phase diagrams for $D_{1}=D_{2}=D=1$, $v_{1}=v_{2}=v=2.5$ for two
different values of $D_{r,1}=D_{r,2}=D_{r}$, and $k/\kappa=1.5$. The magenta
area is associated with a homogeneous phase. The dashed blue line identifies
exceptional points: it separates real and complex oscillating modes, that, if
unstable, could lead to “Spontaneous Oscillations” (blue area). In the
“Unstable I” region (orange area) the conserved mode is unstable and the non-
conserved mode is stable; it becomes unstable in the “Unstable II” phase. (a)
corresponds to small values of $D_{r}$. Assuming $\alpha_{a}>0$ and
$\beta_{a}>0$, the homogeneous phase is stable for
$\bm{\alpha}\cdot\bm{\mu}>0$, i.e., for effective repulsive interaction among
particles. For $\bm{\alpha}\cdot\bm{\mu}<0$ and
$\bm{\beta}\cdot\bm{\Omega}<0$, the colloids experience attractive interaction
and re-orient away from high-density regions: this could lead to a homogeneous
phase whenever the latter dominates the former, to spontaneous oscillations if
they are counterbalanced, or to the creation of asters otherwise. For
$\bm{\beta}\cdot\bm{\Omega}>0$ and $\bm{\alpha}\cdot\bm{\mu}<0$, both
chemotactic drift and alignment tend to aggregate particles, thus leading to
an instability of the density mode. In (b) $D_{r}$ takes larger values. The
main difference from (a) consists in the fact that a locally ordered phase can
arise only for strong enough aligning interaction
$\bm{\beta}\cdot\bm{\Omega}>0$ compared with $D_{r}$.

We will now discuss the instabilities indicated by $\Lambda_{i}$ and how they
can be tuned. We keep in mind that the behavior at small $k$ has the most
predictive power in determining the nonequilibrium steady states exhibited by
the system. The eigenvectors corresponding to $\Lambda_{i}$ contain
information about the combination of fields whose perturbations show
exponential growth. However, the complexity of the problem allows us to make
only the following comment. At vanishing $k$, the system is described by two
conserved modes ($\Lambda_{1,3}$) (a combination of the two density fields)
and two non-conserved modes ($\Lambda_{2,4}$) (combination of polarization
fields). This point will be clear from the Taylor expansion in $k$. An
exponential growth in the conserved modes could lead to active phase
separation. The growth of fluctuations in the non-conserved modes, presumably
leading to a growth in the longitudinal part of $\bm{p}$ leads to a relatively
unexplored type of orientational order called asters [28].

We start by showing that $\Lambda_{1,2}$ is always stable at small $k$. Taylor
expanding $\Lambda_{1,2}$ we find

$\displaystyle\Lambda_{1}(k)=-2D_{r}-k^{2}\left(D-\frac{v^{2}}{6D_{r}}\right)+O(k^{4}),$
(19) $\displaystyle\Lambda_{2}(k)=-k^{2}D^{\rm eff}(0)+O(k^{4}),$

where we introduce the effective diffusivity modified by dynamics of the
polarization field as for an active Brownian particle, $D^{\rm eff}(k)\equiv
D+v^{2}/[3(2D_{r}+Dk^{2})]$ [53, 54, 55, 56]. $D^{\rm eff}(k)$ is always
greater than zero meaning that density perturbations eventually decay. Note
however that the polarization field can show patterning at values of $k\approx
2D_{r}(D-v^{2}/6D_{R})^{-1}$.

Figure 3: Here we present two instances of the $\kappa$-dependence of the
eigenvalues $\Lambda_{3,4}$ in equation (18) for the two-species case. In (a),
corresponding to the small-$v$ limit, as $\kappa$ decreases, the system
displays oscillations at decreasing values of $k$. The parameters are here
chosen to be $D=D_{r}=1$ $D_{c}=0.5$ $v=0.1$ $\rho_{1}=5$, $\rho_{2}=1$,
$\alpha_{1}=1.5$, $\alpha_{2}=2.5$, $\beta_{1}=0.5$, $\beta_{2}=0.25$,
$\mu_{1}=1$, $\mu_{2}=2$, $\Omega_{1}=0.5$, $\Omega_{2}=1.4$. In (b),
corresponding to the large $v$-limit, as the effect of screening increases
oscillations occur at larger values of $k$. The associated parameters are
$D=D_{r}=0.1$ $D_{c}=1$ $v=30$ $\rho_{1}=5$, $\rho_{2}=1$, $\alpha_{1}=1.5$,
$\alpha_{2}=2.5$, $\beta_{1}=2.5$, $\beta_{2}=5.25$, $\mu_{1}=1$, $\mu_{2}=2$,
$\Omega_{1}=5.5$, $\Omega_{2}=5$.

###  III.1 Instabilities in $\Lambda_{3,4}$

The eigenvalues $\Lambda_{3,4}$ depend on all possible pairings of phoretic
mobilities and chemical activity which are four in total.
$\bm{\alpha}\cdot\bm{\mu}$ and $\bm{\beta}\cdot\bm{\mu}$ arise from particle
drift in response to isotropic and anisotropic chemical production
respectively. Similarly, $\bm{\alpha}\cdot\bm{\Omega}$ and
$\bm{\beta}\cdot\bm{\Omega}$ quantify alignment with substrate gradient
whether isotropic and anisotropic respectively. We first show a stability
diagram spanned by $K(k)\bm{\alpha}\cdot\bm{\mu}/D_{c}$ and
$K(k)\bm{\beta}\cdot\bm{\Omega}/D_{c}$; the details about the construction of
the phase diagram are discussed in appendix B. The phase diagram assumes a
fixed value of $k$, however, its structure is qualitatively the same for a
generic choice of the wave vector. Both $\Lambda_{3,4}$ are positive in the
“Unstable II” green area in figure 2, corresponding to an instability in both
the density and the polar fields, presumably a state with phase separation
coupled with aster formation. The “Unstable I” orange area refers to the
instability of the conserved mode only, resulting in phase separation driven
by phoretic activity. These two unstable regions are separated by the line
along which $\Lambda_{3}=0$, and $\Lambda_{4}<0$. Here, $\Lambda_{3,4}$ are
complex numbers for $\Delta_{2}<0$, which can happen only when
$(\bm{\alpha}\cdot\bm{\mu}-\bm{\beta}\cdot\bm{\Omega})<0$. The homogeneous
phase is “Stable” in the magenta region of figure 2 case of incoherent self-
propulsion ($\bm{\beta}\cdot\bm{\Omega}<0$) and chemotactic repulsion
($\bm{\alpha}\cdot\bm{\mu}>0$). The boundary (dashed blue line) between the
green and the blue regions is made of exceptional points: they separate the
unstable complex, associated with “Spontaneous Oscillations”, and unstable
real modes. The blue and the magenta regions in Fig 2 are separated by the
pale yellow line, where $\Lambda_{3,4}$ are purely imaginary: this purely
oscillatory state is strongly affected by the nonlinearities in equations (5)
and (6).

The modes $\Lambda_{3,4}$ possess a complex structure that is visible in the
myriad effects that emerge at different lengthscales. For finite $\kappa$ and
vanishing $k$, we Taylor expand the non-conserving eigenmode as
$\Lambda_{3}=-2D_{r}-k^{2}D^{\prime}+O(k^{4})$, with a modified diffusion
coefficient

$D^{\prime}=\frac{\bar{\rho}\,\bm{\Omega}}{\kappa^{2}D_{c}}\cdot\left(\frac{v}{2D_{r}}\bm{\alpha}-\bm{\beta}\right)-\frac{v^{2}}{6D_{r}}+D.$
(20)

An instability at $k\to 0$ is ruled out by $D_{r}$, which is constrained to be
positive. However, $D^{\prime}$ can change sign signaling an instability at
finite $k\approx\sqrt{|D_{r}/D^{\prime}|}$. At low $k$, the conserved mode is
$\Lambda_{4}=-k^{2}D^{\prime\prime}+O(k^{4})$, where the modified diffusion
constant is

$\displaystyle D^{\prime\prime}=D^{\rm
eff}(0)+\frac{\bar{\rho}\,\bm{\alpha}}{\kappa^{2}D_{c}}\cdot\left(\bm{\mu}-\frac{v}{2D_{r}}\bm{\Omega}\right),$
(21)

and it can change sign leading to active phase separation when the effect of
the phoretic interaction prevails over diffusion. $D^{\prime\prime}$ retains
its positive sign if the combinations $\bm{\alpha}\cdot\bm{\mu}$ and
$\bm{\alpha}\cdot\bm{\Omega}$ are positive and negative respectively -
significant departures can cause $D^{\prime\prime}$ to flip sign and signal an
instability. The first contribution in $D^{\prime\prime}$ arises from
interactions between colloids when each acts as a point source of chemicals.
The interactions are analogous to screened electrostatic ones, where the
positive (respectively, negative) sign of $\mu_{a}$ determines whether the
interactions are repulsive (attractive). This scenario is in contrast with the
expected phenomenology, where colloids interact via long-ranged chemical
fields leading to gravitational collapse or electrostatic screening. [57].
This contribution can be understood as the screened analog of a Keller-Segel-
like interaction, which was reported for the single species case in [45]. The
second term is a combination of the collective turning of the polarization to
point towards a local accumulation of the substrate and consequent drift in
that direction due to self-propulsion; for single-species case, see [45] and
appendix C. To summarize, at the smallest values of $k$, only the conserving
mode can trigger instability.

Complex modes arise when the discriminant in equation (18) becomes negative,
i.e., $\Delta_{2}<0$. As $\Delta_{2}(0)=D_{r}^{2}$ is a positive quantity, it
can change sign only at finite $k$, meaning that the response is oscillatory
only at finite lengthscales. These exceptional points, at which the complex
eigenvalues emerge, are given by the roots of the cubic equation
$\Delta_{2}(k^{2})=0$. They are denoted by the dashed line in the phase
diagram in figure 2. As seen in the expression for $\Delta_{2}$ in equation
(18), if $v$ is negligible, it is approximately a function of $k/\kappa$, and
not just $k$. This means that the value of $k$ where the exceptional points
appear scales with $\kappa$, approaching $k\to 0$ as $\kappa\to 0$, i.e. for
long-range interactions. This effect can be visualized in figure 3 (a), where
we show how reducing the value of $\kappa$, the eigenvalues - in this case,
associated with a stable disordered phase - become complex at smaller values
of $k$. As $\kappa\to 0$, $\Delta_{2}$ approaches a finite value such that the
model predicts (almost) global oscillations with frequency
$\sqrt{|\Delta_{2}|}$, for example at $k=\kappa$

$\Delta_{2}(\kappa)=\bar{\rho}^{2}\left(\frac{\bm{\alpha}\cdot\bm{\mu}-\bm{\beta}\cdot\bm{\Omega}}{4D_{c}}\right)^{2}+D_{r}^{2}+\frac{\bar{\rho}}{2D_{c}}\left[v\bm{\alpha}\cdot\bm{\Omega}-D_{r}\left(\bm{\beta}\cdot\bm{\Omega}+\bm{\alpha}\cdot\bm{\mu}\right)\right]+O(\kappa^{2}).$
(22)

In this regime, complex eigenvalues emerge when the combination
$(\bm{\beta}\cdot\bm{\Omega}+\bm{\alpha}\cdot\bm{\mu})>0$, ensuring that
$\Delta_{2}<0$. Even in the absence of chemotactic alignment ($\Omega_{a}=0$),
we obtain an (unstable) oscillating phase for $\bm{\alpha}\cdot\bm{\mu}$
sufficiently negative. Being $\Lambda_{4}$ associated with the non-conserved
mode in the system, we deduce that an interplay of two fields is essential for
oscillations. The oscillations are spontaneous if the real part of
$\Lambda_{3,4}$ is negative, which is given according to equation (18) by

$\mbox{Re}(\Lambda_{3,4})=-D_{r}-Dk^{2}+\frac{K(k)\bar{\rho}}{D_{c}}(\bm{\alpha}\cdot\bm{\mu}-\bm{\beta}\cdot\bm{\Omega}).$
(23)

The condition $\mbox{Re}(\Lambda_{3,4})=0$ \- the green line in figure 2 \-
corresponds to the threshold of a Hopf bifurcation. Orientational diffusion
stabilizes the homogeneous disordered phase at $k=0$, hence spontaneous
oscillations also occur only at finite $k$ and for
$(\bm{\beta}\cdot\bm{\Omega}-\bm{\alpha}\cdot\bm{\mu})<0$. We expect a regime
that should be dominated by phase separation as well as aster formation. If
the contribution of $D$ and $v$ in $\Lambda_{3,4}$ are non-negligible, the
eigenmodes $\Lambda_{3,4}$ depend both on $k$ and $\kappa$ creating the
possibility of oscillatory instabilities at finite $k$ that are controlled by
$v$ and the combination
$\left(\bm{\alpha}\cdot\bm{\Omega}+{k^{2}}\bm{\beta}\cdot\bm{\mu}/3\right)$.
In this regime of parameters, - as we showcase in figure 3 (b) - as the
screening parameter $\kappa$ decreases, oscillations may emerge at increasing
values of $k$.

As a concluding remark, we mention that the eigenvalues $\Lambda_{3,4}$ are
formally similar to the one-species case; we refer to appendix C for the
details. This similarity becomes a qualitative equivalence whenever it is
possible to factor out the scalar product of the vectorial production rates
and phoretic interactions $\bm{\alpha}$, $\bm{\beta}$, $\bm{\mu}$ and
$\bm{\Omega}$. For instance, if $\alpha_{1}=\alpha_{2}$ and
$\beta_{1}=\beta_{2}$ we have $\bm{\alpha}\cdot\bm{\mu}=\alpha\mu_{\rm eff}$,
$\bm{\alpha}\cdot\bm{\Omega}=\alpha\Omega_{\rm eff}$,
$\bm{\beta}\cdot\bm{\mu}=\beta\mu_{\rm eff}$, and
$\bm{\beta}\cdot\bm{\Omega}=\beta\Omega_{\rm eff}$ such that $\mu_{\rm eff}$
and $\Omega_{\rm eff}$ can be interpreted as one species effective phoretic
couplings.

For the general case $D_{r,1}\neq D_{r,2}$, $D_{1}\neq D_{2}$ and $v_{1}\neq
v_{2}$ all four eigenvalues may be unstable, and the picture becomes richer.
We represent an instance of this behavior in figure 4 (a), which displays the
region of the phase diagram where eigenvalues are all unstable for a given
value of $k$. To fix ideas, if we set $\alpha_{a},\beta_{a}>0$, all the modes
are unstable for sufficiently high and negative values of $\mu_{a}$ and
$\Omega_{a}$: particles are expected to aggregate because of chemotaxis, while
their orientation is expected to point away from the induced high-density
region. These facts, together with diffusion, now occurring over different
timescales for the two species, lead to a strongly unstable phase. In figure
4(b)-(c), we plot the four eigenvalues of the dynamical matrix as a function
of the wave vector $k$ for two specific choices of the parameter values of the
phase diagram in figure 4 (a).

To summarize the main points in this section, we find that phase separation
driven by phoretic interactions emerges as the dominant behavior at the
largest lengthscales. Oscillations, which can also be spontaneous, occur
generically at finite wave numbers. At precisely $\kappa=0$, both the
conserving and the non-conserving modes can be unstable. In this regime, we
truly have two hydrodynamic modes: one corresponding to number conservation,
and the other an order parameter field that can undergo a phase transition and
pick up a finite value.

Figure 4: In panel (a) we show the phase diagram for generic values of the
non-active parameters $D$, $D_{r}$, and $v$. In this case, the system may
display four unstable modes. For instance, in the diagram of panel (a): the
yellow area \raisebox{-.9pt} {2}⃝ corresponds to two complex unstable and real
unstable eigenvalues, while the green one \raisebox{-.9pt} {3}⃝ to four
complex unstable eigenvalues. The violet area \raisebox{-.9pt} {1}⃝ indicates
the region of the parameters where at least one eigenvalue is stable. The
parameters are given by $D_{1}=D_{2}=1$, $D_{r,1}=1$, $D_{r,2}=1.2$,
$v_{1}=10$, $v_{2}=20$, and $k/\kappa=2.5$ . In panels (b) and (c) we
represent the real and imaginary parts of the four eigenvalues of the
dynamical matrix for the parameters corresponding to the bullet and star
symbol in the phase diagram as a function of $k$. The vertical black line
identifies the value of $k/\kappa$ of the phase diagram in the (a) panel.

##  IV Adiabatic elimination of the polarization fields: coupled active
densities

As discussed in the previous section, the dynamical properties of the two
species of Janus particles are described by the conserved slow density mode
$\rho_{a}(\bm{r},t)$ and the longitudinal part of the non-conserved
polarization field ${p}_{\parallel,a}(\bm{r},t)$. Strictly speaking, $\delta
p_{\parallel,a}$ is a slow variable at lengthscales that are small compared to
the screening length $\kappa^{-1}$. While probing the dynamics at the largest
lengthscales, and for finite $\kappa^{-1}$, we can assume that $\delta
p_{\parallel,a}$ relaxes faster than $\delta\rho_{a}$ and express the former
in terms of the latter. Equivalently, the field $\delta
p_{\parallel,a}(\bm{k},t)$ in equation (13) adapts instantaneously to the time
variations of $\delta\rho_{a}(\bm{k},t)$ according to

$\displaystyle\delta p_{\parallel,a}(\bm{k},t)$
$\displaystyle=-ik\sum_{b}\Gamma_{ab}(k)\delta\rho_{b}(\bm{k},t).$ (24)

The matrix $\Gamma$ appearing in equation (24) is defined in terms of those
appearing in equation (14) as
$\Gamma(k)\equiv\left[\mathcal{G}^{rr}(k)\right]^{-1}\mathcal{G}^{rt}(k)$.
Explicitly, its entries are given by

$\displaystyle\Gamma_{ab}(k)\equiv\frac{2}{3\bar{D}_{a}(k)}\left(\frac{v_{a}}{2}\delta_{ab}-\bar{\rho}_{a}\Omega_{a}\Pi_{b}(k)\right),$
(25)

where $\bar{D}_{a}(k)=2D_{r,a}+D_{a}k^{2}$, and the redefined substrate
production rate $\Pi_{a}(k)$ reads

$\Pi_{a}(k)\equiv\frac{1}{\kappa^{2}D_{c}+k^{2}\left[D_{c}-2\left(\bar{\rho}_{1}{\beta_{1}\Omega_{1}}/{\bar{D}_{1}}(k)+\bar{\rho}_{2}{\beta_{2}\Omega_{2}}/{\bar{D}_{2}(k)}\right)/3\right]}\left(\alpha_{a}-k^{2}\frac{v_{a}\beta_{a}}{3\bar{D}_{a}(k)}\right).$
(26)

The matrix $\Gamma$ quantifies the linear response matrix connecting the polar
field to variations of the density, while $\Pi_{a}(k)$ is the effective rate
of substrate production by species $a$.

Figure 5: Phase diagram from the linear stability analysis. The magenta area
identifies the values of $\bar{\rho}_{1}M_{1}\Pi_{1}$ and
$\bar{\rho}_{2}M_{2}\Pi_{2}$ corresponding respectively to a linear stable
regime. The “Unstable I” (orange) region corresponds to the case where the
system presents one unstable and one stable mode. In the “Unstable II” (green)
one, the modes are still real but both unstable. In the correspondence of the
blue region of the phase diagram, “Spontaneous Oscillations” arise in the
system, meaning the two eigenvalues are complex conjugate with positive real
part (unstable). The parabola separates the parameter space between real and
complex eigenvalues.

We now investigate the behavior of the system in this regime, where the only
relevant fields are the densities of the two species of particles. A necessary
condition for this approximation to hold is to have positive eigenvalues of
$\mathcal{G}^{rr}(k)$, thus ensuring the relaxation of the polar field in
equation (13). The enslaving fails when the determinant of $\mathcal{G}^{rr}$
vanishes, which happens when one or both the eigenvalues of $\mathcal{G}^{rr}$
are null. Large enough $D_{r,a}$ rules out this possibility, ensuring the
validity of the enslaving. As discussed in the previous section, the
approximation holds also for sufficiently large $\kappa^{2}$, ruling out aster
condensation. The denominator of $\Pi_{a}(k)$ is proportional to the
determinant of $\mathcal{G}^{rr}$, which is positive whenever the relaxation
approximation holds. Therefore, the sign of $\Pi_{a}(k)$ follows from those of
$\alpha_{a}$, $\beta_{a}$, the wave vector $k$, and their relative amplitude.
Note that in the case of uniform production or consumption of the chemical
$\beta_{a}=0$, the fast relaxation approximation is well defined and the
effective production rate reduces to
$\Pi_{a}(k)=\alpha_{a}/[D_{c}(k^{2}+\kappa^{2})]$, whose sign depends only on
that of $\alpha_{a}$. Substituting the expression of the polar field in
equation (24) in the equation for $\delta\rho_{a}(\bm{k},t)$ in equation (13),
one finds

$\displaystyle\partial_{t}\delta\rho_{a}(\bm{k},t)=k^{2}\sum_{b}\mathcal{G}_{ab}(k)\,\delta\rho_{b}(\bm{k},t),$
(27)

where we have introduced the effective diffusion matrix
$\mathcal{G}(k)\equiv\mathcal{G}^{tt}(k)+\mathcal{G}^{tr}(k)\,\Gamma(k)$. One
can show that the entries of $\mathcal{G}(k)$ can be expressed as

$\mathcal{G}_{ab}(k)=-\delta_{ab}D^{\rm
eff}_{a}(k)-\bar{\rho}_{a}M_{a}(k)\Pi_{b}(k),$ (28)

where, we define an effective mobility $M_{a}$ and an effective diffusivity
$D_{a}^{\rm eff}$ as

$M_{a}(k)\equiv\mu_{a}-\frac{2}{3}\frac{v_{a}\Omega_{a}}{\bar{D}_{a}(k)},\quad
D_{a}^{\rm eff}(k)=D_{a}+\frac{v_{a}^{2}}{3\bar{D}_{a}(k)}.$ (29)

Interestingly enough, the matrix $\mathcal{G}(k)$ inherits the same factorized
structure as its building blocks in equations (14), i.e., a diagonal
contribution coming from diffusion and self-propulsion, and contributions are
multiples of phoretic mobility and chemical production rate as in equation
(15). Note that, while $D_{a}^{\rm eff}$ is always positive, the sign of the
effective mobility $M_{a}$ depends on the relative amplitudes of the phoretic
drift and alignment interactions. By splitting the symmetric and antisymmetric
parts of the dynamical matrix $\mathcal{G}$ as

$\mathcal{G}=-\left(\begin{array}[]{cc}D_{1}^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}&\bar{\rho}_{1}M_{1}\Pi_{2}\\\\[5.0pt]
\bar{\rho}_{2}M_{2}\Pi_{1}&D_{2}^{\rm eff}+\bar{\rho}_{2}M_{2}\Pi_{2}\\\
\end{array}\right)=\left(\begin{array}[]{cc}d_{1}&\chi+\psi\\\\[5.0pt]
\chi-\psi&d_{2}\\\ \end{array}\right),$ (30)

one can identify self-interaction species term $d_{a}=\mathcal{G}_{aa}$, the
symmetric (reciprocal) contribution to the interaction between the two
different species $\chi=(\mathcal{G}_{12}+\mathcal{G}_{21})/2$, and its anti-
symmetric (nonreciprocal) one $\psi=(\mathcal{G}_{12}-\mathcal{G}_{21})/2$. By
its definition, $\chi$ is symmetric under the exchange of the coefficients of
the two species, while $\psi$ reverses sign under the same transformation. At
leading order of a small-$k$ expansion the non-reciprocity parameter $\psi(k)$
reads

$\displaystyle\psi(k)=$ $\displaystyle\,\psi_{0}+k^{2}\psi_{2}+O(k^{2}),$ (31)
$\displaystyle\psi_{0}=$
$\displaystyle-\frac{1}{2\kappa^{2}D_{c}}\left[\left(\bar{\rho}_{1}\mu_{1}\alpha_{2}-\bar{\rho}_{2}\mu_{2}\alpha_{1}\right)-\left(\bar{\rho}_{1}\alpha_{2}\frac{v_{1}\Omega_{1}}{3D_{r,1}}-\bar{\rho}_{2}\alpha_{1}\frac{v_{2}\Omega_{2}}{3D_{r,2}}\right)\right],$
$\displaystyle\psi_{2}=$
$\displaystyle-\frac{1}{12\kappa^{2}D_{c}}\left[\frac{\bar{\rho}_{1}v_{1}\Omega_{1}D_{1}\alpha_{2}}{(D_{r,1})^{2}}-\frac{\bar{\rho}_{2}v_{2}\Omega_{2}D_{2}\alpha_{1}}{(D_{r,2})^{2}}-\left(\bar{\rho}_{1}\mu_{1}\frac{v_{2}\beta_{2}}{D_{r,2}}-\bar{\rho}_{2}\mu_{2}\frac{v_{1}\beta_{1}}{D_{r,1}}\right)+\frac{v_{1}v_{2}}{3D_{r,1}D_{r_{2}}}\left(\bar{\rho}_{1}\Omega_{1}\beta_{2}-\bar{\rho}_{2}\Omega_{2}\beta_{1}\right)\right]$
$\displaystyle-\frac{\psi_{0}}{\kappa^{2}D_{c}}\left(D_{c}-\bar{\rho}_{1}\frac{\beta_{1}\Omega_{1}}{3D_{r,1}}-\bar{\rho}_{2}\frac{\beta_{2}\Omega_{2}}{3D_{r,2}}\right),$

where the first term identifies the macroscopic effective nonreciprocal
interaction between the two different species. The parameter $\psi_{0}$ has
been introduced in [21, 22] in a minimal model for nonreciprocal interaction
between multiple species, called the nonreciprocal Cahn Hilliard (NRCH) model.
The introduction of this term in a phase-separating system of many conserved
densities leads to arrested phase separation, broken spatial parity, and
broken time-reversal symmetry, producing traveling waves and patterns. The
second term couples the two species at the fourth order in gradients and can
be interpreted as a nonreciprocal surface tension [58]. Interestingly enough,
for $\Omega_{a}=\beta_{a}=0$, the expression of $\psi$ simplifies to

$\psi=-\frac{\bar{\rho}_{1}\mu_{1}\alpha_{2}-\bar{\rho}_{2}\mu_{2}\alpha_{1}}{2D_{c}(\kappa^{2}+k^{2})},$
(32)

indicating that uniformly coated colloids interacting via a screened chemical
field generates nonreciprocal couplings and all orders in gradients, the two
nonreciprocal at lowest order in gradients are related simply as
$\psi_{2}=-\psi_{0}\kappa^{-2}$. Non-zero $\beta_{a}$ and $\Omega_{a}$ take us
away from the simple relation between $\psi_{0}$ and $\psi_{2}$, which can now
be tuned independently of one another.

The self interaction $d_{a}$ and the reciprocal interaction $\chi$ can be
expanded similarly as $d_{a}=d_{a,0}+d_{a,2}k^{2}$, and
$\chi=\chi_{0}+\chi_{2}k^{2}$. The coefficient $d_{a,0}$ is the strength of
self-interaction for species $a$, while $\chi_{0}$ is the effective reciprocal
interaction. The terms occurring at higher orders in gradients, namely
$d_{2,a}$, and $\chi_{2}$ are the coefficients for interfacial tension. Recall
that the eigenmodes in the previous section had contributions from symmetric
combinations of the phoretic parameters which could be written as dot products
such as $\bm{\alpha}\cdot\bm{\mu}$. Similar simple relations hold for the
elements of $\mathcal{G}$. We can express $d_{a}$, $\psi$, and $\chi$
compactly by introducing the following matrices

$\bm{\sigma}_{1}=\begin{pmatrix}0&1\\\
-1&0\end{pmatrix},\;\;\;\bm{\sigma}_{2}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\;\;\;\bm{\sigma}_{3}=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},$ (33)

and the vectors of parameters

$\bm{\Omega}^{\prime}=\left(\frac{v_{1}\bar{\rho}_{1}\Omega_{1}}{\bar{\rho}D_{r,1}},\frac{v_{2}\bar{\rho}_{2}\Omega_{2}}{\bar{\rho}D_{r,2}}\right),\;\;\;\;\bm{\beta}^{\prime}=\left(\frac{v_{1}\beta_{1}}{D_{r,1}},\frac{v_{2}\beta_{2}}{D_{r,2}}\right),\;\;\;\;\bm{\Omega}^{\prime\prime}=\left(\frac{\bar{\rho}_{1}D_{1}\Omega_{1}}{\bar{\rho}(D_{r,1})^{2}},\frac{\bar{\rho}_{2}D_{2}\Omega_{2}}{\bar{\rho}(D_{r,2})^{2}}\right),\;\;\text{and}\;\;\bm{\Omega}^{\prime\prime\prime}=\left(\frac{\bar{\rho}_{1}\Omega_{1}}{\bar{\rho}D_{r,1}},\frac{\bar{\rho}_{2}\Omega_{2}}{\bar{\rho}D_{r,2}}\right),$
(34)

and the auxiliary scalar functions

$A_{0}(\mathcal{M})=-\frac{\bar{\rho}}{2\kappa^{2}D_{c}}\left(\bm{\mu}\cdot\mathcal{M}\cdot\bm{\alpha}-\frac{1}{3}\bm{\Omega}^{\prime}\cdot\mathcal{M}\cdot\bm{\alpha}\right)=-\frac{\bar{\rho}}{2\kappa^{2}D_{c}}\left(\bm{\mu}-\frac{1}{3}\bm{\Omega}^{\prime}\right)\cdot\mathcal{M}\cdot\bm{\alpha},$
(35)
$A_{2}(\mathcal{M})=-\frac{\bar{\rho}}{12\kappa^{2}D_{c}}\left(\bm{\Omega}^{\prime\prime}\cdot\mathcal{M}\cdot\bm{\alpha}-\bm{\mu}\cdot\mathcal{M}\cdot\bm{\beta}^{\prime}+\frac{v_{1}v_{2}}{D_{r,1}D_{r_{2}}}\bm{\Omega}\cdot\mathcal{M}\cdot\bm{\beta}\right)-\frac{A_{0}(\mathcal{M})}{\kappa^{2}D_{c}}\left(D_{c}-\bar{\rho}\bm{\beta}\cdot\bm{\Omega}^{\prime\prime\prime}\right).$
(36)

Finally, we can write the coefficients of the model as follows

$\begin{split}\psi_{0}&=A_{0}(\bm{\sigma}_{1}),\\\
\chi_{0}&=A_{0}(\bm{\sigma}_{2}),\\\
d_{1,0}&=-D_{1}^{\rm{eff}}(0)+A_{0}(\mathbb{I}+\bm{\sigma}_{3}),\\\
d_{2,0}&=-D_{2}^{\rm{eff}}(0)+A_{0}(\mathbb{I}-\bm{\sigma}_{3}),\end{split}\quad\quad\begin{split}\psi_{2}&=A_{2}(\bm{\sigma}_{1}),\\\
\chi_{2}&=A_{2}(\bm{\sigma}_{2}),\\\
d_{1,2}&=-D_{1}v_{1}^{2}/[12(D_{r,1})^{2}]+A_{2}(\mathbb{I}+\bm{\sigma}_{3}),\\\
d_{2,2}&=-D_{2}v_{2}^{2}/[12(D_{r,2})^{2}]+A_{2}(\mathbb{I}-\bm{\sigma}_{3}).\end{split}$
(37)

In equation (37), $\psi_{0,2}$ are pseudoscalars while the rest are scalar in
the exchange of the two species.

###  IV.1 Stability diagram

The linear stability of the homogeneous state is determined by the eigenvalues
$\Lambda_{1,2}(k)$ of $\mathcal{G}(k)$ which can be expressed compactly in
terms of the effective mobility $\bar{\rho}_{a}M_{a}$ and the production rate
$\Pi_{a}(k)$ as

$\displaystyle\Lambda_{1,2}$ $\displaystyle=-\left(\frac{D_{1}^{\rm
eff}+D_{2}^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}+\bar{\rho}_{2}M_{2}\Pi_{2}}{2}\pm\sqrt{\Delta}\right)$
(38) $\displaystyle\Delta$ $\displaystyle=\left(\frac{D_{1}^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}-D_{2}^{\rm
eff}-\bar{\rho}_{2}M_{2}\Pi_{2}}{2}\right)^{2}+\bar{\rho}_{1}M_{1}\Pi_{1}\bar{\rho}_{2}M_{2}\Pi_{2}.$

The growth rates of the eigenmodes are then given by $k^{2}\Lambda_{1,2}(k)$
which vanish as $k\to 0$, as a result of number conservation. For purely
reciprocal interactions corresponding to $\psi=0$, $\mathcal{G}$ is symmetric
and the eigenvalues are always real. Another case where the eigenvalues are
always real is for equal effective diffusivities $D_{1}^{\rm eff}=D_{2}^{\rm
eff}=D^{\rm eff}$. In this case, $\mathcal{G}$ is a rank one tensor of the
form

$\mathcal{G}=-D^{\rm
eff}\mathbb{I}-(\bar{\rho}_{1}M_{1},\bar{\rho}_{2}M_{2})^{\rm
T}(\Pi_{1},\Pi_{2}).$ (39)

The eigenmodes are $\Lambda_{1}=-D^{\rm eff}$, $\Lambda_{2}=-(D^{\rm
eff}+\operatorname{tr}{\mathcal{G}})=-(D^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}+\bar{\rho}_{2}M_{2}\Pi_{2})$.
$\Lambda_{1}(k)<0$ is the stable mode, while $\Lambda_{2}(k)$ can be positive
and trigger an instability for $D^{\rm
eff}>-(\bar{\rho}_{1}M_{1}\Pi_{1}+\bar{\rho}_{2}M_{2}\Pi_{2})$. For $D^{\rm
eff}_{1}\neq D^{\rm eff}_{2}$, and considering equation (38), it is clear that
complex eigenvalues can occur only if

$\chi^{2}-\psi^{2}=\bar{\rho}_{1}M_{1}\Pi_{1}\bar{\rho}_{2}M_{2}\Pi_{2}<0,$
(40)

whenever the nonreciprocal coupling exceeds its reciprocal counterpart; this
condition is sufficient for $d_{1}=d_{2}$. Henceforth, without any loss of
generality, we always assume that $\delta D=D_{1}-D_{2}>0$.

Similarly to section III, we first present the phase diagram in the plane of
parameters $\bar{\rho}_{a}M_{a}\Pi_{a}$. In figure 5, the region where the
eigenvalues are complex lies on the convex side of the parabola. The
homogeneous state is stable in the magenta part. Both $\Lambda_{1,2}$ are
unstable in the green region and signal an instability leading to phase
separation. In the orange region, one of the two eigenvalues leads to an
instability. In the purple region, the two eigenvalues are complex and
unstable leading to oscillating densities. The phase diagram has the same
topology for all values of $k$.

We will now discuss how the instabilities appear at different lengthscales,
i.e. as $k$ is varied. Alternatively, if one keeps all the other parameters
fixed and considers the eigenvalues as a function of $k$ only, the system
state can pass through different phases, described as a state curve
parameterized by $k$. Next, we will discuss the instabilities that occur as
$k\to 0$ and at finite $k$ and illustrate them on the phase diagram in figure
5 using state curves.

Figure 6: Type II of instability. Another type of instability that may occur
in our conserved system can be attained whenever the eigenvalues are already
unstable at small values of $k$, i.e., Re$(\Lambda_{2}(0))>0$. In (a) and (b)
we show, respectively, the real and imaginary parts of the most unstable
eigenvalue $\Lambda_{2}$ at the onset of the instability, which is triggered
by variations of the parameter $\mu_{2}$. In (c) the equivalent parametric
plot of the state curve described by the system in the phase diagram is
showcased: the initial point of the trajectory lies in the instability region
and eventually ends in the stable one. The parameters are given by: $D_{c}=1$,
$D_{1}=2$, $D_{2}=1.5$, $D_{r,1}=1$, $D_{r,2}=2$, $v_{1}=1$, $v_{2}=1.5$,
$\bar{\rho}_{1}=\bar{\rho}_{2}=1$, $\alpha_{1}=4$, $\alpha_{2}=0.5$,
$\beta_{1}=0.5$, $\beta_{2}=1$, $\mu_{1}=-1$, $\Omega_{1}=10$, $\Omega_{2}=5$.
As $k$ varies, the phase diagram slightly changes quantitatively but not
qualitatively, thus we draw it for a fixed value of $k$. In this figure, it
has been chosen $k/\kappa=10$.

###  IV.2 Instabilities at vanishing $k$

We expand $\Lambda_{1,2}$ in a Taylor series to obtain the following general
expression at zeroth order in $k$

$\Lambda_{1,2}(0)=\begin{cases}\frac{\operatorname{tr}\mathcal{G}(0)}{2}\mp\sqrt{\Delta(0)}&\text{for
}\Delta(0)>0\\\ \frac{\operatorname{tr}\mathcal{G}(0)}{2}\mp
i\sqrt{-\Delta(0)}&\text{for }\Delta(0)<0\end{cases},$ (41)

where
$\Delta(k)\equiv(\operatorname{tr}\mathcal{G}/2)^{2}-\operatorname{det}\mathcal{G}$.
The coefficients in equation (41) are

$\operatorname{tr}\mathcal{G}(0)=-\frac{1}{2}\left[D_{1}^{\rm
eff}(0)+D_{2}^{\rm
eff}(0)+\frac{\bar{\rho}\,\bm{\alpha}}{D_{c}\kappa^{2}}\cdot\left(\bm{\mu}-\bm{\Omega}\right)\right],$
(42)

and

$\displaystyle\Delta(0)=\left[\operatorname{tr}\mathcal{G}(0)\right]^{2}+\left(\mu_{1}-\frac{v_{1}\Omega_{1}}{3D_{r,1}}\right)\left(\mu_{2}-\frac{v_{2}\Omega_{2}}{3D_{r,2}}\right)\frac{\bar{\rho}_{1}\alpha_{1}\bar{\rho}_{2}\alpha_{2}}{(D_{c}\kappa^{2})^{2}}.$
(43)

A system-wide instability arises whenever Re$(\Lambda_{2})$ is positive,
resulting in bulk phase separation of the two species. At large $k$, the
eigenvalues are stable and real, ensuring the system’s stability at small
length scales. Accordingly, Re$\Lambda_{2}$ attains its maximum value at
$k_{c}$. As the instability is turned off, $k_{c}\rightarrow 0$ and the value
of the associated lengthscale $k_{c}^{-1}$ diverges since the maximum of the
most unstable eigenvalue vanishes, i.e., Re$\Lambda_{2}(k_{c}\rightarrow
0)\rightarrow 0$. This corresponds to a type II of instability according to
the Cross and Hohenberg classification [59]. If at the onset of the
instability $\Lambda_{2}$ is real, it is further classified as stationary,
while in its complex counterpart as oscillatory. We show an instance of this
type of instability in figure 6. In panels (a) and (b) we showcase the real
and imaginary parts of the most unstable mode $\Lambda_{2}$, and how
instability is triggered while varying the $\mu_{2}$ parameter. In this
specific case, we have an instance of stationary instability. However, being
$\Delta(0)$ a nonlinear combination of the phoretic parameters, an oscillatory
instability develops for a large part of the parameter space, e.g., for
$\mu_{a},\Omega_{a}>0$, $\Delta(0)<0$ if

$\operatorname{sign}(\mu_{1}-v_{1}\Omega_{1}/3D_{r,1})\neq\operatorname{sign}(\mu_{2}-v_{2}\Omega_{2}/3D_{r,2}),$
(44)

provided that both species produce chemicals, i.e., $\alpha_{a}>0$. The
associated system-wide oscillations occur with an angular frequency
$\sqrt{|\Delta(0)|}$. Another way to visualize the onset of the instability is
displayed in figure 6 (c), where we show how, by varying $k$, the parameters
$\bar{\rho}_{a}M_{a}(k)\Pi_{a}(k)$ describe a state curve in the phase
diagram, corresponding to the green curve in 6 (a). At $k=0$ the system is the
unstable region of parameters (orange area), and as the value of $k$ increases
it crosses the complex unstable region (blue area), then the complex stable
region (magenta area). Note that, as a consequence of diffusion, at large $k$
the system is always stable, as it can be seen from
$\lim_{k\to\infty}\bar{\rho}_{a}M_{a}(k)\Pi_{a}(k)=0.$ Note that as the value
of $k$ is varied, also the phase diagram changes. However, these changes are
barely perceivable compared to the $\bar{\rho}_{a}M_{a}\Pi_{a}$ ones. This
fact allows us to give a meaningful qualitative representation of the
$k$-parametrized state curve while drawing the diagram at a fixed value of
$k$.

###  IV.3 Finite wavelength instabilities

Instabilities in the system may also occur for Re$(\Lambda_{2}(0))<0$, i.e.,
the species are not separated at the macroscopic scale but produce patterns
with a specific lengthscale. To characterize this type of instability, it is
necessary to expand $\Lambda_{1,2}$ up to $k^{2}$

$\text{Re}(\Lambda_{1,2})=\begin{cases}\frac{\operatorname{tr}\mathcal{G}(0)}{2}\mp\sqrt{\Delta(0)}+\frac{k^{2}}{4}\left[(\operatorname{tr}\mathcal{G})^{\prime\prime}(0)\mp\frac{\Delta^{\prime\prime}(0)}{\sqrt{\Delta(0)}}\right]+O(k^{4})&\text{for
}\Delta(0)>0\\\
\frac{\operatorname{tr}\mathcal{G}(0)}{2}+\frac{k^{2}}{4}(\operatorname{tr}\mathcal{G})^{\prime\prime}(0)+O(k^{4})&\text{for
}\Delta(0)<0\end{cases},$ (45)

where ′ denotes the derivative with respect to $k$. The expressions for terms
that contribute at quadratic order in $k$ are given by

$\displaystyle(\operatorname{tr}\mathcal{G})^{\prime\prime}(0)=-2\left\\{-\frac{D_{1}v_{1}^{2}}{12(D_{r,1})^{2}}-\frac{D_{2}v_{2}^{2}}{12(D_{r,2})^{2}}+\frac{\bar{\rho}\,\bm{\alpha}}{D_{c}\kappa^{2}}\cdot\left[\frac{\bm{\Omega}^{\prime\prime}}{6}-\frac{1}{\kappa^{2}}\left(\bm{\mu}-\frac{\bm{\Omega}^{\prime}}{3}\right)\right]\right\\},$
(46) $\displaystyle\Delta^{\prime\prime}(0)=\left[D_{1}^{\rm
eff}(0)-D_{2}^{\rm
eff}(0)+\frac{\bar{\rho}}{D_{c}\kappa^{2}}\left(\bm{\mu}-\frac{\bm{\Omega}^{\prime}}{3}\right)\cdot\bm{\sigma}_{3}\cdot\bm{\alpha}\right]\left\\{-\frac{D_{1}v_{1}^{2}}{12(D_{r,1})^{2}}-\frac{D_{2}v_{2}^{2}}{12(D_{r,2})^{2}}+\frac{\bar{\rho}\,\bm{\alpha}}{D_{c}\kappa^{2}}\cdot\left[\frac{\bm{\Omega}^{\prime\prime}}{6}-\frac{1}{\kappa^{2}}\left(\bm{\mu}-\frac{\bm{\Omega}^{\prime}}{3}\right)\right]\right\\}$
$\displaystyle\;\;\;\;\;\;\;\;\;\;+2\,\frac{\rho_{1}\alpha_{1}\rho_{2}\alpha_{2}}{(D_{c}\kappa^{2})^{2}}\left\\{\left(\mu_{1}-\frac{v_{1}\Omega_{1}}{3D_{r,1}}\right)\frac{v_{2}\Omega_{2}D_{2}}{6(D_{r,2})^{2}}+\left(\mu_{2}-\frac{v_{2}\Omega_{2}}{3D_{r,2}}\right)\frac{v_{1}\Omega_{1}D_{1}}{6(D_{r,1})^{2}}-\frac{2}{\kappa^{2}}\left(\mu_{1}-\frac{v_{1}\Omega_{1}}{3D_{r,1}}\right)\left(\mu_{2}-\frac{v_{2}\Omega_{2}}{3D_{r,2}}\right)\right\\}.$

Note that the eigenmodes are invariant if the species indices $1$ and $2$ are
swapped. $\Delta^{\prime\prime}(0)$ is invariant as it is a product of two
quantities both of which reverse sign when under the swap
$1\longleftrightarrow 2$. The system is stabilized at the shortest lengths by
diffusive processes. A finite wavelength instability is triggered in the
system if the most unstable mode Re$(\Lambda_{2})$ becomes null at $k=k_{-}$
and then acquires a positive value. It reverses its sign again at $k=k_{+}$.
In this case, there exists a wave number $k_{c}$ intermediate between $k_{-}$
and $k_{+}$ where Re$(\Lambda_{2})$ attains its maximal value, physically
associated with pattern formation at the lengthscale $k_{c}^{-1}$. An
approximate expression for $k_{-}$ can be retrieved from equation (45) as

$k_{-}=\begin{cases}2\sqrt{-\frac{\operatorname{tr}\mathcal{G}(0)/2-\sqrt{\Delta(0)}}{(\operatorname{tr}\mathcal{G})^{\prime\prime}(0)-\Delta^{\prime\prime}(0)/\sqrt{\Delta}(0)}}&\text{for
}\Delta(0)>0\\\\[10.0pt]
\sqrt{-2\frac{\operatorname{tr}\mathcal{G}(0)}{(\operatorname{tr}\mathcal{G})^{\prime\prime}(0)}}&\text{for
}\Delta(0)<0,\end{cases}$ (47)

given that the square root exists.

At arbitrary $k$, this type of behavior is well represented in figure 7, where
we provide a few examples of the onset of the instability. Note that if the
eigenvalues are real at $k=k_{-}\text{ (respectively, }k_{+})$ then
$\Lambda_{2}=0$ and $\Lambda_{1}\leq 0(\geq 0)$, and the state curve crosses
the line corresponding to $\operatorname{det}\mathcal{G}=0$ at $k_{-}\text{
(respectively, }k_{+})$ from its stable (unstable) side. This case is
represented in figures 7 (b) and 8 (a). If the eigenvalues are complex
conjugate at $k_{\pm}$, then $\Lambda_{1/2}(k_{\pm})=\pm
i\text{Im}(\Lambda_{1/2}(k_{\pm}))$ and the state curve crosses the line
corresponding to $\operatorname{tr}\mathcal{G}=0$; this case is shown in
figures 7 (c)-(d) and 8 (b) [60, 58]. As shown in figures 7 (e)-(f) and 8 (c),
a scenario which is a mixture of the two above can also arise: for example, at
small wave numbers the system shows stable oscillations, while pattern
formation appears at finite values of $k>0$. In general, the system parameters
can be tuned in such a way that the local negative maximum
Re$(\Lambda_{2}(k_{c}))$ becomes positive and global, thus leading to
instability. This type of instability is a conserved version of the well-known
Turing or type I instability. Contrarily to the standard Turing instability,
in our model, it can be attained also for $D^{\rm eff}_{1}=D^{\rm eff}_{2}$.
Similarly to the type II instability, it can be stationary or oscillatory.

Figure 7: We represent different instances of the onset of Turing instability,
upon variation of the $\mu_{2}$ parameter. Whenever an instability arises it
is characterized by the real part of the most unstable mode Re$(\Lambda_{2})$
becoming positive. In panel (a), we consider the case where $D_{1}^{\rm
eff}=D_{2}^{\rm eff}$: the eigenvalues are both real and instability may occur
only because of the mode associated with $\Lambda_{2}$. The parameters are
given by: $D_{c}=D_{1}=D_{2}=\bar{\rho}_{1}=\bar{\rho}_{2}=v_{1}=v_{2}=1$,
$D_{r,1}=D_{r,2}=5$, $\alpha_{1}=1.5$, $\alpha_{2}=0.1$, $\beta_{1}=1$,
$\beta_{2}=-0.5$, $\mu_{1}=10$, $\Omega_{1}=5$, $\Omega_{2}=-10$. In (b), the
two effective diffusivities are different and the eigenvalues are still real.
The parameters in panel (b) differing from those in (a) are: $D_{1}=2$,
$D_{2}=1.5$, $D_{r,1}=1$, $D_{r,2}=2$, $v_{2}=1.5$. In panels (c) and (d) we
represent respectively the real and imaginary parts of $\Lambda_{2}$ in the
case where for all values of $k$ the eigenvalues are complex conjugate,
corresponding to an oscillatory Turing instability. The parameters differing
from (b) are: $\alpha_{1}=4$, $\alpha_{2}=0.5$, $\beta_{1}=0.5$,
$\beta_{2}=1$, $\mu_{1}=-10$, $\Omega_{1}=10$, $\Omega_{2}=5$. Similarly, (e)
and (f) correspond to the case where the system presents stable complex
eigenvalues at small values of $k$, while the instability is associated with a
stationary pattern. The parameters from (c) and (d) are $\alpha_{1}=2.5$,
$\alpha_{2}=0.1$, $\beta_{1}=0.4$, $\Omega_{1}=5$, $\Omega_{2}=10$.  Figure 8:
Phase diagram representation of the Turing instability. These three figures
show how, by keeping all system parameters fixed, $k^{2}$Re$(\Lambda_{2})$
describes a state curve as $k$ is varied in the phase diagram. Similarly to
figure 5, the magenta area corresponds to the stable region of the homogeneous
disordered state, while blue and orange refer respectively to complex and real
unstable eigenvalues. The parametric curves correspond in panels in (a), (b),
and (c) respectively to the unstable (green) ones in figures 7 (b), 7 (c)-(d),
and 7 (e)-(f). All three curves start at $k=0$ on the stable side of the phase
diagram, and they end up in the origin, stable point of the phase diagram, as
$\bar{\rho}_{a}M_{a}\Pi_{a}\rightarrow 0$ as $k\rightarrow\infty$, while they
cross to the instability region at intermediate values of $k$. This is the
benchmark of the Turing instability. As in figure 6, the phase diagram changes
very little as $k$ is changed, and for simplicity, we have set $k/\kappa=10$.
Figure 9: In (a) and (c) we display the phase diagram of the effective
interactions: the orange area represents the stable homogeneous phase, which
is divided by the unstable phase by the red dashed line. The nature of the
unstable non-oscillatory regime is captured by the relative sign $s_{M}$ of
the effective mobilities $\bar{\rho}_{a}M_{a}$: the magenta region corresponds
to the separation of the two species ($s_{\rho}(k)<0$), the green one to
aggregation ($s_{\rho}(k)>0$). In (b) and (d), we show the emerging effective
interactions corresponding to (a) and (c) respectively. Specifically, in (b)
the effective interaction between the two type particles for $s_{M}>0$ and for
$s_{M}<0$ in (d). The red arrows indicate the direction of the reciprocal
interaction between particles of the same species, on the other hand, the
black ones refer to nonreciprocal interactions between particles of different
species.

###  IV.4 Effective interaction

To grasp the mechanisms that may lead to these types of instabilities, in the
next section we describe the type of nonreciprocal interactions that may arise
between the two species of Janus colloids. We are now interested in looking at
which of the two system modes, in the unstable regime, dominates the linear
instability and how it affects the growth of the relative concentration of the
two species $\delta\rho_{1}(\bm{k},t)/\delta\rho_{2}(\bm{k},t)$. In
particular, its sign
$s_{\rho}(\bm{k},t)=\operatorname{sign}\left(\delta\rho_{1}(\bm{k},t)/\delta\rho_{2}(\bm{k},t)\right)$
gives information on the type of instability: if negative, there will be a
local depletion of one species in favor of the other one, meaning separation,
otherwise local growth of concentration leads to aggregation.

We show in appendix D that $s_{\rho}(\bm{k},t)$ can be factorized as
$s_{\rho}(\bm{k},t)=s_{M}\,s_{m}(\bm{k},t)$, and it is given by

$s_{M}={\rm sign}(M_{1}/M_{2}),\;\;\;\text{and}\;\;\;s_{m}={\rm
sign}\left(\frac{D_{1}^{\rm eff}+\bar{\rho}_{1}M_{1}\Pi_{1}-D_{2}^{\rm
eff}-\bar{\rho}_{2}M_{2}\Pi_{2}-2\sqrt{\Delta}}{M_{1}\Pi_{1}}\right).$ (48)

It can be easily checked that in the unstable regime of the phase diagram, for
$\bar{\rho}_{2}M_{2}\Pi_{2}(k)>0$ the factor $s_{m}(k)<0$ and one has that the
sign of the instability is opposite that of the ratio of the effective
mobilities $s_{\rho}=-s_{M}$, while for $\bar{\rho}_{2}M_{2}\Pi_{2}(k)<0$ the
type of instability is reversed, i.e., $s_{\rho}=s_{M}$. Therefore,
aggregation of the two-particle species, or equivalently $s_{\rho}>0$, is
expected for $\bar{\rho}_{2}M_{2}\Pi_{2}(k)<0$ whenever $s_{M}>0$ or for
$\bar{\rho}_{2}M_{2}\Pi_{2}(k)>0$ if $s_{M}<0$ (green area in figure 9);
particles separate otherwise (magenta area in figure 9). The type of effective
interaction between the two species depends on the sign of the mobilities
$\bar{\rho}_{a}M_{a}$ and the production rate $\Pi_{a}(k)$. In particular, if
$\bar{\rho}_{a}M_{a}>0$ the particle of species $a$ will move towards regions
where the density of the substrate $c(\bm{r},t)$ decreases, which means that
the particle is attracted by consumers of the substrate with production rate
$\Pi_{b}(k)<0$ and repelled by producers with $\Pi_{b}(k)>0$. The details of
all the possible interactions between the two species are reported in figure
9. Because of the assumption $D_{1}^{\rm eff}>D_{2}^{\rm eff}$, the dynamics
of the first species are faster than that of the second species, and it
responds faster to the presence of chemical gradients, leading to the
prevailing of effective interaction felt by the first species with respect to
the second one. Referring to figure 9, this explains why independently of the
type of effective interaction experienced by the second species, if the first
species is attracted ( respectively repelled) by the second one the system
displays aggregation (respectively separation). In the unstable oscillatory
regime (blue area in figure 9) there is an alternation of depletion and
aggregation in time.

The evolution of these two density perturbations can be characterized by
looking at their phase and amplitude as
$\delta\rho_{1,2}(\bm{k},t)e^{-k^{2}\operatorname{tr}{\mathcal{G}}/2}=A_{1,2}\cos(\sqrt{|\Delta|}k^{2}t+\varphi_{1,2})$.
Albeit the two perturbations are destined to grow exponentially, in this
linear approximation the relative amplitude of the oscillations $A_{1}/A_{2}$
is constant in time, and it depends on the initial value of the perturbation.
On the other hand, the phase difference
$\Delta\varphi=\varphi_{1}-\varphi_{2}$ is independent of the initial
perturbation, and it is given by

$\tan\Delta\varphi=-\frac{2\sqrt{|\Delta|}}{\mathcal{G}_{11}-\mathcal{G}_{22}}=-\sqrt{-1-\frac{4\bar{\rho}_{1}M_{1}\Pi_{1}\bar{\rho}_{2}M_{2}\Pi_{2}}{(\delta
D+\bar{\rho}_{1}M_{1}\Pi_{1}-\bar{\rho}_{2}M_{2}\Pi_{2})^{2}}}.$ (49)

The phase difference, as shown in figure 10, takes value in $(-\pi/2,0)$: in
the unstable region ($\operatorname{tr}\mathcal{G}>0$) the phase difference is
always above $-\pi/2$, while it is in quadrature only in the stable phase for
$\operatorname{tr}\mathcal{G}=0$, compatibly with
$\bar{\rho}_{1}M_{1}\Pi_{1}\in(-\delta D,0)$ and
$\bar{\rho}_{2}M_{2}\Pi_{2}\in(0,\delta D)$.

Figure 10: We showcase the phase difference in equation (49) between the two
species in the oscillating state. (a) displays the phase difference in the
region of complex eigenvalues: on the left of the dashed line, the eigenvalues
are unstable, stable otherwise. (b) represents $\Delta\varphi$ for different
values of $\bar{\rho}_{2}M_{2}\Pi_{2}$ as a function of
$\bar{\rho}_{1}M_{1}\Pi_{1}$. The dashed vertical line separates values of
$\Delta\varphi$ corresponding to unstable eigenvalues (on the left) and to
stable ones (on the right).

##  V Conclusions

In this work, we have shown how effective nonreciprocal interactions arise in
a collection of two types of Janus colloids coupled to the same chemical
substrate. Including two species represents the minimal requirement for
nonreciprocal couplings between number density fields. We refer to figure 11
for a synoptic sketch of our main results. First, we have introduced the
corresponding single-particle dynamics, describing a set of Janus colloids
that can move and re-orient along the gradient of a chemical substrate, which
is produced or consumed by the colloids themselves. From this microscopic
description, we have derived the hydrodynamic equations for the relevant slow
modes, i.e., the particle density and the polarization field for each species,
which capture the collective behavior of the system. We have derived the
corresponding equations for small deviations from the spatially homogenous and
orientationally disordered state, allowing us to establish the linear
stability of this phase. Janus particles with chemical field-mediated
effective interactions are analogous to screened Coulombic systems. The
eigenvalues determining the linear stability of the system assume simple forms
when the translational and rotational diffusion coefficients and self-
propelling velocities are equal. Eigenvalue analysis predicts phase separation
as the most robust behavior at the largest lengthscales with a diffusion
coefficient whose sign is controlled both by the chemotactic drift and the
angular rotation. At the scale of the screening length, phoretic coupling
between number and orientation fields leads to oscillations. Oscillations
appear either through the mechanism of Hopf bifurcation or when the system
crosses an exceptional point. Two pairs of complex eigenvalues appear in the
most generic case, a scenario where two density fields and two orientation
fields undergo oscillations.

For large rotational diffusivity, the polar fields simply follow the density
gradients. In this scenario, the polarization degrees of freedom can be
expressed in terms of the density by an adiabatic approximation. The resulting
equations for the two density fields are linearized around the homogeneous
solution yield expressions for the interaction coefficients featuring wave-
vector-dependent activity and mobility coefficients. The interplay among the
associated effects leads to a nonreciprocal interaction between particles of
the same or different species leading to aggregation or separation phenomenon.
We have found different types of short-wavelength, stationary, or oscillatory
instabilities [61, 31]. In particular, contrary to what happens in standard
Turing instability, we get such instability even in the case of equal
effective diffusivities, as a consequence of the wave-vector-dependence of the
phoretic parameters. The analysis provides a direct link to the nonreciprocal
Cahn-Hilliard model (NRCH) [21, 22], and suggests that nonreciprocity should
be incorporated in the surface tension to approach a more complete theoretical
framework for scalar active densities. Recent papers have elucidated that NRCH
serves as a minimal model for known systems such as active-passive mixtures,
mass-conserving reaction-diffusion systems, and active gels [62, 61]. To the
best of our knowledge, our paper is the first work that starts from the
microscopic model of a chemically active swimmer including self-propulsion and
orientational dynamics to enumerate the various contributions to effective
intra-species and inter-species interactions (both reciprocal and
nonreciprocal) in terms of single-particle phoretic or enzymatic activity [63,
64] and mobilities thus providing several routes to realizing the NRCH. In
general, the speed of self-propulsion and chemical activity could depend on
the number densities thus providing a route observing the effect of
nonlinearity in nonreciprocal interactions [25].

An explicit manifestation of non-reciprocity is to enforce reactive couplings
between thermodynamic fluxes that should not be so coupled. Gradients of
chemical potentials should be coupled dissipatively with symmetric
coefficients [65], while velocity fields and density are advectively coupled
through coefficients of the same magnitude [66, 67, 68, 69]. Our system
presents two scenarios when Onsager’s principle is violated in both forms –
cross-couplings between densities of different species and coupling between
the longitudinal component of the polarization with density.

Taking cues from the analysis, it is important to explore the full dynamical
behavior of the system through a solution of the equations presented here or
in agent-based simulations of the microscopic model. A condensation of the
longitudinal component of the polarity only, a state called asters in [45]
deserves a thorough study examining questions such as long-range correlations
[70, 71, 72]. We expect a proliferation of defects for a single species,
somewhat similar to and yet distinct (i.e. occurring through a different
mechanism) from the defects observed in a Malthusian flock [73]. For two
species, where the stability analysis shows all modes to be unstable, we
speculate a state with interacting defects. Our work can be generalized in
several ways – to multicomponent mixtures interacting with several substrates
[20], coupling the mixture to a momentum-conserving fluid [74, 75], dynamics
at an interfaces and close to boundaries [76, 77, 78], and entropy production
[79]. Several aspects of our work can be generalized to other versions of
tactic systems - whether it is phototaxis [80] or quorum sensing [81].
Finally, our work illustrates that the physics of active mixtures represents a
rich area of research and presents many predictions that can be tested in
experiments [82, 38, 83].

Figure 11: Synoptic diagram of the main results of the paper. We schematize
our work’s logical flow and results. We have started from two sets of coupled
nonlinear differential equations, each describing the time evolution of the
density and polarization field for a given species of Janus colloids. Then, we
have discussed their linearization around a homogeneous disordered state. We
have split the analysis into two different parameters’ regimes: i) we have
first considered equal diffusivities and self-propelling velocities for the
two species, ii) then the limit of fast relaxation of the polarization field.

##  Appendix A Coarse-graining of the Langevin dynamics

In this section, we show how to derive equations (5) and
$\eqref{eq:ptruncated}$. We start by considering the evolution of
$\mathcal{P}_{a}(\bm{r},\bm{n},t)$, the probability density of a particle of
species $a\in\\{1,\dots,n\\}$ to be at position $\bm{r}$ and orientation
$\bm{n}$ at time $t$, defined as

$\mathcal{P}_{a}(\bm{r},\bm{n},t)=\Bigg{\langle}\sum_{i=1}^{N_{a}}\delta(\bm{r}-\bm{r}_{a,i})\delta(\bm{n}-\bm{n}_{a,i})\Bigg{\rangle},$
(50)

for a set of $N_{a}$ Janus colloids of species $a$, where the average is taken
with respect to the noise contribution to the dynamics. Then, considering that
the process in equation (1) follows Stratonovich convention of stochastic
calculus, it can be easily shown that $\mathcal{P}_{a}(\bm{r},\bm{n},t)$
satisfies the following Fokker-Planck equation

$\displaystyle\partial_{t}\mathcal{P}_{a}(\bm{r},\bm{n},t)=$
$\displaystyle-\nabla\cdot\left\\{\left[v_{a}\bm{n}-\mu_{a}\nabla
c\right]\mathcal{P}_{a}(\bm{r},\bm{n},t)\right\\}+D_{a}\nabla^{2}\mathcal{P}_{a}(\bm{r},\bm{n},t)$
(51) $\displaystyle-\mathcal{R}\cdot\left[\Omega_{a}\left(\bm{n}\times\nabla
c\right)\mathcal{P}_{a}(\bm{r},\bm{n},t)\right]+D_{r,a}\mathcal{R}^{2}\,\mathcal{P}_{a}(\bm{r},\bm{n},t),$

where $\mathcal{R}\equiv\bm{n}\times\nabla_{n}$ is the orientational gradient
operator. The first line on the right-hand side of equation (51) describes the
contribution to the probability flux due to drift and diffusion of the
particle position, whereas the second line to alignment interaction and
diffusion of its orientation.

In order to find an equation for $\rho_{a}$ and $\bm{p}_{a}$, we restrict our
analysis to the first and second moments of the orientation $\bm{n}$ by
closing the corresponding hierarchy of infinite many equations for the moments
generated from equation (51). We start by integrating equation (51) with
respect to $\bm{n}$, which leads to the time evolution of density of the
particles
$\rho_{a}(\bm{r},t)\equiv\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\mathcal{P}_{a}(\bm{r},\bm{n},t)$
of species $a$ in equation (5). Similarly, one can calculate the dynamics of
the polarization field
$\bm{p}_{a}(\bm{r},t)\equiv\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\bm{n}\,\mathcal{P}_{a}(\bm{r},\bm{n},t)$
as

$\partial_{t}\bm{p}_{a}=-\nabla\cdot\left[v_{a}\mathbb{Q}_{a}-\mu_{a}\bm{p}_{a}\nabla
c\right]-\frac{v_{a}}{3}\nabla\rho_{a}+\Omega_{a}\left[\frac{2}{3}\rho_{a}\nabla
c-\nabla c\cdot\mathbb{Q}_{a}\right]+(D_{a}\nabla^{2}-2D_{r,a})\bm{p}_{a},$
(52)

where we have introduced the nematic tensor
$\mathbb{Q}_{a}(\bm{r},t)=\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\left(\bm{n}\bm{n}-\mathbb{I}/3\right)\mathcal{P}_{a}(\bm{r},\bm{n},t)$.

In the evaluation of equation (52), which is obtained by multiplying by
$\bm{n}$ equation (51) we have calculated the following non trivial integrals:

  * i)

The first contribution that we consider is the $i$-component of

$\displaystyle-\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\,n_{i}\,\mathcal{R}\cdot\left[(\bm{n}\times\nabla
c)\,\mathcal{P}_{a}(\bm{r},\bm{n},t)\right]$
$\displaystyle=\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\,(\mathcal{R}_{j}n_{i})\,(\bm{n}\times\nabla
c)_{j}\,\mathcal{P}_{a}(\bm{r},\bm{n},t)$ (53)
$\displaystyle=-(\partial_{m}c)\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\,n_{k}n_{l}(\delta_{il}\delta_{km}-\delta_{im}\delta_{kl})\,\mathcal{P}_{a}(\bm{r},\bm{n},t)$
$\displaystyle=-(\partial_{j}c)\,\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,(n_{i}n_{j}-\delta_{ij})\mathcal{P}_{a}(\bm{r},\bm{n},t)$
$\displaystyle=\left[\frac{2}{3}\rho_{a}(\bm{r},t)\nabla c-\nabla
c\cdot\mathbb{Q}_{a}(\bm{r},t)\right]_{i}.$

On the right hand side of the first line we exploit the fact that
$\mathcal{R}_{i}$ satisfies the typical properties of gradient operator (that
is, chain differentiation rule and hence integration by part), and the
relation $\mathcal{R}_{i}n_{j}=-\epsilon_{ijk}n_{k}$. The second line is
obtained by contraction of the Levi Civita symbol
$\epsilon_{jik}\epsilon_{jlm}=\delta_{il}\delta_{km}-\delta_{im}\delta_{kl}$.
The last equality follows from the definition of the particle density
$\rho_{a}$ and the 3-dimensional nematic tensor $\mathbb{Q}_{a}$.

  * ii)

The second non-trivial term contributing to the $\bm{p}_{a}$ dynamics is given
by the one associated with angular diffusion, given by

$\displaystyle\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\,n_{i}\,\mathcal{R}^{2}\,\mathcal{P}_{a}(\bm{r},\bm{n},t)$
$\displaystyle=\int_{|\bm{n}|=1}\mathrm{d}\bm{n}\,\,\mathcal{P}_{a}(\bm{r},\bm{n},t)\,\mathcal{R}^{2}n_{i}=-2p_{a,i},$
(54)

that simply follows from the relation $\mathcal{R}^{2}n_{i}=-2n_{i}$.

Closure in the moment expansion can be attained by considering the case where
the nematic order is negligible $\mathbb{Q}_{a}=0$. This truncation of the
hierarchy of $\bm{n}$-moments in (51) simplifies equation (52) for $\bm{p}$ to
(6).

###  A.1  Microscopic origin of phoretic interactions

The phoretic couplings $\mu_{a}$, $\Omega_{a}$, and the production rates
$\alpha_{a}$, and $\beta_{a}$ can be expressed in terms of microscopic
parameters describing the geometric distribution of the mobility $\mu^{(a)}$
and activity $\alpha^{(a)}$ on the surface of spherical Janus particles. If we
restrict to a mobility $\mu^{(a)}$ and activity $\alpha^{(a)}$ that is axis-
symmetric with respect to $\bm{n}_{a}$, they are parametrized along the
surface of the Janus particle only via $\cos\theta$, where $\theta$ denotes
the angle with respect to the symmetry axis. It is then convenient to expand
$\mu^{(a)}(\cos\theta)$ and $\alpha^{(a)}(\cos\theta)$ in Legendre polynomials
according to

$\mu^{(a)}(\cos\theta)=\sum_{m=0}^{\infty}P_{m}(\cos\theta)\mu^{(a)}_{m},\quad\quad\alpha^{(a)}(\cos\theta)=\sum_{m=0}^{\infty}P_{m}(\cos\theta)\alpha^{(a)}_{m},$
(55)

where $P_{m}(\cos\theta)$ is the $m-$th degree Legendre polynomials, and the
coefficients $\mu^{(a)}_{m}$ and $\alpha^{(a)}_{m}$ are given by

$\mu^{(a)}_{m}=\left(m+\frac{1}{2}\right)\int_{0}^{\pi}\mathrm{d}\theta\,\sin\theta\,P_{m}(\cos\theta)\mu^{(a)}(\cos\theta),\quad\quad\alpha^{(a)}_{m}=\left(m+\frac{1}{2}\right)\int_{0}^{\pi}\mathrm{d}\theta\,\sin\theta\,P_{m}(\cos\theta)\alpha^{(a)}(\cos\theta).$
(56)

It can be shown that the velocity $\bm{v}_{a}$ and angular velocity
$\bm{\omega}_{a}$ [44, 50] due to phoretic interaction with the chemical
substrate are given by

$\bm{v}_{a}(\bm{r},t)=-\mu^{(a)}_{0}\nabla
c(\bm{r},t)+\frac{3}{10}\mu^{(a)}_{2}\left(\bm{n}_{a}\bm{n}_{a}-\frac{\mathbb{I}}{3}\right)\cdot\nabla
c(\bm{r},t),\quad\quad\bm{\omega}_{a}=-\frac{3\mu^{(a)}_{1}}{4R_{a}}\bm{n}_{a}\times\nabla
c(\bm{r},t),$ (57)

where $R_{a}$ is the radius of the colloid. From last equation we immediately
read the coefficients $\mu_{a}=\mu^{(a)}_{0}$ and
$\Omega_{a}=-3\mu^{(a)}_{1}/(4R_{a})$ appearing in equation (1), while
$\mu^{(a)}_{2}=0$ if we assume hemispherically coated Janus colloids.

###  A.2 Equation for the chemical substrate

The set of dynamical equations is completed by the evolution of the substrate
density field $c(\bm{r},t)$, which is given by

$\partial_{t}c-D_{c}(\nabla^{2}-\kappa^{2})c=\Pi(\bm{r},t)=\sum_{a=1}^{n}\left[\alpha_{a}\rho_{a}-\beta_{a}\nabla\cdot\bm{p}_{a}+O(R_{a}^{5})\right],$
(58)

where $R_{a}$ is the radius of the $a-$species Janus particle, and $\Pi$ the
local production rate of chemical substrate.

The average production rate $\Pi(\bm{r},t)$ can be expressed by integrating
the contribution coming from the local production rate
$\alpha^{(a)}(\cos\theta_{i_{a}})$, parametrized by $\theta_{i_{a}}$, the
angle between a point on the surface of each $i_{a}\in\\{1,\dots,N_{a}\\}$
Janus colloid and its axis, i.e.,

$\displaystyle\Pi(\bm{r},t)$
$\displaystyle=\Bigg{\langle}\sum_{a=1}^{n}\sum_{i_{a}=1}^{N_{a}}\int_{|\bm{R}_{i_{a}}|=R_{a}}\delta(\bm{r}-\bm{r}_{i_{a}}-\bm{R}_{i_{a}})\alpha^{(a)}(\cos\theta_{i_{a}})\Bigg{\rangle}$
(59)
$\displaystyle=\Bigg{\langle}\sum_{a=1}^{n}R_{a}^{2}\sum_{i_{a}=1}^{N_{a}}\sum_{m=0}^{\infty}\alpha^{(a)}_{m}\int_{0}^{\pi}\mathrm{d}\theta_{i_{a}}\,\sin\theta_{i_{a}}\int_{0}^{2\pi}\mathrm{d}\varphi_{i_{a}}\,P_{m}(\cos\theta_{i_{a}})\left[1-R_{a}\hat{\bm{r}}_{i_{a}}\cdot\nabla+O(R_{a}^{2})\right]\delta(\bm{r}-\bm{r}_{i_{a}})\Bigg{\rangle}$
$\displaystyle=\sum_{a=1}^{n}4\pi
R_{a}^{2}\Bigg{\langle}\sum_{i_{a}=1}^{N_{a}}\sum_{m=0}^{\infty}\alpha^{(a)}_{m}\left\\{\delta_{m,0}-\frac{R_{a}}{3}\delta_{m,1}(\bm{n}_{i_{a}}\cdot\nabla)+O(R_{a}^{2})\right\\}\delta(\bm{r}-\bm{r}_{i_{a}})\Bigg{\rangle}$
$\displaystyle=\sum_{a=1}^{n}4\pi
R_{a}^{2}\left[\alpha^{(a)}_{0}\,\rho_{a}(\bm{r},t)-R_{a}\frac{\alpha^{(a)}_{1}}{3}\,\nabla\cdot\bm{p}_{a}(\bm{r},t)+O(R_{a}^{2})\right],$

where $\bm{r}$ is a point on the surface of the particle, $\bm{r}_{i_{a}}$ the
location of the center of the $i_{a}$-th particle of species $a$, and
$\bm{R}_{i_{a}}=\bm{r}-\bm{r}_{i_{a}}$ is a point on the particle in the
particle’s reference frame. In the second line, we expand the activity
$\alpha^{(a)}$ in Legendre polynomials according to equation (55), we Taylor-
expand the Dirac delta for $R_{a}\ll r$, and we make the surface integral
explicit. The last two lines follow from integration and the definition of
$\rho_{a}$ and $\bm{p}_{a}$. From the above equation, we read

$\alpha_{a}=4\pi R_{a}^{2}\alpha^{(a)}_{0},\quad\quad\beta_{a}=\frac{4}{3}\pi
R_{a}^{3}\alpha^{(a)}_{1}.$ (60)

##  Appendix B Construction of the phase diagram

We now discuss the steps that we follow to construct the stability phase
diagram presented in the main text. We consider a two-dimensional linear with
the associated two eigenvalues

$\displaystyle\Lambda_{1,2}=\frac{\operatorname{tr}\mathcal{G}}{2}\mp\sqrt{\Delta},\;\;\;\text{with}\;\;\;\Delta=\left(\frac{\operatorname{tr}\mathcal{G}}{2}\right)^{2}-\operatorname{det}\mathcal{G}.$
(61)

We now summarise the results of the eigenvalue analysis in a stability
diagram. Generically, we can distinguish three significant regions denoted by
$R_{\Delta,1,2}$ delimited by curves $C_{\Delta,1,2}$:

$\begin{split}&R_{\Delta}:\Delta\leq 0,\\\
&R_{1}:\operatorname{tr}\mathcal{G}(k)=\Lambda_{1}+\Lambda_{2}<0,\\\
&R_{2}:\operatorname{det}\mathcal{G}(k)=\Lambda_{1}\Lambda_{2}>0,\end{split}\quad\quad\begin{split}&C_{\Delta}:\Delta=0,\\\
&C_{1}:\operatorname{tr}\mathcal{G}=0,\\\
&C_{2}:\operatorname{det}\mathcal{G}=0.\end{split}$ (62)

$R_{\Delta,1,2}$ identify the regions with real eigenvalue, positive trace of
$\mathcal{G}$ and positive determinant $\rm{det}\mathcal{G}$ respectively. The
intersection of $R_{1}$ and $R_{2}$ identifies the stability region of the
system. The region $R_{\Delta}$ corresponds to the region of complex
eigenvalues. The boundary $C_{\Delta}$ is the line of exceptional points,
where the eigenvalues are equal. Note that $R_{\Delta}$ is always contained in
$R_{2}$, implying that by crossing the portion of $C_{1}$ belonging to
$R_{\Delta}$ the stability of the system changes. To illustrate one
application of the considerations above, as an example, we construct in detail
the phase diagram shown in figure 5.

We can easily determine the topology of these regions in the plane of the
parameters $(\bar{\rho}_{1}M_{1}\Pi_{1},\bar{\rho}_{2}M_{2}\Pi_{2})$ while
keeping $D_{2}^{\rm eff}$ and $\delta D>0$ fixed. It can be simply checked
from equation (38) for $\Delta$, that this curve is a parabola defined for
$\bar{\rho}_{1}M_{1}\Pi_{1}\leq 0$ with symmetry axis
$\bar{\rho}_{2}M_{2}\Pi_{2}+\bar{\rho}_{1}M_{1}\Pi_{1}=0$ and vertex $(-\delta
D/4,\delta D/4)$. It is represented in figure 5 by the blue dashed curve: the
interior of the parabola corresponds to complex eigenvalues, while its
complementary $R_{\Delta}$ to the region of real eigenvalues.

Along the curve $C_{1}$ defined by the line
$\bar{\rho}_{2}M_{2}\Pi_{2}+\bar{\rho}_{1}M_{1}\Pi_{1}+D_{1}^{\rm
eff}+D_{2}^{\rm eff}=0$, the eigenvalues
$\Lambda_{1,2}=\pm\sqrt{\bar{\rho}_{1}M_{1}\Pi_{1}\bar{\rho}_{2}M_{2}\Pi_{2}}$
are equal and opposite. A complex conjugate pair of $\Lambda_{1,2}$ change
their sign on crossing $C_{1}$. Thus $C_{1}$ lying in $R_{\Delta}$ represents
points where Hopf bifurcation occurs dividing it into two regions - one where
oscillations grow and the other where they decay. A pair of real $\Lambda_{1}$
are equal and opposite on the curve $C_{1}$, which means it lies in the region
where one of the eigenmodes is unstable. Moreover, $C_{1}$ is parallel to the
axis of symmetry of $C_{\Delta}$, such that they intersect only at one point
$P_{\star}=\left(-(D_{1}^{\rm eff})^{2}/\delta D,(D_{2}^{\rm eff})^{2}/\delta
D\right)$ where the eigenvalues are both null (blue star in figure 5). In
figure 5 this former branch of $C_{1}$ is represented by the pale green semi-
line that originates from $P_{\star}$.

The expression of the curve $C_{2}$ is given by line
$\bar{\rho}_{1}M_{1}\Pi_{1}/D_{1}^{\rm
eff}+\bar{\rho}_{2}M_{2}\Pi_{2}/D_{2}^{\rm eff}+1=0$ and its corresponding
eigenvalues read $\Lambda_{1}=\operatorname{tr}\mathcal{G}$ and
$\Lambda_{2}=0$. Note that $C_{2}$ is tangent to $C_{\Delta}$ at $P_{\star}$
where also $\Lambda_{1}=0$. Thus, as in the case of $C_{1}$, the point
$P_{\star}$ splits $C_{2}$ into two semi-lines with different behavior: one
for $\bar{\rho}_{1}M_{1}\Pi_{1}>(D_{1}^{\rm eff})^{2}/\delta D$ where
$\Lambda_{1}$ is stable (red semi-line in figure 5) and the other that lays in
the (real) instability region. Indeed, for
$\bar{\rho}_{1}M_{1}\Pi_{1}>(D_{1}^{\rm eff})^{2}/\delta D$ above $C_{2}$ we
have real and stable eigenvalues, that become unstable below.

Similar considerations allow us to obtain the phase diagrams in figure 2.

##  Appendix C  Details of single species dynamics

Already at its linear description in equation (9), it is a very hard task to
deal with the complexity of the multi-species dynamics. For this reason, it is
useful to get some insights from the single component case. Indeed, as we have
anticipated in section III, it allows us to qualitatively understand certain
simplified regimes of the multi-species case. For the one species case, the
linearized dynamics of a perturbation to the homogeneous solution
$\rho(\bm{r},t)=\bar{\rho}$, $\bm{p}(\bm{r},t)=0$, and
$c(\bm{r},t)=\bar{c}=\alpha\bar{\rho}/(D_{c}\kappa^{2})$ is then given by

$\displaystyle\partial_{t}\delta\rho(\bm{k},t)=-k^{2}\left[\mu\bar{\rho}\,\delta
c(\bm{k},t)+D\,\delta\rho(\bm{k},t)\right]-ikv\,\delta
p_{\parallel}(\bm{k},t),$ (63) $\displaystyle\partial_{t}\delta
p_{\parallel}(\bm{k},t)=-\frac{ik}{3}\left[v\,\delta\rho(\bm{k},t)-2\Omega\bar{\rho}\,\delta
c(\bm{k},t)\right]-(2D_{r}+Dk^{2})\,\delta p_{\parallel}(\bm{k},t),$
$\displaystyle\delta
c(\bm{k},t)=\frac{\alpha\delta\rho(\bm{k},t)-ik\beta\delta
p_{\parallel}(\bm{k},t)}{D_{c}\left(\kappa^{2}+k^{2}\right)},\quad\partial_{t}\delta\bm{p}_{\perp}(\bm{k},t)=-(2D_{r}+Dk^{2})\delta\bm{p}_{\perp}(\bm{k},t).$

Note that, to ensure the physical requirement $\bar{c}>0$ of positive
substrate density, this description entails a positive $\alpha>0$ and
$\kappa\neq 0$: the substrate is on average created by the Janus colloids, and
its dynamics has to be screened. Equation (63) tells us that perturbation to
the homogeneous density profile $\delta\rho(\bm{k},t)$ is influenced only by
the component along $\bm{k}$ of the polar field, while its orthogonal
contribution $\delta\bm{p}_{\perp}(\bm{k},t)$ is exponentially suppressed over
time due to angular diffusion. For this reason, we restrict our analysis to
the evolution of $\delta\rho$ and $\delta p_{\parallel}$.

Accordingly, we can express the linearized dynamics of the system as

$\partial_{t}(\delta\rho(\bm{k},t),\delta
p_{\parallel}(\bm{k},t))^{T}=\mathcal{G}(\bm{k})(\delta\rho(\bm{k},t),\delta
p_{\parallel}(\bm{k},t))^{T},$ (64)

where we identify the one-species dynamical matrix as

$\mathcal{G}(\bm{k})=-\left(\begin{array}[]{cc}Dk^{2}+K(k)\frac{\bar{\rho}\mu\alpha}{D_{c}}&ik\left(v-K(k)\frac{\bar{\rho}\mu\beta}{D_{c}}\right)\\\\[7.0pt]
ik\left(\frac{v}{3}-K(k)\frac{2\bar{\rho}\Omega\alpha}{3k^{2}D_{c}}\right)&2D_{r}+Dk^{2}-K(k)\frac{2\bar{\rho}\Omega\beta}{3D_{c}}\\\
\end{array}\right),$ (65)

with $K(k)\equiv k^{2}/(k^{2}+\kappa^{2})$. We recall that, at the single
particle level, the phoretic drift $\mu$ and alignment interaction $\Omega$
determine how the Janus particles respond to the gradient distribution of the
substrate $\nabla c$: for positive values of $\mu$ the particles escape from
high $c$-concentration regions, while they point towards high concentration
region for $\mu<0$; for $\Omega>0$ particles align along $\nabla c$ (high
concentration) and anti-align otherwise. For $\beta>0$ more substrate
particles are produced in the orientation $\bm{n}$ direction, they are
consumed in the opposite case.

Figure 12: In panel (a), we show all the possible channels of effective
interactions in $\mathcal{G}$ due to phoretic effects. On the top left we
represent the effect of phoresis on the effective diffusivity of the
linearized dynamics corresponding to $\mathcal{G}_{11}$. The first column
identifies the type of perturbation to the homogenous phase: in this case a
local collection of Janus colloids with random orientation $\bm{n}$ (black
arrow) producing the chemical substrate in a homogeneous way (yellow area) in
a concentration higher than $\bar{\rho}$. In the second and third columns, the
effect of the chemotactic drift $\mu\nabla c$ is shown, whose orientation is
represented by the blue arrows. Similarly, in the top right corner, we account
for the advective channel due to chemotactic drift and asymmetric production
chemical described by the parameter $\beta$ (green area). On the bottom left,
we display the effect of $\alpha$ and the chemotactic alignment (red arrows)
to the pressure channel. Finally, the bottom right panel sketches the effect
of the effective interactions on rotational diffusion. In panel (b), we show
the regions of the phase diagram - similar to the one in figure 2 \- where the
elements of the dynamical matrix $\mathcal{G}_{ij}$ lead to a change of sign
of the diffusivity and velocity positive terms at the single phoretic channel
level, that is, when they are interpreted as the response coefficient of the
$i$-field to the solely $j$-field perturbation. Here we have defined the
auxiliary parameters $\sigma=\alpha/\beta$ and
$\Delta_{\sigma}(k)=2D_{r}-v[\sigma+k^{2}/(3\sigma)]$, that together with the
combination of parameters on the axes of the phase diagram, uniquely
characterize the system.

Therefore, we can interpret the $\mathcal{G}_{ij}$ element of the dynamical
matrix $\mathcal{G}$ as the response coefficient of the $i$-field to a
$j$-field small perturbation, where $i=1$ and $i=2$ identify respectively the
density and polar field fluctuations. Note that the phoretic contribution to
$\mathcal{G}$ appears only via one of the two production rates $\alpha$ and
$\beta$, and only one of the chemotactic interactions $\mu$ and $\Omega$. This
is a consequence of the linearized dynamics and it allows us to consider one
type of phoretic effect at the time to the $i$-field, depending on the nature
of the $j$-field perturbations to the homogeneous phase. To better understand
the effect of these response coefficients we refer to figure 12, where we
display the effect of these (linearized) effective interactions on the
stability of the homogeneous phase:

  * •

The first element $\mathcal{G}_{11}$ measures how small density fluctuations
are amplified or suppressed while considering no perturbation to its
disordered orientational component. This channel of interaction is associated
with the effect of the homogeneous production of chemicals, parametrized by
$\alpha$, and the chemotactic drift $\mu$. This diffusion channel (shown on
the top left of figure 12(a)) stabilizes the homogeneous phase for $\mu>0$
because of chemorepulsive interactions, whereas it can lead to aggregation in
the case of effectively chemoattractive ones for $\mu<0$. In figure 12 (b) we
mark in yellow the region of parameters where $\mathcal{G}_{11}$ becomes
negative: phoretic interaction may lead to a negative effective diffusivity.

  * •

The second channel, described by $\mathcal{G}_{12}$, tells us how a small
perturbation to the orientational disorder can affect the density in terms of
$\beta\mu$. To fix ideas we set $\beta>0$ (top right of figure 12 (a)),
corresponding to the asymmetric production of chemicals along the self-
propelling direction $\bm{n}$. If we consider as a perturbation of the polar
field the local alignment of a set of Janus colloids, we get that the effect
of chemorepulsion tends to reduce the effective self-propelling velocity,
whereas it increases for $\mu<0$. We refer to this channel of interaction as
the advective one, since here phoretic interactions contribute to the
advective part of the density dynamics. The red region in figure 12 (b)
corresponds to $\mathcal{G}_{12}/(ik)$ negative, i.e. when the effective self-
propelling velocity becomes negative.

  * •

The third channel of interaction, associated to $\mathcal{G}_{21}$ (bottom
left in figure 12 (a)), can be interpreted as a pressure term, in analogy with
the case of Toner Tu type of equations [72, 70, 71]. In particular, self-
propulsion acts as a mechanism to restore orientational disorder due to
possible local perturbation of the density field. The phoretic contribution to
this mechanism is represented by $\Omega\alpha$, that is to the homogenous
production of chemicals and the alignment interaction. If $\Omega>0$ Janus
particles point towards high concentration regions of chemical leading to
larger values of $\delta p_{\parallel}$ and possibly to aggregation of
particles. On the other hand, if $\Omega<0$ particles point away from high
concentration regions, thus destroying local order. Also in this case phoretic
interaction may lead to a negative value of the net velocity
$\mathcal{G}_{21}/(ik)$ (green area in figure 12 (b)). Indeed, in the limit of
fast relaxation of the polar field, e.g., $D_{r}+Dk^{2}\gg
2K(k)\Omega\beta/(3D_{c})$ ($\mathcal{G}_{22}<0$), we get $\delta
p_{\parallel}(\bm{k},t)\simeq-(\mathcal{G}_{21}/\mathcal{G}_{22})\delta\rho(\bm{k},t)$
and
$\partial_{t}\delta\rho(\bm{k},t)=(\operatorname{det}{\mathcal{G}}/\mathcal{G}_{22})\delta\rho(\bm{k},t)$.
If (for simplicity) we set $\beta=0$, it is then apparent that the stability
of the effective diffusivity
$\operatorname{det}{\mathcal{G}}/\mathcal{G}_{22}=\mathcal{G}_{11}-v[vk^{2}/3-K(k)2\Omega\alpha/(3D_{c})]/(2D_{r}+Dk^{2})$
depends crucially on the sign and intensity of $\mathcal{G}_{21}$.

  * •

Finally, the term $\mathcal{G}_{22}$ identifies the contribution of alignment
interaction to the effective orientational diffusivity. Namely, it carries
information about how the polar field is affected by local orientational
order. As for $\mathcal{G}_{12}$ we first set $\beta>0$ (bottom right of
figure 12 (a)): in the case of locally ordered particles, for $\Omega>0$ the
particles are most likely be already aligned to high concentration regions
thus creating aggregation, alignment, and instability in the homogeneous
phase; for $\Omega<0$ the particles start to rotate since chemical is produced
in correspondence of $\bm{n}$ while the same $\bm{n}$ tries to escape high
concentration regions, thus stabilizing disordered configurations. In general,
as $\beta\Omega>0$ phoretic interactions favor local orientational order,
while they stabilize the homogeneous phase for $\beta\Omega<0$. The blue
region in figure 12 (b) denotes the region of parameters where the effective
orientational diffusivity becomes negative because of phoretic interactions.

It is worth noting that all the phoretic contributions in $\mathcal{G}$ are
proportional to $K(k)\bar{\rho}/D_{c}$: their contributions become stronger
for a higher concentration of particles and increasing values of substrate
screening length $1/\kappa$. The region of parameters corresponding to the
instability of each of these channels is highlighted in figure 12 (b). Then,
the full phase diagram (the light blue shaded curves in figure 12 (b)) is not
just given by the superposition of the four individual response coefficients
$\mathcal{G}_{ij}$, but from their interplay. This mixture of different
channels is captured by the eigenvalues of the dynamical matrix
$\mathcal{G}(k)$. They can be expressed via the trace
$\operatorname{tr}{\mathcal{G}}$ and determinant
$\operatorname{det}{\mathcal{G}}$ of $\mathcal{G}$ as

$\displaystyle\Lambda_{1/2}(k)$
$\displaystyle=\frac{\operatorname{tr}\mathcal{G}(k)\mp\sqrt{\left(\operatorname{tr}\mathcal{G}(k)\right)^{2}-4\operatorname{det}\mathcal{G}(k)}}{2}$
(66)
$\displaystyle=-\left[D_{r}+k^{2}D+\frac{K(k)}{2}\frac{\bar{\rho}}{D_{c}}\left(\mu\alpha-\frac{2}{3}\Omega\beta\right)\pm\sqrt{\Delta(k)}\right],$

with discriminant $\Delta(k)$ explicitly given by

$\displaystyle\Delta(k)=$
$\displaystyle\left[D_{r}+k^{2}D+\frac{K(k)}{2}\frac{\bar{\rho}}{D_{c}}\left(\mu\alpha-\frac{2}{3}\Omega\beta\right)\right]^{2}$
(67)
$\displaystyle-\left(k^{2}D+K(k)\frac{\bar{\rho}\mu\alpha}{D_{c}}\right)\left(2D_{r}+Dk^{2}-K(k)\frac{2\bar{\rho}\Omega\beta}{3D_{c}}\right)-\left(v-K(k)\frac{\bar{\rho}\mu\beta}{D_{c}}\right)\left(\frac{k^{2}v}{3}-K(k)\frac{2\bar{\rho}\Omega\alpha}{3D_{c}}\right).$

Note that the presence of the conserved field density $\delta\rho(\bm{k},t)$
implies that one of the eigenvalues must vanish at $k=0$. This property is
satisfied by $\Lambda_{2}(k)$ and reveals an important feature of the
unscreened ($\kappa=0$) case, i.e., in the limit $k\rightarrow 0$ and
$\kappa\rightarrow 0$ do not commute. Indeed, if one takes the limit
$\kappa\rightarrow 0$ before setting $k=0$, $K(k)=1$ becomes constant,
inconsistently with $\Lambda_{2}(0)=0$. The first eigenvalue $\Lambda_{1}(k)$
is associated with a non-conserved mode that at $k=0$ reduces to
$\Lambda_{1}(0)=-2D_{r}$, reflecting the fact that angular diffusion
stabilizes the homogeneous disordered phase. At small wave numbers and
screened interactions $\kappa>0$, the eigenvalues behave as

$\displaystyle\Lambda_{1}(k)$
$\displaystyle=-2D_{r}-k^{2}\left[\frac{\bar{\rho}\Omega}{3D_{c}\kappa^{2}}\left(-2\beta+v\alpha\right)-\frac{v^{2}}{6D_{r}}+D\right]+O(k^{4}),$
(68) $\displaystyle\Lambda_{2}(k)$
$\displaystyle=-k^{2}\left[\frac{\bar{\rho}\alpha}{D_{c}\kappa^{2}}\left(\mu-\frac{v\Omega}{3D_{r}}\right)+D_{\rm
eff}(0)\right]+O(k^{4}).$

It as apparent from equation (68) that for small values of $k$ the first
eigenvalue $\Lambda_{1}(k)$ is stable due to the angular diffusion, while the
eigenvalue $\Lambda_{2}(k)$ associated with the conserved field is stable for

$\frac{\bar{\rho}\alpha}{D_{c}}\left(\mu-\frac{v\Omega}{3D_{r}}\right)+\kappa^{2}D_{\rm
eff}(0)\geq 0.$ (69)

The last relation states that, at the macroscopic scale, the homogeneous phase
is stable if, even in the case of effective attractive interaction among the
particles, i.e., $\alpha(\mu-v\Omega/(3D_{r}))<0$, the traslational and
orientational noise prevails. For fully unscreened interactions
$k\gg\kappa\simeq 0$ and $K(k)\simeq 1$, at leading order in $\kappa/k$ the
eigenvalues become

$\Lambda_{1,2}=-\left[D_{r}+\frac{\bar{\rho}}{2D_{c}}\left(\mu\alpha-\frac{2}{3}\Omega\beta\right)\pm\sqrt{\left[D_{r}+\frac{\bar{\rho}}{2D_{c}}\left(\mu\alpha-\frac{2}{3}\Omega\beta\right)\right]^{2}-\frac{\bar{\rho}\mu\alpha}{D_{c}}\left(2D_{r}-\frac{2\bar{\rho}\Omega\beta}{3D_{c}}\right)-\left(v-\frac{\bar{\rho}\mu\beta}{D_{c}}\right)\frac{2\bar{\rho}\Omega\alpha}{3D_{c}}}\right],$
(70)

which implies that, in the case of long-range interactions among the Janus
colloids, at a large enough scale the system presents two non-conserved modes
that can be both stabilized by strong enough orientational disorder.
Physically, being the number of particles a conserved quantity, in the limit
$k\rightarrow 0$ there must always be a conserved mode that vanishes as
$k^{2}$, signaling that the unscreened regime holds only for $\kappa\ll k$ and
$k$ very small. Notably, the very same behavior has been found in the two
species case in section III.1 for equal diffusivity and self-propelling
velocity between the two species.

##  Appendix D  Effective interaction

Here we complement the calculations missing in section IV.4. We are interested
in solving the linear system in equation (27) with initial conditions
$\delta\rho(\bm{k},0)$ in the unstable regime. In particular, we need an
expression for the ratio $\delta\rho_{1}(\bm{k},t)/\delta\rho_{2}(\bm{k},t)$
in the large $t$ limit. This is captured by the ratio of the elements of the
eigenvector corresponding to the eigenvalue $\Lambda_{2}$, associated with the
most unstable mode which dominates the dynamics. The eigenvectors are given by

$\hat{e}_{1,2}(k)=\left(\begin{array}[]{c}\frac{D_{1}^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}-D_{2}^{\rm
eff}-\bar{\rho}_{2}M_{2}\Pi_{2}}{2}\pm\sqrt{\Delta}\\\\[5.0pt]
\bar{\rho}_{2}M_{2}\Pi_{1}\\\ \end{array}\right).$ (71)

In the non-oscillatory regime, the relative amplitude of the two density
perturbations is given by

$\displaystyle\frac{\delta\rho_{1}(\bm{k},t)}{\delta\rho_{2}(\bm{k},t)}$
$\displaystyle=\frac{\bar{\rho}_{1}M_{1}}{\bar{\rho}_{2}M_{2}}\frac{D_{1}^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}-D_{2}^{\rm
eff}-\bar{\rho}_{2}M_{2}\Pi_{2}-2\sqrt{\Delta}}{2\bar{\rho}_{1}M_{1}\Pi_{1}},$
(72)

at leading order in $e^{k^{2}\Lambda_{2}t}$.

In the unstable oscillatory regime (blue area in figure 9) there is an
alternation of depletion and aggregation in time according to

$\displaystyle\frac{\delta\rho_{1}(\bm{k},t)}{\delta\rho_{2}(\bm{k},t)}=$ (73)
$\displaystyle\frac{\sin\left(\sqrt{|\Delta|}\,k^{2}t\right)\left[\delta\rho_{1}(\bm{k},0)(-D_{1}^{\rm
eff}+D_{2}^{\rm
eff}-\bar{\rho}_{1}M_{1}\Pi_{1}+\bar{\rho}_{2}M_{2}\Pi_{2})-\delta\rho_{2}(\bm{k},0)2\bar{\rho}_{1}M_{1}\Pi_{2}\right]+2\sqrt{|\Delta|}\,\delta\rho_{1}(\bm{k},0)\cos\left(\sqrt{|\Delta|}\,k^{2}t\right)}{\sin\left(\sqrt{|\Delta|}\,k^{2}t\right)\left[\delta\rho_{2}(\bm{k},0)(D_{1}^{\rm
eff}-D_{2}^{\rm
eff}+\bar{\rho}_{1}M_{1}\Pi_{1}-\bar{\rho}_{2}M_{2}\Pi_{2})-\delta\rho_{1}(\bm{k},0)2\bar{\rho}_{2}M_{2}\Pi_{1}\right]+2\sqrt{|\Delta|}\,\delta\rho_{2}(\bm{k},0)\cos\left(\sqrt{|\Delta|}\,k^{2}t\right)},$

where the common exponentially growing factor
$e^{k^{2}\operatorname{tr}{\mathcal{G}}/2}$ cancel in the ratio. To better
characterize this (linearly) oscillating phase we look at the rescaled
variables
$\delta\rho_{1,2}(\bm{k},t)e^{-k^{2}\operatorname{tr}{\mathcal{G}}/2}$, whose
time evolution describes an ellipse. The evolution of these two periodic
trajectories can characterized by looking at their phase and amplitude as
$\delta\rho_{1,2}(\bm{k},t)e^{-k^{2}\operatorname{tr}{\mathcal{G}}/2}=A_{1,2}\cos(\sqrt{|\Delta|}k^{2}t+\varphi_{1,2})$,
where

$\displaystyle\tan\varphi_{1}$
$\displaystyle=\frac{\delta\rho_{1}(\bm{k},0)(\mathcal{G}_{11}-\mathcal{G}_{22})+\delta\rho_{2}(\bm{k},0)2\mathcal{G}_{12}}{2\sqrt{|\Delta|}\,\delta\rho_{1}(\bm{k},0)},$
(74) $\displaystyle\tan\varphi_{2}$
$\displaystyle=-\frac{\delta\rho_{2}(\bm{k},0)(\mathcal{G}_{11}-\mathcal{G}_{22})-\delta\rho_{1}(\bm{k},0)2\mathcal{G}_{21}}{2\sqrt{|\Delta|}\,\delta\rho_{2}(\bm{k},0)},$
$\displaystyle A_{1}$
$\displaystyle=2\sqrt{\mathcal{G}_{12}\left[(\delta\rho_{1}(\bm{k},0))^{2}\mathcal{G}_{21}+\delta\rho_{1}(\bm{k},0)\delta\rho_{2}(\bm{k},0)(\mathcal{G}_{11}-\mathcal{G}_{22})+(\delta\rho_{2}(\bm{k},0))^{2}\mathcal{G}_{12}\right]},$
$\displaystyle A_{2}$
$\displaystyle=2\sqrt{\mathcal{G}_{21}\left[(\delta\rho_{1}(\bm{k},0))^{2}\mathcal{G}_{21}-\delta\rho_{1}(\bm{k},0)\delta\rho_{2}(\bm{k},0)(\mathcal{G}_{11}-\mathcal{G}_{22})+(\delta\rho_{2}(\bm{k},0))^{2}\mathcal{G}_{12}\right]},$

where we recall that complex eigenvalues exist only for
$\bar{\rho}_{1}M_{1}\Pi_{1}<0$ and $\bar{\rho}_{2}M_{2}\Pi_{2}>0$.

## References

  * 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, J. Phys.: Condens. Matter 32, 193001 (2020).
  * Bowick  _et al._ [2022] M. J. Bowick, N. Fakhri, M. C. Marchetti, and S. Ramaswamy, Phys. Rev. X 12, 010501 (2022).
  * Golestanian [2012] R. Golestanian, Phys. Rev. Lett. 108, 038303 (2012).
  * Golestanian [2022a] R. Golestanian, in  _Active Matter and Nonequilibrium Statistical Physics: Lecture Notes of the Les Houches Summer School: Volume 112, September 2018_ (Oxford University Press, 2022).
  * Helbing and Molnár [1995] D. Helbing and P. Molnár, Phys. Rev. E 51, 4282 (1995).
  * Uchida and Golestanian [2010] N. Uchida and R. Golestanian, Phys. Rev. Lett. 104, 178103 (2010).
  * Hickey  _et al._ [2023] D. J. Hickey, R. Golestanian, and A. Vilfan, PNAS 120, e2307279120 (2023).
  * Drescher  _et al._ [2009] K. Drescher, K. C. Leptos, I. Tuval, T. Ishikawa, T. J. Pedley, and R. E. Goldstein, Phys. Rev. Lett. 102, 168101 (2009).
  * Ziepke  _et al._ [2022] A. Ziepke, I. Maryshev, I. S. Aranson, and E. Frey, Nat. Commun. 13, 6727 (2022).
  * Loos and Klapp [2020] S. A. M. Loos and S. H. L. Klapp, New J. Phys. 22, 123051 (2020).
  * Osat and Golestanian [2023] S. Osat and R. Golestanian, Nat. Nanotechnol. 18, 79 (2023).
  * Dinelli  _et al._ [2023] A. Dinelli, J. O’Byrne, A. Curatolo, Y. Zhao, P. Sollich, and J. Tailleur, Nat. Commun. 14, 7035 (2023).
  * Soto and Golestanian [2014] R. Soto and R. Golestanian, Phys. Rev. Lett. 112, 068301 (2014).
  * Soto and Golestanian [2015] R. Soto and R. Golestanian, Phys. Rev. E 91, 052304 (2015).
  * Nasouri and Golestanian [2020] B. Nasouri and R. Golestanian, Phys. Rev. Lett. 124, 168003 (2020).
  * Cocconi  _et al._ [2023] L. Cocconi, H. Alston, and T. Bertrand, Phys. Rev. Res. 5, 043032 (2023).
  * Agudo-Canalejo and Golestanian [2019] J. Agudo-Canalejo and R. Golestanian, Phys. Rev. Lett. 123, 018101 (2019).
  * Ouazan-Reboul  _et al._ [2021] V. Ouazan-Reboul, J. Agudo-Canalejo, and R. Golestanian, Eur. Phys. J. E 44, 1 (2021).
  * Ouazan-Reboul  _et al._ [2023a] V. Ouazan-Reboul, J. Agudo-Canalejo, and R. Golestanian, Nat. Commun. 14, 10.1038/s41467-023-40241-w (2023a).
  * Ouazan-Reboul  _et al._ [2023b] V. Ouazan-Reboul, R. Golestanian, and J. Agudo-Canalejo, Phys. Rev. Lett. 131, 128301 (2023b).
  * Saha  _et al._ [2020] S. Saha, J. Agudo-Canalejo, and R. Golestanian, Phys. Rev. X 10, 041009 (2020).
  * You  _et al._ [2020] Z. You, A. Baskaran, and M. C. Marchetti, PNAS 117, 19767 (2020).
  * Frohoff-Hülsmann  _et al._ [2021a] T. Frohoff-Hülsmann, J. Wrembel, and U. Thiele, Phys. Rev. E 103, 042602 (2021a).
  * Frohoff-Hülsmann and Thiele [2023a] T. Frohoff-Hülsmann and U. Thiele, Phys. Rev. Lett. 131, 107201 (2023a).
  * Saha and Golestanian [2022] S. Saha and R. Golestanian, Effervescent waves in a binary mixture with non-reciprocal couplings (2022), arXiv:2208.14985 [cond-mat.stat-mech] .
  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Nature 592, 363 (2021).
  * Knežević  _et al._ [2022] M. Knežević, T. Welker, and H. Stark, Sci. Rep. 12, 19437 (2022).
  * Dadhichi  _et al._ [2020] L. P. Dadhichi, J. Kethapelli, R. Chajwa, S. Ramaswamy, and A. Maitra, Phys. Rev. E 101, 052601 (2020).
  * Loos  _et al._ [2023] S. A. M. Loos, S. H. L. Klapp, and T. Martynec, Phys. Rev. Lett. 130, 198301 (2023).
  * Cavagna  _et al._ [2017] A. Cavagna, I. Giardina, A. Jelic, S. Melillo, L. Parisi, E. Silvestri, and M. Viale, Phys. Rev. Lett. 118, 138003 (2017).
  * Duan  _et al._ [2023a] Y. Duan, J. Agudo-Canalejo, R. Golestanian, and B. Mahault, Phys. Rev. Lett. 131, 148301 (2023a).
  * Sinha and Chaudhuri [2024] A. Sinha and D. Chaudhuri, Soft Matter 20, 788–795 (2024).
  * Scheibner  _et al._ [2020] C. Scheibner, A. Souslov, D. Banerjee, P. Surówka, W. T. M. Irvine, and V. Vitelli, Nat. Phys. 16, 475 (2020).
  * Poncet and Bartolo [2022] A. Poncet and D. Bartolo, Phys. Rev. Lett. 128, 048002 (2022).
  * Shankar  _et al._ [2022] S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, Nat. Rev. Phys. 4, 380 (2022).
  * Golestanian  _et al._ [2005] R. Golestanian, T. B. Liverpool, and A. Ajdari, Phys. Rev. Lett. 94, 220801 (2005).
  * Howse  _et al._ [2007] J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Phys. Rev. Lett. 99, 048102 (2007).
  * Sharan  _et al._ [2023] P. Sharan, A. Daddi-Moussa-Ider, J. Agudo-Canalejo, R. Golestanian, and J. Simmchen, Small 19, 2300817 (2023).
  * Lavergne  _et al._ [2019] F. A. Lavergne, H. Wendehenne, T. Bäuerle, and C. Bechinger, Science 364, 70 (2019).
  * Durve  _et al._ [2018] M. Durve, A. Saha, and A. Sayeed, Eur. Phys. J. E 41, 49 (2018).
  * Ivlev  _et al._ [2015] A. V. Ivlev, J. Bartnick, M. Heinen, C.-R. Du, V. Nosenko, and H. Löwen, Phys. Rev. X 5, 011035 (2015).
  * Anderson [1989] J. L. Anderson, Annu. Rev. Fluid Mech. 21, 61 (1989).
  * Rückner and Kapral [2007] G. Rückner and R. Kapral, Phys. Rev. Lett. 98, 150603 (2007).
  * Golestanian  _et al._ [2007] R. Golestanian, T. Liverpool, and A. Ajdari, New J. Phys. 9, 126 (2007).
  * Saha  _et al._ [2014] S. Saha, R. Golestanian, and S. Ramaswamy, Phys. Rev. E 89, 062316 (2014).
  * Golestanian [2022b] R. Golestanian, in  _Active matter and nonequilibrium statistical physics_ , Lecture Notes of the 2018 Les Houches Summer School (Oxford University Press, London, England, 2022).
  * Liebchen  _et al._ [2015] B. Liebchen, D. Marenduzzo, I. Pagonabarraga, and M. E. Cates, Phys. Rev. Lett. 115, 258301 (2015).
  * Liebchen  _et al._ [2017] B. Liebchen, D. Marenduzzo, and M. E. Cates, Phys. Rev. Lett. 118, 268001 (2017).
  * Saha  _et al._ [2019] S. Saha, S. Ramaswamy, and R. Golestanian, New J. Phys. 21, 063006 (2019).
  * Kanso and Michelin [2019] E. Kanso and S. Michelin, J. Chem. Phys. 150, 044902 (2019).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ (Springer, 2009).
  * Pavliotis [2014] G. A. Pavliotis,  _Stochastic Processes and Applications_ (Springer, 2014).
  * Romanczuk  _et al._ [2012] P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Eur. Phys. J. Spec. Top. 202, 1 (2012).
  * Stenhammar  _et al._ [2014] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E. Cates, Soft matter 10, 1489 (2014).
  * Siebert  _et al._ [2018] J. T. Siebert, F. Dittrich, F. Schmid, K. Binder, T. Speck, and P. Virnau, Phys. Rev. E 98, 030601 (2018).
  * Solon  _et al._ [2015] A. P. Solon, M. E. Cates, and J. Tailleur, Eur. Phys. J. Spec. Top. 224, 1231 (2015).
  * Cohen and Golestanian [2014] J. A. Cohen and R. Golestanian, Phys. Rev. Lett. 112, 068302 (2014).
  * Duan  _et al._ [2023b] Y. Duan, J. Agudo-Canalejo, R. Golestanian, and B. Mahault, Phys. Rev. Lett. 131, 148301 (2023b).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
  * Frohoff-Hülsmann  _et al._ [2021b] T. Frohoff-Hülsmann, J. Wrembel, and U. Thiele, Phys. Rev. E 103, 042602 (2021b).
  * Frohoff-Hülsmann and Thiele [2023b] T. Frohoff-Hülsmann and U. Thiele, Phys. Rev. Lett. 131, 107201 (2023b).
  * Brauns and Marchetti [2023] F. Brauns and M. C. Marchetti, Non-reciprocal pattern formation of conserved fields (2023), arXiv:2306.08868 [nlin.PS] .
  * Agudo-Canalejo  _et al._ [2018a] J. Agudo-Canalejo, P. Illien, and R. Golestanian, Nano Lett. 18, 2711 (2018a).
  * Agudo-Canalejo  _et al._ [2018b] J. Agudo-Canalejo, T. Adeleke-Larodo, P. Illien, and R. Golestanian, Acc. Chem. Res. 51, 2365 (2018b).
  * De Groot and Mazur [2013] S. R. De Groot and P. Mazur,  _Non-equilibrium thermodynamics_ (Courier Corporation, 2013).
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Prawar Dadhichi  _et al._ [2018] L. Prawar Dadhichi, A. Maitra, and S. Ramaswamy, J. Stat. Mech. Theory Exp. 2018, 123201 (2018).
  * Das  _et al._ [2002] J. Das, M. Rao, and S. Ramaswamy, EPL 60, 418 (2002).
  * Das  _et al._ [2004] J. Das, M. Rao, and S. Ramaswamy, Nonequilibrium steady states of the isotropic classical magnet (2004), arXiv:cond-mat/0404071 [cond-mat.soft] .
  * Toner and Tu [1998] J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998).
  * Toner and Tu [1995] J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995).
  * Toner  _et al._ [2005] J. Toner, Y. Tu, and S. Ramaswamy, Ann. Phys. 318, 170 (2005).
  * Besse  _et al._ [2022] M. Besse, H. Chaté, and A. Solon, Phys. Rev. Lett. 129, 268003 (2022).
  * Lushi  _et al._ [2018] E. Lushi, R. E. Goldstein, and M. J. Shelley, Phys. Rev. E 98, 052411 (2018).
  * Vinze and Michelin [2024] P. Vinze and S. Michelin, Phys. Rev. Fluids 9, 014202 (2024).
  * Malgaretti  _et al._ [2016] P. Malgaretti, M. N. Popescu, and S. Dietrich, Soft Matter 12, 4007 (2016).
  * Uspal  _et al._ [2015] W. E. Uspal, M. N. Popescu, S. Dietrich, and M. Tasinkevych, Soft Matter 11, 434 (2015).
  * Popescu  _et al._ [2017] M. N. Popescu, W. E. Uspal, and S. Dietrich, J. Phys. Condens. Matter 29, 134001 (2017).
  * Bebon  _et al._ [2024] R. Bebon, J. F. Robinson, and T. Speck, Thermodynamics of active matter: Tracking dissipation across scales (2024), arXiv:2401.02252 [cond-mat.soft] .
  * Lozano  _et al._ [2016] C. Lozano, B. ten Hagen, H. Löwen, and C. Bechinger, Nat. Commun. 7, 12828 (2016).
  * O’Byrne and Tailleur [2020] J. O’Byrne and J. Tailleur, Phys. Rev. Lett. 125, 208003 (2020).
  * Wang and Simmchen [2023] L. Wang and J. Simmchen, Soft Matter 19, 540 (2023).
  * Sánchez  _et al._ [2015] S. Sánchez, L. Soler, and J. Katuri, Angew. Chem. Int. Ed. 54, 1414 (2015).