# Non-reciprocal interactions enhance heterogeneity

Timoteo Carletti<EMAIL_ADDRESS>Department of Mathematics and
naXys, Namur Institute for Complex Systems, University of Namur, Rue Grafé 2,
B5000 Namur, Belgium  Riccardo Muolo  Department of Mathematics and naXys,
Namur Institute for Complex Systems, University of Namur, Rue Grafé 2, B5000
Namur, Belgium

###### Abstract

We study a process of pattern formation for a generic model of species
anchored to the nodes of a network where local reactions take place, and that
experience non-reciprocal non-local long-range interactions, encoded by the
network directed links. By assuming the system to exhibit a stable homogeneous
equilibrium whenever only local interactions are considered, we prove that
such equilibrium can turn unstable once suitable non-reciprocal non-local
long-range interactions are allowed for. Stated differently we propose
sufficient conditions allowing for patterns to emerge by using a non-symmetric
coupling, while initial perturbations about the homogenous equilibrium always
fade away by assuming reciprocal coupling, namely the latter is stable. The
instability, precursor of the emerging spatio-temporal patterns, can be traced
back, via a linear stability analysis, to the complex spectrum of an
interaction non-symmetric Laplace operator. The proposed theory is then
applied to several paradigmatic dynamical models largely used in the
literature to study the emergence of patterns or synchronisation. Taken
together, our results pave the way for the understanding of the many and
heterogeneous patterns of complexity found in ecological, chemical or physical
systems composed by interacting parts, once no diffusion takes place.

##  I Introduction

We live in an interconnected world Newman (2010); Barabási (2016) where
complex patterns Pastor-Satorras and Vespignani (2010) spontaneously emerge
from the intricate web of nonlinear interactions existing between the basic
units by which the system under study is made of Nicolis and Prigogine (1977).
These emergent structures can be found in the synchronised activity of
neurones, resulting from the exchange of electrochemical signals via synapses
Pikovsky  _et al._ (2001); Arenas  _et al._ (2008), as well as the geometric
visual hallucinations product of the retinocortical map linking the retina and
the striate cortex Bressloff  _et al._ (2001); Bressloff and D (2003). Such
self-organised structures can also manifest in groups of fireflies that flash
at unison, each one by observing the behaviour of their close neighbours
Strogatz (2003). Or they can materialise as the striped or spotted motifs on
the skin of the zebrafish Danio rerio, due to the long-range interactions
between melanophores and xantophores, without requiring diffusion nor any kind
of cell motion Nakamasu  _et al._ (2009); Bullara and De Decker (2015).

The latter phenomena, as well as many other ones modelled within a similar
framework of local reactions and long-range interactions without displacement
of the reacting species, cannot be ascribed to a Turing instability Turing
(1952), a widely used paradigm for pattern formation. Indeed the latter
requires a diffusive process, whereas the common key factor linking the above
examples is the immobility of the reacting species and the existence of a web
of long-range interactions due to some signal propagation.

In this work we set the theoretical basis for the understanding of self-
organisation in systems without diffusing species where asymmetric long-range
interactions play a pivotal role. Network science Newman (2010); Barabási
(2016) provides a natural framework where to study such phenomena. Indeed
local interactions can be described by a dynamical system evolving on each
node, that represents a portion of physical space, say a natural habitat or a
chemical reactor, large enough to contain enough species to describe their
evolution with a nonlinear ordinary differential equation. On the other hand,
non-local long-range interactions can be modelled by links connecting the
nodes of the network allowing thus species anchored to each distinct node to
communicate by exchanging some generic signals. Let us observe that our
definition of local and long-range dynamics can differ from the one sometimes
used in the literature, that considers local dynamics to be related to the in-
node processes but also to the interactions with species sitting in nodes at
distance one, i.e., first neighbourhood, while long-range interactions are
used to describe far away interacting nodes, i.e., distances larger than two.
For this reason in this work we name the latter process non-local long-range
interaction; let us notice however that when it will be clear from the
context, we will simply use long-range interactions. Let us finally stress
that the latter do not involve displacement of the reacting species, for this
reason the theory hereby developed applies to systems without diffusion. An
interesting example is the web of chemical light-triggered reactions obtained
by connecting reactors where the local concentration of chemicals, i.e., in
each node, determines the amount of light to be put onto connected distant
nodes, i.e., long-range interaction. Chemicals do not leave any nodes but they
interact at distance Gizynski and Gorecki (2017); Voorsluijs  _et al._ (2017);
Parrilla-Gutierrez  _et al._ (2020).

In the aforementioned examples, the system converges to an homogeneous
solution, being stationary or time varying, once only local interactions are
taken into account, i.e., long-range ones are silenced, while heterogeneous
solutions spontaneously emerge in presence of a suitable web of interactions
among the units. Let us observe that the existence of an attracting
homogeneous solution for the decoupled system is a natural and largely assumed
working hypothesis, see for instance Turing (1952) in the framework of Turing
instability or Pecora and Carroll (1998) for synchronisation.

A preliminary result in this direction has been proposed in Cencetti  _et al._
(2020) with the assumption of reciprocal non-local interactions and by
treating the latter in a mean-field setting. There, authors have been able to
prove that the stable homogeneous solution, existing once the long-range
interactions have been silenced, can turn unstable by introducing a suitable
symmetric non-local coupling and eventually lead to the emergence of a
spatially (and temporally) dependent solution. The aim of the present work is
to make one step forward by studying the general case of non-reciprocal
interactions and their impact on the system outcome. Indeed, the interactions
existing among the constituting units are often not symmetric; this is the
case of plant-animal mutualistic networks Bascompte  _et al._ (2003); Abramson
_et al._ (2011), the specific example of sheeps and deers Hollingdale  _et
al._ (2018) or the case of olfactory receptor neurones in the Drosophila
antenna Zhang  _et al._ (2019), just to mention a few.

Anticipating our conclusion, we claim that the diversity and heterogeneity of
patterns observed in Nature, being associated to spatial or temporal non-
homogeneous states, is enhanced by non-reciprocal long-range interactions. The
onset of the instability, precursor of the pattern, can be detected with a
linear stability analysis, providing a condition on the complex spectrum of a
non-symmetric consensus Laplacian operator resulting from the mean-field
ansatz as we will hereafter explain. The proposed framework is general enough
to cover systems of any dimension $d\geq 2$ displaying a fixed-point or a
limit cycle homogeneous solution, once we silence the long-range interactions,
by proving the existence of a suitable non-reciprocal web of interactions
driving the emergence of patterns while the latter can never exist in the case
of symmetric support. To emphasise the relevance of the proposed theory to
many and different research domains, we applied the developed theory to
several paradigmatic models largely used in the literature to study the
emergence of patterns or synchronisation. In conclusion, the proposed
mechanism provides the way for alternative routes to pattern formation, beyond
the Turing one Turing (1952); Nakao and Mikhailov (2010), and to the emergence
of desynchronised states Pecora and Carroll (1998), suitable for all phenomena
where diffusion is not the main driver for the heterogeneity of the complex
patterns seen in Nature, opening thus new possibilities for modelling
ecological, chemical and physical interacting systems, endowed with non-
reciprocal couplings and without diffusing species.

##  II The model

Let us consider a dynamical system composed by $n$ identical parts and assume
the $d$-dimensional vector,
$\mathbf{x}^{(i)}(t)=(x^{(i)}_{1}(t),\dots,x^{(i)}_{d}(t))^{\top}$, to
represent the state of the $i$-th copy, for $i=1,\dots,n$. The isolated
systems are described by an ordinary differential equation 111Let us observe
that the theory hereby developed could be straightforwardly adapted to the
case of discrete time, where, i.e., the ODE is replaced with a finite
differences system., resulting from the assumption of a “well stirred”
distribution of species inside each node

$\frac{d\mathbf{x}^{(i)}}{dt}=\mathbf{f}(\mathbf{x}^{(i)})\quad\forall
i=1,\dots,n\,,$ (1)

where $\mathbf{f}$ is a generic nonlinear function responsible for the local
interactions, i.e., depending on the species anchored to the same node. Let us
now allow each system to possess non-local interactions, i.e., the growth rate
of the system anchored at node $i$ is influenced by some nonlinear function of
the amount of species in distant nodes. Moreover we assume the latter to be
described by

$\frac{d\mathbf{x}^{(i)}}{dt}=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})\quad\forall
i=1,\dots,n\,,$ (2)

where $A_{ij}$ is the (possibly weighted) non-symmetric adjacency matrix
encoding the long-range interactions, i.e., $A_{ij}=1$ if and only if node $j$
influences node $i$. Let $k_{i}^{(in)}=\sum_{j}A_{ij}$ be the in-degree of
node $i$ and observe that we allow $A_{ii}=1$, thus the in-degree takes into
account also the possible self-loops. Finally,
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})$ is a nonlinear function that
describes the effect of the $j$–th system on the $i$-th one. Moreover we
require
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$,
namely the self-interaction is represented by the original nonlinear function
$\mathbf{f}$ describing the evolution of the isolated systems. This is a
natural assumption allowing to recover the initial system (1) composed by
isolated units, i.e., once we set $\mathbf{A}=\mathbf{I}_{n}$, being the
latter the $n\times n$ identify matrix 222Let us observe that this condition
has been already used in Cencetti  _et al._ (2020) and it is related to the
one adopted in Gambuzza  _et al._ (2021).. Let us observe that the right hand
side of (2) can be rewritten as the average of the interactions perceived by
node $i$,
$\langle\mathbf{F}(\mathbf{x}^{(i)},\cdot)\rangle=\sum_{j}A_{ij}\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})/k_{i}^{(in)}$,
hence describing the mean-field ansatz. Finally, let us stress that species
cannot move across nodes: the long-range interactions are thus mediated by
some generic signals, as we have proposed in the introduction and we will
discuss in the following. The model (2) describes thus local interactions
coupled with long-range ones without the need of diffusion, this is in
contrast with the fundamental assumption of Turing instability Turing (1952);
Nakao and Mikhailov (2010) where diffusion plays a key role, or in the
synchronisation based on diffusive-like couplings Pecora and Carroll (1998);
Pecora  _et al._ (1997) (see Fig. 1 for a schematic representation of the
proposed model).

Figure 1:  A schematic visual representation of the model (2). Nodes (large
white circles) contain two species (blue and red dots). We then focus on the
growth rate of the “blue” species in the $i$-th node, that results from two
terms; the first one is the local reaction $\mathbf{f}$, represented by the
light-blue and light-red curved arrows, involving only “blue” and “red”
species indexed by $i$ (green oval). The second contribution arises from the
long-range interactions $\mathbf{F}$, represented as violet triangles pointing
from each nearby node to the $i$-th one. The average of these terms impacts
the growth rate of the “blue” specie in node $i$ (orange oval).

By defining the matrix

$\mathcal{L}_{ij}=\frac{A_{ij}}{k^{(in)}_{i}}-\delta_{ij}\,,$ (3)

and by using the property of the function $\mathbf{F}$, we can rewrite Eq. (2)
as

$\frac{d\mathbf{x}^{(i)}}{dt}=\mathbf{f}(\mathbf{x}^{(i)})+\sum_{j}\mathcal{L}_{ij}\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})\quad\forall
i=1,\dots,n\,.$ (4)

Once the long-range interactions are silenced and each node interacts only
with itself, then $\mathcal{L}=0$ and thus Eq. (4) reduces to (1). Let us
observe that in the case of reciprocal interactions, the above defined matrix
$\mathcal{L}$ corresponds to the consensus Laplace operator Ghosh  _et al._
(2014); Lambiotte  _et al._ (2011); Olfati-Saber  _et al._ (2007); Krause
(2009), also named reactive Laplace matrix in Cencetti  _et al._ (2020), whose
spectrum is real and non-positive. In the case under scrutiny, involving non-
reciprocal interactions, the spectrum is generally complex but one can prove
that $\Lambda^{(1)}=0$ is still an eigenvalue associated to the uniform
eigenvector $\phi^{(1)}\sim(1,\dots,1)^{\top}$. Moreover the Gershgorin circle
theorem Geršgorin (1931) allows to prove that the real part of the spectrum of
$\mathcal{L}$ is contained in the strip $[-2,0]$ in the complex plane, hence
$\mathcal{L}$ is stable having all the eigenvalues but $0$ with a negative
real part. Observe that by using a symmetric network the eigenvalues are real
and non-positive, even if the consensus Laplace is not symmetric (see Appendix
A). Let us stress that the operator $\mathcal{L}$ results from the mean-field
assumption and not from any kind of diffusive-like process as in the cases
where the combinatorial Laplace operator arises because of the Fick law.

Let us assume $\mathbf{x}^{(i)}(t)=\mathbf{s}(t)$, $i=1,\dots,n$, to be a
solution of the initial system (1); then because of the above hypothesis on
$\mathbf{F}$ and of the definition of $k_{i}^{(in)}$, it is also a spatially,
i.e., node dependent solution of Eq. (4). To study the bifurcation of patchy
solutions from the stable homogeneous one, $\mathbf{x}^{(i)}(t)=\mathbf{s}(t)$
for all $i$, we consider a node dependent perturbation about the latter
solution, i.e., $\mathbf{u}^{(i)}(t)=\mathbf{x}^{(i)}(t)-\mathbf{s}(t)$, whose
evolution can be studied by inserting it into Eq. (4) and then keeping only
first order terms, assumed to be small enough. Hence we obtain for all
$i=1,\dots,n$

$\frac{d\mathbf{u}^{(i)}}{dt}=\left[\mathbf{J}_{1}(\mathbf{s}(t))+\mathbf{J}_{2}(\mathbf{s}(t))\right]\mathbf{u}^{(i)}+\sum_{j}\mathcal{L}_{ij}\mathbf{J}_{2}(\mathbf{s}(t))\mathbf{u}^{(j)}\,,$
(5)

where we have introduced the Jacobian matrices
$\mathbf{J}_{1}(\mathbf{s})=\partial_{\mathbf{x}_{1}}\mathbf{F}(\mathbf{x}_{1},\mathbf{x}_{2})\rvert_{(\mathbf{s},\mathbf{s})}$,
i.e., the derivatives are computed with respect to the first group of
variables, and
$\mathbf{J}_{2}(\mathbf{s})=\partial_{\mathbf{x}_{2}}\mathbf{F}(\mathbf{x}_{1},\mathbf{x}_{2})\rvert_{(\mathbf{s},\mathbf{s})}$,
i.e., the derivatives are performed with respect to the second group of
variables. In both cases the derivatives are evaluated on the reference
solution $\mathbf{s}(t)$.

This latter equation encodes $n$ linear systems involving matrices with size
$d\times d$. To progress with the analytical understanding, we assume the
existence of an eigenbasis 333Let us observe that this hypothesis can be
relaxed (see Appendix B) and one could obtain similar results by invoking the
Jordan block decomposition as done in Nishikawa and Motter (2006) in the case
of synchronisation. For sake of definitiveness we preferred to present our
results under this restrictive assumption, but allowing for straightforward
analysis. of right eigenvectors $\phi^{(\alpha)}$, $\alpha=1,\dots,n$ for
$\mathcal{L}$ and then, following the ideas pioneered by Turing (1952); Pecora
and Carroll (1998), we project the former equation onto each eigendirection,
i.e.,
$\mathbf{u}^{(i)}=\sum_{\alpha}\mathbf{u}^{(\alpha)}\phi^{(\alpha)}_{i}$, to
eventually obtain (see Appendix B)

$\frac{d\mathbf{u}^{(\alpha)}}{dt}=\left[\mathbf{J}(\mathbf{s}(t))+\mathbf{J}_{2}(\mathbf{s}(t))\Lambda^{(\alpha)}\right]\mathbf{u}^{(\alpha)}\quad\forall\alpha=1,\dots,n\,,$
(6)

where $\mathbf{J}=\mathbf{J}_{1}+\mathbf{J}_{2}$ and $\Lambda^{(\alpha)}$ is
the eigenvalue relative to the eigenvector $\phi^{(\alpha)}$. The above
equation enables us to infer the (in)stability of the homogeneous solution,
$\mathbf{s}(t)$, by studying the Master Stability Function Pecora and Carroll
(1998); Huang  _et al._ (2009), namely the largest Lyapunov exponent of Eq.
(6), considered as a function of $\Lambda^{(\alpha)}$.

To make one step further in the study of the problem, let us hypothesise that
each isolated system converges to the same stationary point, i.e., the stable
homogeneous solution is time independent, $\mathbf{s}(t)=\mathbf{s}_{0}$.
Hence Eq. (6) rewrites for all $\alpha=1,\dots,n$, as

$\frac{d\mathbf{u}^{(\alpha)}}{dt}=\left[\mathbf{J}(\mathbf{s}_{0})+\mathbf{J}_{2}(\mathbf{s}_{0})\Lambda^{(\alpha)}\right]\mathbf{u}^{(\alpha)}:=\mathbf{J}^{(\alpha)}\mathbf{u}^{(\alpha)}\,,$
(7)

Where the constant matrix $\mathbf{J}^{(\alpha)}$ has been defined by the
latter equation. The homogeneous solution will prove unstable to spatial
dependent perturbations if (at least) one eigenmode $\hat{\alpha}$ exists for
which the largest real part of the eigenvalues, $\lambda_{i}$, of
$\mathbf{J}^{(\hat{\alpha})}$ is positive, the latter being known in the
literature with the name of dispersion relation, hereby denoted by
$\rho_{\alpha}=\max_{i=1,\dots,d}\Re\lambda_{i}(\Lambda^{(\alpha)})$, where we
emphasised the dependence of the latter on the spectrum of the Laplace matrix.
Let us observe that the positivity of the dispersion relation initiates the
instability but it does not determine the final outcome of the system.
Convergence to other homogeneous solutions is possible if $\mathbf{f}$
possesses several zeros, namely if the system (1) exhibits multiple
equilibria. However, this kind of solution can arise with lower probability
with respect to a patchy one; indeed, assuming to act on a parameter in Eq.
(4), then a homogeneous zero of the latter generically bifurcates in an
heterogenous one.

##  III Results

For sake of pedagogy and to be able to determine closed and manageable
analytical formulas, let us assume the local systems to be $2$ dimensional,
i.e., $d=2$ in Eq. (1). Let us however emphasise that, as we will discuss
later on, our results go beyond this simplified framework and we can indeed
prove that non-reciprocal interactions can drive the emergence of heterogenous
patterns in any dimension $d\geq 2$, even when this is not possible using a
symmetric coupling. Assuming thus $d=2$, being the eigenvalues of
$\mathbf{J}^{(\alpha)}$ the solutions of the second order equation

$\lambda_{i}^{2}-\mathrm{tr}\mathbf{J}^{(\alpha)}\lambda_{i}+\det\mathbf{J}^{(\alpha)}=0\,,$
(8)

we can adapt to the present case the analysis done in Asllani  _et al._
(2014a, 2020) and express the condition for the onset of instability, i.e.,
$\rho_{\hat{\alpha}}>0$, as follows

$\exists\hat{\alpha}>1\text{ st
}(\Im\Lambda^{(\hat{\alpha})})^{2}S_{2}(\Re\Lambda^{(\hat{\alpha})})\leq-
S_{1}(\Re\Lambda^{(\hat{\alpha})})\,,$ (9)

where $S_{2}(\xi)$, resp. $S_{1}(\xi)$, is a second, resp. fourth, degree
polynomial in $\xi$ (see Appendix B for more details and the explicit form of
$S_{1}$ and $S_{2}$ in terms of the model parameters).

To illustrate the potential of the theory let us consider several examples of
dynamical systems largely used in the literature as paradigmatic models for
synchronisation or patterns emergence.

###  III.1 The Brusselator model

Let us consider the Brusselator model Prigogine and Nicolis (1967); Prigogine
and Lefever (1968), often invoked in the literature as a paradigm nonlinear
reaction scheme for studying self-organised phenomena, synchronisation Muolo
_et al._ (2021), Turing patterns Nakao and Mikhailov (2010); Asllani  _et al._
(2014b, 2016); Petit  _et al._ (2017); Asllani  _et al._ (2020) and
oscillation death Koseska  _et al._ (2013); Lucas  _et al._ (2018). The key
feature of the model is the presence of two species, reacting via a cubic
nonlinearity

$\begin{dcases}\frac{du}{dt}&=1-(b+1)u+cu^{2}v\\\ \frac{dv}{dt}&=bu-
cu^{2}v\,,\end{dcases}$ (10)

where $b>0$ and $c>0$ act as tunable model parameters. One can easily realise
the existence of a unique equilibrium $u^{*}=1$ and $v^{*}=b/c$, that results
stable if the Jacobian of the reaction part evaluated on the equilibrium
$\mathbf{J}_{\mathrm{Bxl}}=\left(\begin{smallmatrix}b-1&c\\\
-b&-c\end{smallmatrix}\right)$, has a negative trace,
$\mathrm{tr}\mathbf{J}_{\mathrm{Bxl}}=b-c<1$, and a positive determinant
$\det\mathbf{J}_{\mathrm{Bxl}}=c>0$.

The model can be cast in the framework presented in the previous section by
setting $\mathbf{x}^{(i)}=(u_{i},v_{i})$ and
$\mathbf{f}(\mathbf{x}^{(i)})=(1-(b+1)u_{i}+cu_{i}^{2}v_{i},bu_{i}-cu_{i}^{2}v_{i})$,
for $i=1,\dots,n$ to denote the $n$ isolated systems. For sake of concreteness
let us consider the coupling given by

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=(1-(b+1)u_{i}+cu_{i}^{2}v_{i},bu_{j}-cu_{i}^{2}v_{i})\,,$

where $\mathbf{x}^{(j)}=(u_{j},v_{j})$. Such function clearly satisfies the
constraint
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$,
hence Eq. (2) becomes (see Appendix C)

$\begin{dcases}\frac{du_{i}}{dt}&=1-(b+1)u_{i}+cu_{i}^{2}v_{i}\\\
\frac{dv_{i}}{dt}&=bu_{i}-cu_{i}^{2}v_{i}+b\sum_{j}\mathcal{L}_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,,$ (11)

namely a set of $n$ Brusselator models (10) coupled via long-range connections
in the second variable.

Given the above coupling we can explicitly compute the polynomials
$S_{1}(\xi)$ and $S_{2}(\xi)$ as a function of the model parameters (see
Appendix C) and characterise the instability region defined by (9) as reported
in panel a) of Fig. 2 where we show in the complex plane
$(\Re\Lambda,\Im\Lambda)$ the regions for which the instability condition is
satisfied (grey), for a given set of parameters. Patterns do emerge if there
exists at least one eigenvalue $\Lambda^{(\hat{\alpha})}$ belonging to this
region. For sake of simplicity, we hereby assume the non-reciprocal
interactions to be described by a directed Erdős-Rényi network made of $50$
nodes and the probability to create a directed link to be $0.05$. The
symmetric coupling is obtained by considering all the existing links to be
reciprocal ones and, as already noticed, the eigenvalues are negative real
numbers. In conclusion, if the model parameters shape an instability region
that does not intersect the real negative axis (see panel a) in Fig. 2), then
only an asymmetric coupling can drive the instability and the ensuing
(oscillatory) patterns (see panel b) in Fig. 2), while this is impossible for
any web of reciprocal long-range interactions (see panel c) in Fig. 2) and the
system solution converges toward the homogenous solution 444Throughout the
work, the numerical simulations have been performed by initialising the system
$\delta$-close to the homogeneous equilibrium by drawing uniformly random
perturbations in $(-\delta,\delta)$. Then its time evolution has been
numerically simulated using a $4$-th order Runge-Kutta method over a time span
of the order of $-\log\delta/\max_{\alpha}\rho_{\alpha}$, namely sufficiently
long to (possibly) increase the $\delta$-perturbation up to a macroscopic
size..

Figure 2: Instability region and patterns for the Brusselator model. In panel
a) we report the region of the complex plane $(\Re\Lambda,\Im\Lambda)$ for
which the instability condition is satisfied (grey). For the chosen parameters
values ($b=4.3$ and $c=5.0$) we can observe that the instability region does
not intersect the real axis and thus only non-reciprocal coupling can exhibit
complex eigenvalues (white dots) entering into the instability region and thus
initiate the pattern as shown in panel b) where we report $u_{i}(t)$ vs $t$
starting from initial conditions close to the stable equilibrium, $u_{*}=1$.
Any symmetric coupling determines real eigenvalue (black dots in panel a))
that cannot give rise to the instability, as shown in panel c) where we report
$u_{i}(t)$ vs $t$ starting from the same initial conditions used in panel b).
The underlying coupling is obtained with a directed Erdős-Rényi network with
$n=50$ nodes and a probability for a direct link to exist between two nodes is
$p=0.05$.

###  III.2 The Mimura-Murray model

The second model we consider is the Mimura-Murray system Mimura and Murray
(1978); Nakao and Mikhailov (2010); Challenger  _et al._ (2015); Asllani  _et
al._ (2015). It also involves two reacting species, that can be associated to
prey and predator, whose densities denoted by $u$ and $v$ evolve via the
nonlinear system of ODEs:

$\begin{dcases}\frac{du}{dt}&=u\left(\frac{a+bu-u^{2}}{c}-v\right)\\\
\frac{dv}{dt}&=v\left(u-(1+dv)\right)\,,\end{dcases}$ (12)

where $a,b,c$ and $d$ are positive parameters. The model possesses $6$
equilibria, whose stability and positivity depend on the value of the chosen
parameters. We here focus on the fixed point $(u_{*},v_{*})$

$u_{*}=1+\frac{bd-2d-c+\sqrt{\Delta}}{2d}\text{ and
}v_{*}=\frac{bd-2d-c+\sqrt{\Delta}}{2d^{2}}\text{ where
}\Delta=(bd-2d-c)^{2}+4d^{2}(a+b-1)\,,$ (13)

and assume $a=35$, $b=15$, $c=20$ and $d=0.4$ which in turn implies
$(u_{*},v_{*})\sim(2.28,3.20)$. The Jacobian matrix evaluated at the fixed
point reads
$\mathbf{J}_{\mathrm{MM}}\sim\left(\begin{smallmatrix}1.19&-2.28\\\
3.20&-1.28\end{smallmatrix}\right)$, hence, $\det\mathbf{J}_{\mathrm{MM}}>0$
and $\mathrm{tr}\mathbf{J}_{\mathrm{MM}}<0$ and the fixed point is a stable
equilibrium.

Let us consider the following long-range coupling

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=\left(u_{i}\left(\frac{a+bu_{j}-u_{i}^{2}}{c}-v_{j}\right),v_{i}\left(u_{j}-(1+dv_{i})\right)\right)\,,$

where $\mathbf{x}^{(i)}=(u_{i},v_{i})$ and $\mathbf{x}^{(j)}=(u_{j},v_{j})$.
Let us observe that the present coupling contains nonlinear terms in the
variables $u_{i}v_{j}$ and $u_{j}v_{i}$, while in the previous section we
dealt with a linear term. Let
$\mathbf{f}(\mathbf{x}^{(i)})=\left(u_{i}\left(\frac{a+bu_{i}-u_{i}^{2}}{c}-v_{i}\right),v_{i}\left(u_{i}-(1+dv_{i})\right)\right)$,
then
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$,
and we are thus in the framework of Eq. (2) that now can be rewritten as (see
Appendix D)

$\begin{dcases}\frac{du_{i}}{dt}&=u_{i}\left(\frac{a+bu_{i}-u_{i}^{2}}{c}-v_{i}\right)+\frac{b}{c}u_{i}\sum_{j}\mathcal{L}_{ij}u_{j}-u_{i}\sum_{j}\mathcal{L}_{ij}v_{j}\\\
\frac{dv_{i}}{dt}&=v_{i}\left(u_{i}-(1+dv_{i})\right)+v_{i}\sum_{j}\mathcal{L}_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,.$ (14)

One can realise that we are dealing with $n$ Mimura-Murray models (12) coupled
via nonlinear long-range connections in both variables.

As explained in the Appendix D one can numerically determine the instability
region by evaluating the polynomials $S_{1}(\xi)$ and $S_{2}(\xi)$ as a
function of the model parameters (see panel a) of Fig. 3). We can observe that
the instability region does not intersect the negative real axis and thus only
complex eigenvalues, associated to non-symmetric coupling can enter into such
region. This is the case reported in panel a) of Fig. 3 and confirmed by the
time evolution of the density $u_{i}(t)$ (panel b) in the case of non
symmetric long-range interactions and panel c) in the symmetric case). The
non-reciprocal interactions are obtained with a directed Erdős-Rényi network
made of $50$ nodes and the probability to create a directed link to be $0.05$.

Figure 3: Instability region and patterns for the Mimura-Murray model. In
panel a) we report the region of the complex plane $(\Re\Lambda,\Im\Lambda)$
for which the instability condition is satisfied (grey). For the chosen
parameters values ($a=35$, $b=15$, $c=20$ and $d=2/5$) we can observe that the
instability region does not intersect the real axis and thus only non-
reciprocal coupling can exhibit complex eigenvalues (white dots) entering into
the instability region and thus initiate the pattern as shown in panel b)
where we report $u_{i}(t)$ vs $t$ starting from initial conditions close to
the stable equilibrium, $u_{*}\sim 2.28$, $v_{*}\sim 3.20$. Any symmetric
coupling determines real eigenvalue (black dots in panel a)) that cannot give
rise to the instability, as shown in panel c) where we report $u_{i}(t)$ vs
$t$ starting from the same initial conditions used in panel b). The underlying
coupling is obtained with a directed Erdős-Rényi network with $n=50$ nodes and
a probability for a direct link to exist between two nodes is $p=0.05$.

###  III.3 The Volterra model

By remaining in the framework of ecological models, we consider an application
dealing with a Volterra model McKane and Newman (2005) modified by the
assumption of the presence of non-reciprocal interactions among animals. For a
sake of pedagogy, we hereby decided to present our results with a simplified
and abstract model that describes the interactions of abstract prey and
predators in an ecological setting, let us however observe that our theory
could be modified as to include other phenomena, e.g., harvesting Layek and
Pati (2019); Pati  _et al._ (2020). The model we are interesting in results
thus be described by:

$\begin{dcases}\frac{dx}{dt}=-dx+c_{1}xy\\\ \frac{dy}{dt}=ry-
sy^{2}-c_{2}xy\,.\end{dcases}$ (15)

Here $x$ denotes the concentration of predators, while $y$ stands for the
preys, and all the parameters are assumed to be positive. The Volterra model
(15) admits a nontrivial fixed-point, $x^{*}=\frac{c_{1}r-sd}{c_{1}c_{2}}$,
$y^{*}=\frac{d}{c_{1}}$, which is positive and stable, provided that
$c_{1}r-sd>0$. Condition hereby assumed to hold true.

Following the above presented scheme, let us now consider the long-range
coupling

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=\left(-dx_{i}+ac_{1}x_{j}y_{i}+(1-a)c_{1}x_{i}y_{j},ry_{i}-sy_{i}^{2}-c_{2}x_{j}y_{i}\right)\,,$

where $\mathbf{x}^{(i)}=(x_{i},y_{i})$ and $\mathbf{x}^{(j)}=(x_{j},y_{j})$.
The parameter $a$ belongs to $[0,1]$ and weights the contribution arising from
the nonlinear terms $x_{i}y_{j}$ and $x_{j}y_{i}$; such term disappears in the
local dynamics. Let
$\mathbf{f}(\mathbf{x}^{(i)})=\left(-dx_{i}+c_{1}x_{i}y_{i},ry_{i}-sy_{i}^{2}-c_{2}x_{i}y_{i}\right)$,
then
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$.
In Appendix E we show that the main equation (2) can now be written as

$\begin{dcases}\frac{dx_{i}}{dt}&=-dx_{i}+c_{1}x_{i}y_{i}+ac_{1}y_{i}\sum_{j}\mathcal{L}_{ij}x_{j}+(1-a)c_{1}x_{i}\sum_{j}\mathcal{L}_{ij}y_{j}\\\
\frac{dy_{i}}{dt}&=ry_{i}-sy_{i}^{2}-c_{2}x_{i}y_{i}-c_{2}y_{i}\sum_{j}\mathcal{L}_{ij}x_{j}\,,\end{dcases}$
(16)

where we can recognise the $n$ copies of the isolated Volterra system (15) and
the coupling due to the long-range interactions.

We can thus compute the explicit form of the polynomials $S_{1}(\xi)$ and
$S_{2}(\xi)$ as a function of the model parameters (see Appendix E) and
characterise the instability region defined by (9) as reported in Fig. 4 where
we show in the complex plane $(\Re\Lambda,\Im\Lambda)$ the regions for which
the instability condition is satisfied (grey), for a given set of parameters.
By using again a directed Erdős-Rényi network made of $50$ nodes and the
probability to create a directed link to be $0.5$ to describe the non-
reciprocal interactions, we can show the existence of eigenvalues
$\Lambda^{(\hat{\alpha})}$ (white dots in panel a) of Fig. 4) belonging to
this region, associate thus to a patchy solution (panel b) of Fig. 4). The
symmetric coupling being associated to a real spectrum is unable to drive the
system away from the homogenous equilibrium (panel c) of Fig. 4). Because of
the peculiar shape of the instability region, one can have parameters values
for which the instability region intersects the real axis (see Fig. 7 in
Appendix E). In this case also a symmetric coupling can trigger the
instability. From this example one can draw a general conclusion: if the
instability is possible by using symmetric long-range interactions, then the
same holds true for non-reciprocal ones; on the other hand, patterns resulting
from an instability due to non-symmetric interactions can never emerge if the
long-range interactions do reciprocate.

Figure 4: Instability region and patterns for the Volterra model. In panel a)
we report the region of the complex plane $(\Re\Lambda,\Im\Lambda)$ for which
the instability condition is satisfied (grey). For the chosen parameters
values ($c_{1}=2$, $c_{2}=13$, $r=1$, $s=1$, $d=0.4$ and $a=0.05$) we can
observe that the instability region does not intersect the real axis and thus
only non-reciprocal coupling can exhibit complex eigenvalues (white dots)
entering into the instability region and thus initiate the pattern as shown in
panel b) where we report $y_{i}(t)$ vs $t$ starting from initial conditions
close to the stable equilibrium. Any symmetric coupling determines real
eigenvalue (black dots in panel a)) that cannot give rise to the instability
as reported in the panel c) where we show $y_{i}(t)$ vs $t$. The underlying
coupling is a directed Erdős-Rényi network with $n=50$ nodes and a probability
for a direct link to exist between two nodes is $p=0.5$.

###  III.4 The FitzHugh-Nagumo model

The next system we consider is the FitzHugh-Nagumo model FitzHugh (1961);
Nagumo  _et al._ (1962); Rinzel and Keller (1973), a nonlinear system often
used as paradigm for the study of the emergence of Turing patterns Ohta
(1989); Shoji  _et al._ (2007); Carletti and Nakao (2020); Siebert  _et al._
(2020); Carletti and Nakao (2020) as well as for synchronisation phenomena
Vragovic  _et al._ (2006); Aqil  _et al._ (2012). The model has been conceived
in the framework of neuroscience as a schematisation of an electric impulse
propagating through an axon, for this reason, we believe that it would make a
suitable example for the long-range interactions we are describing in this
work.

The FitzHugh-Nagumo model can be described by the system of ODE

$\begin{dcases}\frac{du}{dt}&=\mu u-u^{3}-v\\\ \frac{dv}{dt}&=\gamma(u-\beta
v)\,,\end{dcases}$ (17)

where the parameters $\mu$, $\gamma$ and $\beta$ are assumed to be positive.
We will hereby focus on its behaviour close to the fixed point
$(u_{*},v_{*})=(0,0)$. The linear stability analysis ensures stability of the
latter under the conditions $\mu<\gamma\beta$ and $\mu\beta<1$. Let us observe
that, once such conditions are not met, the system undergoes a supercritical
Hopf-Andronov bifurcation Strogatz (2014): the equilibrium point becomes
unstable giving birth to a limit cycle solution. In this study we will limit
ourselves to the former case, leaving the oscillating case for a future work.
Let us also observe that this model is not a kinetic one, since the $-v$ term
appearing in the rate evolution for $u$, expresses a negative cross-effect
Szili and Tóth (1993); we can thus claim that pattern formation finds
applications beyond morphogenesis and chemical frameworks.

By setting $\mathbf{x}^{(i)}=(u_{i},v_{i})$ and
$\mathbf{f}(\mathbf{x}^{(i)})=(\mu u_{i}-u_{i}^{3}-v_{i},\gamma(u_{i}-\beta
v_{i}))$, we can cast the model in the framework presented in the previous
section by using the following coupling

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=(\mu
u_{i}-u_{i}^{3}-v_{j},\gamma(u_{j}-v_{i}))\,,$

where $\mathbf{x}^{(j)}=(u_{j},v_{j})$, that clearly satisfies
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$.
The main equation (2) becomes (see Appendix F)

$\begin{dcases}\frac{du_{i}}{dt}&=\mu
u_{i}-u_{i}^{3}-v_{i}-\sum_{j}\mathcal{L}_{ij}v_{j}\\\
\frac{dv_{i}}{dt}&=\gamma(u_{i}-\beta
v_{i})+\gamma\sum_{j}\mathcal{L}_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,,$ (18)

namely a set of $n$ FitzHugh-Nagumo models (17) coupled via long-range
connections in both variables.

We can then compute the polynomials $S_{1}(\xi)$ and $S_{2}(\xi)$ (see
Appendix F) and determine the instability region (grey) as reported in the
panel a) of Fig. 5. By using a directed Erdős-Rényi network made of $30$ nodes
and the probability to create a directed link to be $0.1$ to encode the non-
reciprocal interactions, we can show the existence of eigenvalues (white dots)
lying in the instability region and thus allowing for the instability onset
(panel b) of Fig. 5). Let observe that because of the shape of the instability
region, exhibiting a non-empty intersection with the negative real axis, there
can also be symmetric couplings for which patterns do emerge (panel c) of Fig.
5).

Figure 5: Instability region and patterns for the FitzHugh-Nagumo model. In
panel a) we report the region of the complex plane $(\Re\Lambda,\Im\Lambda)$
for which the instability condition is satisfied (grey). For the chosen
parameters values ($\alpha=1.5$, $\beta=0$, $\gamma=2.5$ and $\mu=0.01$) we
can observe that the instability region intersects the real axis and thus both
reciprocal (black dots) and non-reciprocal coupling (white dots) can exhibit
eigenvalues entering into the instability region and thus initiate the pattern
as shown in panels b) and c) where we report $u_{i}(t)$ vs $t$ starting from
initial conditions close to the stable equilibrium, $u_{*}=0$ for the
asymmetric coupling and the symmetric one. The underlying coupling is obtained
with a directed Erdős-Rényi network with $n=30$ nodes and a probability for a
direct link to exist between two nodes is $p=0.1$.

###  III.5 The Stuart-Landau model

As anticipated, the proposed method goes beyond the framework above presented
dealing with stationary homogeneous solution, but can be extended to study
systems exhibiting an oscillating behaviour, being the latter a regular or a
chaotic one. To support this claim, let us study for sake of definitiveness
the Stuart-Landau system (SL) Stuart and DiPrima (1978); Kuramoto (1984), as
paradigmatic model of nonlinear oscillators. The complex amplitude $w$ evolves
thus in time according to

$\frac{dw}{dt}=\sigma w-\beta|w|^{2}w\,,$ (19)

where $\sigma=\sigma_{\Re}+i\sigma_{\Im}$ and $\beta=\beta_{\Re}+i\beta_{\Im}$
are complex parameters, and we denote by $z_{\Re}$, resp. $z_{\Im}$, the real,
resp. imaginary, part of the complex number $z$. Let us observe that the SL
determines the normal form for a generic system close to a supercritical Hopf-
bifurcation, hence the results hereby presented are more general than the
specific model explored. Let us thus consider a system made of $n$ identical
SL oscillators and set the parameters such that each isolated oscillator
converges to the same limit cycle
$w_{LC}(t)=\sqrt{\sigma_{\Re}/\beta_{\Re}}e^{i\omega t}$,
$\omega=\sigma_{\Im}-\beta_{\Im}\sigma_{\Re}/\beta_{\Re}$. We are then
interested in determining the conditions responsible for the persistence of
this synchronous behaviour once the SL oscillators are allowed to interact
through non-reciprocal long-range interactions or, if on the contrary, an
instability sets up and drives the whole system to a new (spatially possibly)
heterogeneous state. Assume again the long-range interactions to be modelled
through a mean field ansatz (2), for instance

$\frac{dw_{j}}{dt}=\frac{\sigma}{k^{(in)}_{j}}\sum_{\ell}A_{j\ell}w_{\ell}-\beta
w_{j}|w_{j}|^{2}=\sigma w_{j}-\beta
w_{j}|w_{j}|^{2}+\sigma\sum_{\ell}\mathcal{L}_{j\ell}w_{\ell}\,,$ (20)

where $w_{j}$ is the complex state variable of the $j$-th SL system and
$A_{j\ell}$ encodes the non-reciprocal coupling. Such system admits an
homogeneous stable limit cycle solution if $\sigma_{\Re}>0$ and
$\beta_{\Re}>0$. We can then show the existence of non-reciprocal couplings
able to trigger the instability by destabilising the limit cycle solution,
eventually driving the system toward a new heterogeneous wavy solution.

Indeed, according to the theory hereby developed, we can always determine
model parameters allowing for an instability region in the complex plane (see
Appendix G), that does not intersect the real axis; the spectrum of a
reciprocal web of long-range interactions could thus never belong to the
instability region (black dots in Fig. 6) and any perturbation fades away. On
the other hand, the complex spectrum associated to non-reciprocal couplings
could intersect the instability region (white dots in Fig. 6), driving thus
the system toward the formation of patterns. We can thus state a claim similar
to the ones made above, namely non-reciprocal long-range interactions can more
easily drive the system toward a desynchronised state.

Figure 6: Instability region and patterns for the Stuart-Landau model. The
region of instability in the complex plane $(\Re\Lambda,\Im\Lambda)$ is
reported (grey) in panel a). We fixed the model parameters to the values
$\sigma_{\Re}=1$, $\sigma_{\Im}=4.3$, $\beta_{\Re}=1$ and $\beta_{\Im}=1$,
resulting into an instability region that does not intersect the real axis,
thus only an asymmetric coupling (white dots) can posses a complex spectrum
entering the instability region and initiate an heterogenous pattern as shown
in panel b) where we report $\Re w_{j}(t)$ vs time. Any symmetric coupling
determines real eigenvalue (black dots in panel a)) and thus it cannot
initiate the instability, in consequence all the units do synchronise (data
not shown). The underlying coupling is a directed Erdős-Rényi network with
$n=40$ nodes and a probability for a direct link to exist between two nodes is
$p=0.08$.

###  III.6 A sufficient condition for $d\geq 2$

As already anticipated, the $d=2$ dimensional systems have been presented for
sake of pedagogy, because they allow to obtain explicit and analytical
conditions for the onset of heterogeneous solutions, where the role of the
non-reciprocal links is clearly understood. We are now able to show that
similar results hold true in any dimensions, $d\geq 2$, and for generic
systems beyond the paradigmatic examples above presented. We will indeed prove
that given a system built using reciprocal interactions for which the
reference solution is stable and thus any small enough perturbation about it,
will fade away, then one can always find non-reciprocal interactions capable
to destabilise the reference solution, provided a condition on the model
parameters holds true, and can thus drive the system toward a new heterogenous
asymptotic state.

The starting point is thus the general system given by (4) and Eq. (7)
obtained by linearising the dynamics about the stationary reference solution,
$\mathbf{s}_{0}$, and then by projecting onto the Laplace eigenbasis,
$\phi^{(\alpha)}$, $\alpha=1,\dots,n$. The characteristic polynomial, whose
roots determine the stability feature of the reference solution, is now given
by:

$p_{d}(\lambda):=(-1)^{d}\left[a_{0}\lambda^{d}+a_{1}\lambda^{d-1}+a_{2}\lambda^{d-2}+\dots+a_{d}\right]\,,$
(21)

where the coefficients $a_{j}$, $j=0,\dots,d$, depend on the model parameters
and on the eigenvalues $\Lambda^{(\alpha)}$. A straightforward computation
allows to obtain for instance

$a_{0}=1\,,\quad a_{1}=-\mathrm{tr}\mathbf{J}^{(\alpha)}\quad\text{and}\quad
a_{d}=\det\mathbf{J}^{(\alpha)}\,,$

with
$\mathbf{J}^{(\alpha)}=\mathbf{J}(\mathbf{s}_{0})+\mathbf{J}_{2}(\mathbf{s}_{0})\Lambda^{(\alpha)}$;
it is thus clear that Eq. (8) is a particular case of the latter. A more
cumbersome computation (see Appendix H) allows to obtain an explicit formula
for $a_{2}$ that will be needed in the following. In the general
$d$-dimensional case, there are not explicit formulas for the roots of the
polynomial $p_{d}$ and even if they exist, their use to understand the role of
the model parameters and the impact of the eigenvalues $\Lambda^{(\alpha)}$
will be hopeless. We have thus to resort to a different approach based on the
Routh-Hurwitz stability criterion Routh (1877); Hurwitz (1895); Barnett
(1983), allowing to prove the (in)stability feature of a real coefficients
polynomial. In particular we will use the fact that a necessary condition to
have stable polynomial, i.e., a polynomial whose roots have negative real
parts, is that all its coefficients exhibit the same sign, in the present case
they should be positive being $a_{0}=1$.

Before to proceed further, we observe that the coefficients $a_{j}$ of the
polynomial $p_{d}$ are not (in general) real numbers, because of their
dependence on the complex $\Lambda^{(\alpha)}$, one cannot thus directly apply
the above criterium. To overcome this issue we introduce the polynomial
$q(\lambda)$ of degree $2d$ given by:

$q(\lambda):=p_{d}(\lambda)\overline{p_{d}(\overline{\lambda})}=\lambda^{2d}+b_{1}\lambda^{2d-1}+b_{2}\lambda^{2d-2}+\dots+b_{2d}\,,$

where $\overline{z}$ denotes the complex conjugate of the complex number $z$.
One can prove (see Appendix H) that the roots of this polynomial have the same
real parts of the roots of $p_{d}$ (counted twice) and moreover its
coefficients become real numbers, whose expression can be related to the
coefficients of $p_{d}$.

The proof of our statement proceeds thus as follows. Assume to have a system
whose nodes do interact through a symmetric ensemble of links and fix the
model parameters in such a way the reference solution $\mathbf{s}_{0}$ is
stable, with respect to heterogeneous perturbations. Hence for all
$\alpha=1,\dots,n$, all the roots $\lambda_{i}$, $i=1,\dots,d$, of $p_{d}$
have negative real parts and so do the ones of $q$. Invoking the necessary
condition of the Routh-Hurwitz stability criterion we can conclude that the
coefficients of $q$ are positive numbers whenever we set
$\Im\Lambda^{(\alpha)}=0$, for all $\alpha$, because the Laplace spectrum is
real (see Appendix A):

$b_{j}\rvert_{\Im\Lambda^{(\alpha)}=0}>0\quad\forall j=1,\dots,2d\,.$

We can then prove (we refer to Appendix H for more details) that if the model
parameters satisfy the condition

$\mathrm{tr}\mathbf{J}^{2}_{2}<0\,,$ (22)

then we can always find a non-reciprocal coupling such that the imaginary part
of the spectrum does satisfy the condition

$\exists\hat{\alpha}>1\text{
s.t.}\left(\Im\Lambda^{(\hat{\alpha})}\right)^{2}>\frac{\left(2\Re
a_{2}+\left[\mathrm{tr}\mathbf{J}^{(\hat{\alpha})}\right]^{2}\right)\rvert_{\Im\Lambda^{(\hat{\alpha})}=0}}{|\mathrm{tr}\mathbf{J}_{2}^{2}|}\,,$
(23)

ensuring that the coefficient $b_{2}$ is negative and thus the polynomial $q$
is not stable. It admits thus at least one root, the one associated to
$\hat{\alpha}$, whose real part is positive, and so does $p_{d}$. In
conclusion the reference solution $\mathbf{s}_{0}$ is unstable, any initial
perturbation is amplified and the system moves away from the homogeneous
solution. Let us observe that we can always build a web of non-reciprocal
interactions satisfying the latter conditions, by using the ideas developed in
Nicoletti  _et al._ (2021) to generate asymmetric Laplacian matrices, and thus
directed networks, that possess prescribed spectra.

##  IV Conclusion

In conclusion, we have proposed and analysed a mechanism for pattern
formation, being the latter stationary or time dependent, rooted on the
presence of non-reciprocal long-range interactions, to show that self-
organisation can manifest without the presence of diffusion. Long-range non-
reciprocal interactions, modelled with a mean-field scheme, return a consensus
Laplace operator, whose complex spectrum sets the conditions for the
instability and the ensuing spatial or temporal patterns. The intrinsic
asymmetries within the system, described by non-reciprocal interactions, are
the key factor in the disruption of the homogeneous solution and thus the
driver for the diversity of patterns of complexity one can observe in Nature.
The theory has been applied to several paradigmatic models used in the
literature to study the emergence of patterns or synchronisation. The proposed
mechanism complements thus the Turing one Turing (1952); Nakao and Mikhailov
(2010) and it is suitable for all phenomena where patterns emergence is not
driven by a diffusive process but immobile species interact through non-
reciprocal long-range couplings.

Acknowledgements
We thank Luca Gallo and Jean-François de Kemmeter for discussions. RM is
supported by an FRIA-FNRS PhD fellowship, Grant FC 33443, funded by the
Walloon region.

## References

  * Newman (2010) M. E. Newman,  _Networks: An Introduction_ (Oxford University Press, Oxford, 2010).
  * Barabási (2016) A.-L. Barabási,  _Network science_ (Cambridge university press, 2016).
  * Pastor-Satorras and Vespignani (2010) R. Pastor-Satorras and A. Vespignani, Nature Physics 6, 480 (2010).
  * Nicolis and Prigogine (1977) G. Nicolis and I. Prigogine,  _Self-organization in nonequiibrium systems: From dissipative structures to order through fluctuations_ (J. Wiley and Sons, 1977).
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization_ (Cambridge University Press, Cambridge, UK, 2001).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics reports 469, 93 (2008).
  * Bressloff  _et al._ (2001) P. C. Bressloff, C. J. D, G. Martin, T. P. J, and W. M. C, Phil. Trans. R. Soc. Lond. B 356, 299 (2001).
  * Bressloff and D (2003) P. C. Bressloff and C. J. D, J. of Physiology Paris 97, 221 (2003).
  * Strogatz (2003) S. Strogatz,  _Sync : The Emerging Science of Spontaneous Order_ (Hyperion, 2003).
  * Nakamasu  _et al._ (2009) A. Nakamasu, G. Takahashi, A. Kanbe, and S. Kondo, Proceedings of the National Academy of Sciences 106, 8429 (2009).
  * Bullara and De Decker (2015) D. Bullara and Y. De Decker, Nature Communications 6, 6971 (2015).
  * Turing (1952) A. M. Turing, Phil. Trans. R. Soc. Lond. B 237, 37 (1952).
  * Gizynski and Gorecki (2017) K. Gizynski and J. Gorecki, Phys. Chem. Chem. Phys 19, 6519 (2017).
  * Voorsluijs  _et al._ (2017) V. Voorsluijs, I. G. Kevrekidis, and Y. De Decker, Phys. Chem. Chem. Phys. 19, 22528 (2017).
  * Parrilla-Gutierrez  _et al._ (2020) J. M. Parrilla-Gutierrez, A. Sharma, S. Tsuda, G. J. T. Cooper, G. Aragon-Camarasa, K. Donkers, and L. Cronin, Nature Communications 11, 1442 (2020).
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Phys. Rev. Letters 80, 2109 (1998).
  * Cencetti  _et al._ (2020) G. Cencetti, F. Battiston, T. Carletti, and D. Fanelli, Chaos, Solitons and Fractals 134, 109707 (2020).
  * Bascompte  _et al._ (2003) J. Bascompte, P. Jordano, C. J. Melián, and J. M. Olesen, Proceedings of the National Academy of Sciences 100, 9383 (2003).
  * Abramson  _et al._ (2011) G. Abramson, C. A. Trejo Soto, and L. Oña, PLoS ONE 6, e21028 (2011).
  * Hollingdale  _et al._ (2018) E. Hollingdale, F. J. Pérez-Barbería, and D. McPetrie Walker, PLoS ONE 13, e0208202 (2018).
  * Zhang  _et al._ (2019) Y. Zhang, T. K. Tsang, E. A. Bushong, L.-A. Chu, A.-S. Chiang, M. H. Ellisman, J. Reingruber, and C.-Y. Su, Nature Communications 10, 1560 (2019).
  * Nakao and Mikhailov (2010) H. Nakao and A. S. Mikhailov, Nature Physics 6, 544 (2010).
  * Pecora  _et al._ (1997) L. M. Pecora, T. L. Carroll, G. A. Johnson, D. J. Mar, and J. F. Heagy, Chaos: An Interdisciplinary Journal of Nonlinear Science 7, 520 (1997).
  * Ghosh  _et al._ (2014) R. Ghosh, S.-h. Teng, K. Lerman, and X. Yan, in  _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_ (ACM, 2014) pp. 1406–1415.
  * Lambiotte  _et al._ (2011) R. Lambiotte, R. Sinatra, J.-C. Delvenne, T. S. Evans, M. Barahona, and V. Latora, Physical Review E 84, 017102 (2011).
  * Olfati-Saber  _et al._ (2007) R. Olfati-Saber, J. A. Fax, and R. M. Murray, Proceedings of the IEEE 95, 215 (2007).
  * Krause (2009) U. Krause, Elemente der Mathematik 64, 1 (2009).
  * Geršgorin (1931) S. Geršgorin, Bulletin de l’Académie des Sciences de l’URSS. Classe des sciences mathématiques et na , 749 (1931).
  * Huang  _et al._ (2009) L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, Phys. Rev. E 80, 036204 (2009).
  * Asllani  _et al._ (2014a) M. Asllani, J. D. Challenger, F. S. Pavone, L. Sacconi, and D. Fanelli, Nature Communications 5, 4517 (2014a).
  * Asllani  _et al._ (2020) M. Asllani, T. Carletti, D. Fanelli, and P. K. Maini, The European Physical Journal B 93, 1 (2020).
  * Prigogine and Nicolis (1967) I. Prigogine and G. Nicolis, J. Chem. Phys. 46, 3542 (1967).
  * Prigogine and Lefever (1968) I. Prigogine and R. Lefever, J. Chem. Phys. 48, 1695 (1968).
  * Muolo  _et al._ (2021) R. Muolo, T. Carletti, J. P. Gleeson, and M. Asllani, Entropy 23, 36 (2021).
  * Asllani  _et al._ (2014b) M. Asllani, D. M. Busiello, T. Carletti, D. Fanelli, and G. Planchon, Phys. Rev. E 90, 042814 (2014b).
  * Asllani  _et al._ (2016) M. Asllani, T. Carletti, and D. Fanelli, Eur. Phys. J. B , 89 (2016).
  * Petit  _et al._ (2017) J. Petit, B. Lauwens, D. Fanelli, and T. Carletti, Phys. Rev. Letters 119, 148301 (2017).
  * Koseska  _et al._ (2013) A. Koseska, E. Volkov, and J. Kurths, Physics reports 531, 173 (2013).
  * Lucas  _et al._ (2018) M. Lucas, D. Fanelli, T. Carletti, and J. Petit, EPL (Europhysics Letters) 121, 58008 (2018).
  * Mimura and Murray (1978) M. Mimura and J. D. Murray, J. Theor. Biol. 75, 249 (1978).
  * Challenger  _et al._ (2015) J. Challenger, D. Fanelli, and R. Burioni, Phys. Rev. E 92, 022818 (2015).
  * Asllani  _et al._ (2015) M. Asllani, D. M. Busiello, T. Carletti, D. Fanelli, and G. Planchon, Scientific Reports 5, 12927 (2015).
  * McKane and Newman (2005) A. J. McKane and T. J. Newman, Physical review letters 94, 218102 (2005).
  * Layek and Pati (2019) G. Layek and N. Pati, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 093104 (2019).
  * Pati  _et al._ (2020) N. Pati, G. Layek, and N. Pal, Chaos, Solitons and Fractals 140, 110184 (2020).
  * FitzHugh (1961) R. FitzHugh, Biophys. J. 1, 445 (1961).
  * Nagumo  _et al._ (1962) J. Nagumo, S. Arimoto, and S. Yoshizawa, Proc. IRE 50, 2061 (1962).
  * Rinzel and Keller (1973) J. Rinzel and J. B. Keller, Biophys. J. 13, 1313 (1973).
  * Ohta (1989) T. Ohta, Prog. Theoret. Phys. Suppl. 99, 425 (1989).
  * Shoji  _et al._ (2007) H. Shoji, K. Yamada, D. Ueyama, and T. Ohta, Phys. Rev. E 75, 046212 (2007).
  * Carletti and Nakao (2020) T. Carletti and H. Nakao, Phys. Rev. E 101, 022203 (2020).
  * Siebert  _et al._ (2020) B. A. Siebert, C. L. Hall, J. P. Gleeson, and M. Asllani, Phys. Rev. E 102, 052306 (2020).
  * Vragovic  _et al._ (2006) I. Vragovic, E. Louis, and A. Díaz-Guilera, Europhys. Letters 76, 780 (2006).
  * Aqil  _et al._ (2012) M. Aqil, K.-S. Hong, and M.-Y. Jeong, Communications in Nonlinear Science and Numerical Simulation 17, 1615 (2012).
  * Strogatz (2014) S. Strogatz,  _Nonlinear Dynamics and Chaos_ (Avalon Publishing, 2014).
  * Szili and Tóth (1993) L. Szili and J. Tóth, Phys. Rev. E 48, 183 (1993).
  * Stuart and DiPrima (1978) J. T. Stuart and R. C. DiPrima, Proc R Soc Lond A 362, 27 (1978).
  * Kuramoto (1984) Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Springer-Verlag, New York, 1984).
  * Routh (1877) E. J. Routh,  _Stability of a given state of motion_ (MacMillan and Co, London, 1877).
  * Hurwitz (1895) A. Hurwitz, Math. Ann. 46, 273 (1895).
  * Barnett (1983) S. Barnett,  _Polynomials and Linear Control Systems_ (Dekker, New York, 1983).
  * Nicoletti  _et al._ (2021) S. Nicoletti, T. Carletti, D. Fanelli, G. Battistelli, and L. Chisci, Journal of complex networks 2, 015004 (2021).
  * Gambuzza  _et al._ (2021) L. V. Gambuzza, F. Di Patti, G. Luca, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, Nature Communication 12, 1255 (2021).
  * Nishikawa and Motter (2006) T. Nishikawa and A. E. Motter, Physical Review E 73, 065106R (2006).

##  Appendix A About the spectrum of the reactive Laplacian in case of
reciprocal interactions.

Let us assume the network of interactions to be symmetric, $A_{ij}=A_{ji}$,
for all $i$ and $j$, and let us recall the definition of the reactive
Laplacian $\mathcal{L}_{ij}=A_{ij}/k_{i}-\delta_{ij}$, where
$k_{i}=\sum_{j}A_{ij}$, is the node degree. Let us introduce the symmetric
Laplace matrix $\mathcal{L}^{sym}_{ij}=A_{ij}/\sqrt{k_{i}k_{j}}-\delta_{ij}$
and observe that its eigenvalues are real and negative but the null one. Let
$\mathbf{D}$ be the diagonal matrix with the nodes degree on the diagonal,
then

$\mathcal{L}=\mathbf{D}^{-1}\mathbf{A}-\mathbf{I}=\mathbf{D}^{-1/2}\left[\mathbf{D}^{-1/2}\mathbf{A}-\mathbf{I}\right]\mathbf{D}^{1/2}=\mathbf{D}^{-1/2}\mathcal{L}^{sym}\mathbf{D}^{1/2}\,.$

Namely $\mathcal{L}$ and $\mathcal{L}^{sym}$ are similar matrices and thus
exhibit the same set of eigenvalues.

##  Appendix B Linear stability analysis

Let us consider the homogeneous reference solution $\mathbf{s}(t)$ and a
spatially dependent perturbation about the latter,
$\mathbf{x}^{(i)}(t)=\mathbf{s}(t)+\mathbf{u}^{(i)}(t)$, then by inserting
this information into Eq. (2) (main text) and by retaining only the linear
terms in $\mathbf{u}^{(i)}$ we obtain

$\displaystyle\frac{d\mathbf{u}^{(i)}}{dt}=\frac{d\mathbf{x}^{(i)}}{dt}-\frac{d\mathbf{s}}{dt}$
$\displaystyle=$
$\displaystyle\frac{1}{k^{(in)}_{i}}{\sum_{j}A_{ij}\mathbf{F}(\mathbf{s}+\mathbf{u}^{(i)},\mathbf{s}+\mathbf{u}^{(j)})}-\mathbf{f}(\mathbf{s})$
$\displaystyle=$
$\displaystyle\frac{1}{k^{(in)}_{i}}{\sum_{j}A_{ij}\left(\sum_{\ell}\partial_{x_{\ell}^{(i)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(i)}+\sum_{\ell}\partial_{{x}_{\ell}^{(j)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(j)}\right)}$
$\displaystyle=$
$\displaystyle\sum_{\ell}\partial_{x_{\ell}^{(i)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(i)}+\frac{1}{k^{(in)}_{i}}{\sum_{j}A_{ij}\sum_{\ell}\partial_{{x}_{\ell}^{(j)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(j)}}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{1}\mathbf{u}^{(i)}+\frac{1}{k^{(in)}_{i}}{\sum_{j}A_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}}\quad\forall
i=1,\dots,n\,,$

where we recall that
$\mathbf{J}_{1}=\partial_{\mathbf{x}_{1}}\mathbf{F}(\mathbf{x}_{1},\mathbf{x}_{2})_{(\mathbf{s},\mathbf{s})}$
and
$\mathbf{J}_{2}=\partial_{\mathbf{x}_{2}}\mathbf{F}(\mathbf{x}_{1},\mathbf{x}_{2})_{(\mathbf{s},\mathbf{s})}$.
By observing that $\mathbf{f}(\mathbf{x})=\mathbf{F}(\mathbf{x},\mathbf{x})$
one can prove that
$\partial_{\mathbf{x}}\mathbf{f}:=\mathbf{J}=\mathbf{J}_{1}+\mathbf{J}_{2}$.
Hence by slightly rewriting the previous equation, we obtain

$\displaystyle\frac{d\mathbf{u}^{(i)}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mathbf{J}_{1}\mathbf{u}^{(i)}+\frac{1}{k^{(in)}_{i}}{\sum_{j}A_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}}=\mathbf{J}_{1}\mathbf{u}^{(i)}+\mathbf{J}_{2}\mathbf{u}^{(i)}+\sum_{j}\left(\frac{A_{ij}}{k^{(in)}_{i}}-\delta_{ij}\right)\mathbf{J}_{2}\mathbf{u}^{(j)}$
(24) $\displaystyle=$
$\displaystyle\left(\mathbf{J}_{1}+\mathbf{J}_{2}\right)\mathbf{u}^{(i)}+\sum_{j}\mathcal{L}_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}$
$\displaystyle=$
$\displaystyle\mathbf{J}\mathbf{u}^{(i)}+\sum_{j}\mathcal{L}_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}\,,$

where we used the matrix $\mathcal{L}$, given by Eq. (3) (main text).

By introducing the $n\times d$ vector
$\mathbf{u}=((\mathbf{u}^{(1)})^{\top},\dots,(\mathbf{u}^{(n)})^{\top})^{\top}$,
we can eventually rewrite the latter equation in a compact form as:

$\frac{d\mathbf{u}}{dt}(t)=\left[\mathbf{I}_{n}\otimes\left(\mathbf{J}_{1}+\mathbf{J}_{2}\right)+\mathcal{L}\otimes\mathbf{J}_{2}\right]\mathbf{u}=\left[\mathbf{I}_{n}\otimes\mathbf{J}+\mathcal{L}\otimes\mathbf{J}_{2}\right]\mathbf{u}\,,$
(25)

where $\otimes$ denotes the Kronecker product of matrices. In this way we
emphasise the role of the Jacobian of the isolated system
$\mathbf{I}_{n}\otimes\mathbf{J}$ and the one for the coupling part,
$\mathcal{L}\otimes\mathbf{J}_{2}$, that vanishes once we assume
$\mathbf{A}=\mathbf{I}_{n}$.

Eq. (24) is a linear system involving matrices with size $nd\times nd$. One
can reduce the complexity of the latter by assuming the existence of an
eigenbasis of $\mathcal{L}$, $\phi^{(\alpha)}$, $\alpha=1,\dots,n$, with
associated eigenvalue $-2<\Re\Lambda^{(\alpha)}\leq 0$. Then by rewriting
$\mathbf{u}^{(i)}=\sum_{\beta}\mathbf{u}^{(\beta)}\phi^{(\beta)}_{i}$ and
inserting the latter into (24) we obtain

$\displaystyle\sum_{\beta}\frac{d\mathbf{u}^{(\beta)}}{dt}\phi_{i}^{(\beta)}$
$\displaystyle=$
$\displaystyle\sum_{\beta}\mathbf{J}(\mathbf{s}(t))\mathbf{u}^{(\beta)}\phi_{i}^{(\beta)}+\sum_{j}\sum_{\beta}\mathcal{L}_{ij}\mathbf{J}_{2}(\mathbf{s}(t))\mathbf{u}^{(\beta)}\phi_{j}^{(\beta)}$
$\displaystyle=$
$\displaystyle\sum_{\beta}\mathbf{J}(\mathbf{s}(t))\mathbf{u}^{(\beta)}\phi_{i}^{(\beta)}+\sum_{\beta}\Lambda^{(\beta)}\mathbf{J}_{2}(\mathbf{s}(t))\mathbf{u}^{(\beta)}\phi_{i}^{(\beta)}\,.$

By multiplying by the left eigenvectors $\psi_{i}^{(\alpha)}$ and summing over
the index $i$, recall that
$\sum_{i}\psi_{i}^{(\alpha)}\phi_{i}^{(\beta)}=\delta_{\alpha\beta}$, we can
obtain Eq. (6) (main text), namely

$\frac{d\mathbf{u}^{(\alpha)}}{dt}=\mathbf{J}(\mathbf{s}(t))\mathbf{u}^{(\alpha)}+\Lambda^{(\alpha)}\mathbf{J}_{2}(\mathbf{s}(t))\mathbf{u}^{(\alpha)}\quad\forall\alpha=1,\dots,n\,.$

Assuming to deal with $d=2$ dimensional systems and with a stationary
reference solution, i.e., $\mathbf{s}(t)=\mathbf{s}_{0}$ for all $t$; one can
then realise that the eigenvalues of $\mathbf{J}^{(\alpha)}$, $\lambda_{i}$,
are the solutions of the second order equation

$\lambda_{i}^{2}-\mathrm{tr}\mathbf{J}^{(\alpha)}\lambda_{i}+\det\mathbf{J}^{(\alpha)}=0\quad\forall
i=1,\dots,n\,.$

Note that if on the other hand, we were interested in studying $d>2$ dynamical
systems, then we should consider roots of $d$-degree polynomial equations for
which there is no general closed form solution and thus one has to recur to
numerics to determine the instability region or to another ideas as we have
shown in the main text and detailed in Appendix H.

Back to the $d=2$ case, one can express the real part of the above roots as:

$\Re\lambda_{i}=\frac{1}{2}\left(\Re\mathrm{tr}\mathbf{J}^{(\alpha)}+\gamma\right)\,,$

where

$\gamma=\sqrt{\frac{A+\sqrt{A^{2}+B^{2}}}{2}}\,,A=\left(\Re\mathrm{tr}\mathbf{J}^{(\alpha)}\right)^{2}-\left(\Im\mathrm{tr}\mathbf{J}^{(\alpha)}\right)^{2}-4\Re\mathrm{det}\mathbf{J}^{(\alpha)}\text{
and
}B=2\Re\mathrm{tr}\mathbf{J}^{(\alpha)}\Im\mathrm{tr}\mathbf{J}^{(\alpha)}-4\Im\mathrm{det}\mathbf{J}^{(\alpha)}\,.$

A straightforward but lengthy computation allows to rewrite the condition for
instability Eq. (9) (main text), in terms of two polynomials, $S_{2}$ of
second degree and $S_{1}$ of fourth degree. More precisely,
$S_{2}(\xi)=c_{2,2}\xi^{2}+c_{2,1}\xi+c_{2,0}$ with coefficients

$\displaystyle c_{2,2}$ $\displaystyle=$
$\displaystyle-\det\mathbf{J}_{2}(4\det\mathbf{J}_{2}-\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2})$
$\displaystyle c_{2,1}$ $\displaystyle=$
$\displaystyle-\Delta_{1}(4\det\mathbf{J}_{2}-\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2})$
(26) $\displaystyle c_{2,0}$ $\displaystyle=$
$\displaystyle-\Delta_{1}^{2}+\Delta_{1}\mathrm{tr}\mathbf{J}\mathrm{tr}\mathbf{J}_{2}-\det\mathbf{J}_{2}\left(\mathrm{tr}\mathbf{J}\right)^{2}\,,$

and
$S_{1}(\xi)=c_{1,4}\xi^{4}+c_{1,3}\xi^{3}+c_{1,2}\xi^{2}+c_{1,1}\xi+c_{1,0}$
with coefficients

$\displaystyle c_{1,4}$ $\displaystyle=$
$\displaystyle\det\mathbf{J}_{2}\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2}$
$\displaystyle c_{1,3}$ $\displaystyle=$
$\displaystyle\mathrm{tr}\mathbf{J}_{2}(\Delta_{1}\mathrm{tr}\mathbf{J}_{2}+2\det\mathbf{J}_{2}\mathrm{tr}\mathbf{J})$
$\displaystyle c_{1,2}$ $\displaystyle=$
$\displaystyle\det\mathbf{J}\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2}+2\Delta_{1}\mathrm{tr}\mathbf{J}\mathrm{tr}\mathbf{J}_{2}+\det\mathbf{J}_{2}\left(\mathrm{tr}\mathbf{J}\right)^{2}$
(27) $\displaystyle c_{1,1}$ $\displaystyle=$
$\displaystyle\mathrm{tr}\mathbf{J}(2\det\mathbf{J}\mathrm{tr}\mathbf{J}_{2}+\Delta_{1}\mathrm{tr}\mathbf{J})$
$\displaystyle c_{1,0}$ $\displaystyle=$
$\displaystyle\det\mathbf{J}\left(\mathrm{tr}\mathbf{J}\right)^{2}\,,$

where we introduced
$\Delta_{1}=J_{2,11}J_{22}-J_{2,12}J_{21}-J_{2,21}J_{12}+J_{2,22}J_{11}$.

Let us now consider the general case in which the existence of an eigenbasis
of $\mathcal{L}$ is not guaranteed, namely the Laplace matrix is defective. We
can then invoke the Jordan canonical form theorem to determine an invertible
$n\times n$ matrix $\mathbf{P}$ such that

$\mathbf{P}^{-1}\mathcal{L}\mathbf{P}=\mathbf{B}=\mathrm{diag}(\mathbf{B}_{1},\dots,\mathbf{B}_{\ell})\,,$

where the $\mathbf{B}_{j}$ is the $m_{j}\times m_{j}$ Jordan block,
$m_{1}+\dots+m_{\ell}=n$,

$\mathbf{B}_{j}=\left(\begin{matrix}\Lambda^{(j)}&&&\\\ 1&\Lambda^{(j)}&&\\\
&\ddots&\ddots&\\\ &&1&\Lambda^{(j)}\end{matrix}\right)\,,$

being $\Lambda^{(j)}$ the $j$-th eigenvalue of $\mathcal{L}$. Because
$\Lambda^{(1)}=0$ we also have $\mathbf{B}_{1}=0$.

Let us consider again Eq. (25). By defining
$\mathbf{Q}=\mathbf{P}\otimes\mathbf{I}_{d}$ and
$\mathbf{v}=\mathbf{Q}^{-1}\mathbf{u}$ we get

$\frac{d\mathbf{v}}{dt}(t)=\mathbf{Q}^{-1}\frac{d\mathbf{u}}{dt}(t)=(\mathbf{P}^{-1}\otimes\mathbf{I}_{d})\left[\mathbf{I}_{n}\otimes\mathbf{J}+\mathcal{L}\otimes\mathbf{J}_{2}\right](\mathbf{P}\otimes\mathbf{I}_{d})\mathbf{Q}^{-1}\mathbf{u}=\left[\mathbf{I}_{n}\otimes\mathbf{J}+\mathbf{B}\otimes\mathbf{J}_{2}\right]\mathbf{v}\,.$
(28)

The vector $\mathbf{v}$ inherits the Jordan decomposition, hence we can write
$\mathbf{v}=((\mathbf{v}^{(1)})^{\top},\dots,(\mathbf{v}^{(\ell)})^{\top})^{\top}$,
where $\mathbf{v}^{(j)}$ is a $(d\times m_{j})$-dimensional vector. Eq. (28)
can thus be rewritten as

$\displaystyle\frac{d\mathbf{v}^{(1)}}{dt}(t)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\mathbf{J}&&\\\ &\ddots&\\\
&&\mathbf{J}\end{matrix}\right)\mathbf{v}^{(1)}\,,$ (29)
$\displaystyle\frac{d\mathbf{v}^{(j)}}{dt}(t)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\mathbf{J}&&\\\ &\ddots&\\\
&&\mathbf{J}\end{matrix}\right)\mathbf{v}^{(j)}+\left(\begin{matrix}\Lambda^{(j)}\mathbf{J}_{2}&&&\\\
\mathbf{J}_{2}&\Lambda^{(j)}\mathbf{J}_{2}&&\\\ &\ddots&\ddots&\\\
&&\mathbf{J}_{2}&\Lambda^{(j)}\mathbf{J}_{2}\end{matrix}\right)\mathbf{v}^{(j)}\quad\forall
j=2,\dots,\ell\,.$ (30)

The first Eq. (29) returns again the stability condition for the isolated
systems, which is a working assumption for the developed framework. Let us now
consider the generic Eq. (30) and by writing
$\mathbf{v}^{(j)}=((\mathbf{v}^{(j)}_{1})^{\top},\dots,(\mathbf{v}^{(j)}_{m_{j}})^{\top})^{\top}$,
where $\mathbf{v}^{(j)}_{i}\in\mathbf{R}^{d}$ for all $i=1,\dots,m_{j}$, we
get

$\displaystyle\frac{d\mathbf{v}^{(j)}_{1}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mathbf{J}\mathbf{v}^{(j)}_{1}+\Lambda^{(j)}\mathbf{J}_{2}\mathbf{v}^{(j)}_{1}\,,$
(31) $\displaystyle\frac{d\mathbf{v}^{(j)}_{2}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mathbf{J}\mathbf{v}^{(j)}_{2}+\Lambda^{(j)}\mathbf{J}_{2}\mathbf{v}^{(j)}_{2}+\mathbf{J}_{2}\mathbf{v}^{(j)}_{1}\,,$
$\displaystyle\vdots$ $\displaystyle\frac{d\mathbf{v}^{(j)}_{m_{j}}}{dt}(t)$
$\displaystyle=$
$\displaystyle\mathbf{J}\mathbf{v}^{(j)}_{m_{j}}+\Lambda^{(j)}\mathbf{J}_{2}\mathbf{v}^{(j)}_{m_{j}}+\mathbf{J}_{2}\mathbf{v}^{(j)}_{m_{j}-1}\,.$
(33)

To simplify the following analysis, we will assume to deal with a stationary
reference solution $\mathbf{s}(t)=\mathbf{s}_{0}$, hence all the involved
matrices are constant ones. Hence the first Eq. (31) is the analogous of Eq.
(6) (main text) and one can determine a condition on $\Lambda^{(j)}$ to make
the projection $\mathbf{v}_{1}^{(j)}$ unstable; in the case $d=2$ this
accounts to perform the analysis above presented returning the condition (9),
where $\hat{\alpha}$ is replaced by $j$.

Let us now consider the second Eq. (B) and observe that it is composed by two
terms, the first one involves the same matrix of the first equation,
$\mathbf{J}+\Lambda^{(j)}\mathbf{J}_{2}$, while the second one depends on the
projection $\mathbf{v}_{1}^{(j)}(t)$. If, for the choice of $\Lambda^{(j)}$
the matrix is unstable and thus $\mathbf{v}_{1}^{(j)}(t)$ has an exponential
growth, then the same is true for $\mathbf{v}_{2}^{(j)}(t)$. By considering
the remaining equations and by exploiting the peculiar shape of the system,
this allows to prove that if the first Eq. (31) returns an unstable solution,
then all the solutions $\mathbf{v}_{i}^{(j)}(t)$ are unstable as well.

In conclusion if the Laplace matrix $\mathcal{L}$ is defective, one can check
the instability condition on the available eigenvalues and conclude about the
emergence of patterns solely based on this information. Let us observe that
this is a sufficient condition, indeed it can happen that the matrix
$\mathbf{J}+\Lambda^{(j)}\mathbf{J}_{2}$ is stable, but the presence of Jordan
blocks introduces a transient (polynomial) growth in the linear regime that
results strong enough to limit the validity of the linear approximation. Thus
the nonlinear system could exhibit orbits departing from the homogeneous
reference solution, only infinitesimal perturbations will be attracted to the
latter. Stated differently, the stability basin of the reference solution
shrinks considerably in presence of defective Laplace matrix, the solution is
thus stable but finite perturbations can be amplified.

##  Appendix C Analysis of the Brusselator model

Let us consider the two species Brusselator model,
$\mathbf{f}(\mathbf{x}^{(i)})=(1-(b+1)u_{i}+cu_{i}^{2}v_{i},bu_{i}-cu_{i}^{2}v_{i})$,
where $\mathbf{x}^{(i)}=(u_{i},v_{i})$, and the non-local coupling given by

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=(1-(b+1)u_{i}+cu_{i}^{2}v_{i},bu_{j}-cu_{i}^{2}v_{i})\,,$

where $\mathbf{x}^{(j)}=(u_{j},v_{j})$. Eq. (2) can thus be rewritten as

$\begin{dcases}\frac{du_{i}}{dt}&=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}\left(1-(b+1)u_{i}+cu_{i}^{2}v_{i}\right)\\\
\frac{dv_{i}}{dt}&=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}\left(bu_{j}-cu_{i}^{2}v_{i}\right)\end{dcases}\quad\forall
i=1,\dots,n\,,$

by recalling the definition of $k_{i}^{(in)}=\sum_{j}A_{ij}$ and by observing
that the first equation does not depend on the index $j$, it can be
straightforwardly simplified to return

$\frac{du_{i}}{dt}=1-(b+1)u_{i}+cu_{i}^{2}v_{i}\quad\forall i=1,\dots,n\,.$

Let us now consider the second equation and isolate the terms not depending on
the index $j$

$\frac{dv_{i}}{dt}=-cu_{i}^{2}v_{i}+\frac{b}{k_{i}^{(in)}}\sum_{j}A_{ij}u_{j}=-cu_{i}^{2}v_{i}+b\sum_{j}\left(\frac{A_{ij}}{k_{i}^{(in)}}-\delta_{ij}\right)u_{j}+bu_{i}=bu_{i}-cu_{i}^{2}v_{i}+b\sum_{j}\mathcal{L}_{ij}u_{j}\quad\forall
i=1,\dots,n\,.$

In conclusion we obtain

$\begin{dcases}\frac{du_{i}}{dt}&=1-(b+1)u_{i}+cu_{i}^{2}v_{i}\\\
\frac{dv_{i}}{dt}&=bu_{i}-cu_{i}^{2}v_{i}+b\sum_{j}\mathcal{L}_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,,$

From the definitions of $\mathbf{f}$ and $\mathbf{F}$ we can compute the
associated Jacobian matrices:

$\mathbf{J}_{\mathrm{Bxl}}=\left(\begin{matrix}b-1&c\\\
-b&-c\end{matrix}\right)\text{ and }\mathbf{J}_{2}=\left(\begin{matrix}0&0\\\
b&0\end{matrix}\right)\,,$

and thus, by using the formulas (B) and (B) for the coefficients $c_{1,i}$,
$i=0,1,2,3,4$, and $c_{2,i}$, $i=0,1,2$, we obtain the polynomials $S_{1}$ and
$S_{2}$:

$S_{1}(\xi)=c(b-c-1)^{2}(1-b\xi)\text{ and }S_{2}(\xi)=-b^{2}c^{2}\,.$

The condition the complex eigenvalues $\Lambda^{(\alpha)}$ of $\mathcal{L}$
have to satisfy to induce an instability is eventually given by

$(\Im\Lambda^{(\alpha)})^{2}\geq\frac{(b-c-1)^{2}}{cb^{2}}\left(1-b\Re\Lambda^{(\alpha)}\right)\,.$

##  Appendix D Analysis of the Mimura-Murray model

Let us consider the Mimura-Murray model with reaction terms,
$\mathbf{f}(\mathbf{x}^{(i)})=\left(u_{i}\left(\frac{a+bu_{i}-u_{i}^{2}}{c}-v_{i}\right),v_{i}\left(u_{i}-(1+dv_{i})\right)\right)$,
where $\mathbf{x}^{(i)}=(u_{i},v_{i})$. The chosen non-local coupling is

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=\left(u_{i}\left(\frac{a+bu_{j}-u_{i}^{2}}{c}-v_{j}\right),v_{i}\left(u_{j}-(1+dv_{i})\right)\right)\,.$

From Eq. (2) we obtain

$\begin{dcases}\frac{du_{i}}{dt}&=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}u_{i}\left(\frac{a+bu_{j}-u_{i}^{2}}{c}-v_{j}\right)\\\
\frac{dv_{i}}{dt}&=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}v_{i}\left(u_{j}-(1+dv_{i})\right)\end{dcases}\quad\forall
i=1,\dots,n\,,$

By using again the definition of $k_{i}^{(in)}$ and by isolating the terms
independent from $j$

$\begin{dcases}\frac{du_{i}}{dt}&=u_{i}\frac{a-u_{i}^{2}}{c}+\frac{u_{i}}{k_{i}^{(in)}}\sum_{j}A_{ij}\left(\frac{u_{j}}{c}-v_{j}\right)\\\
\frac{dv_{i}}{dt}&=-v_{i}(1+dv_{i})+\frac{v_{i}}{k_{i}^{(in)}}\sum_{j}A_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,.$

By adding and removing suitable terms we eventually obtain

$\begin{dcases}\frac{du_{i}}{dt}&=u_{i}\left(\frac{a+bu_{i}-u_{i}^{2}}{c}-v_{i}\right)+\frac{b}{c}u_{i}\sum_{j}\mathcal{L}_{ij}u_{j}-u_{i}\sum_{j}\mathcal{L}_{ij}v_{j}\\\
\frac{dv_{i}}{dt}&=v_{i}\left(u_{i}-(1+dv_{i})\right)+v_{i}\sum_{j}\mathcal{L}_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,.$

The Jacobian matrices follow directly from the definitions of the functions
$\mathbf{f}$ and $\mathbf{F}$:

$\mathbf{J}_{\mathrm{MM}}=\left(\begin{matrix}u_{*}\frac{b-2u_{*}}{c}&-u_{*}\\\
v_{*}&-dv_{*}\end{matrix}\right)\text{ and
}\mathbf{J}_{2}=\left(\begin{matrix}\frac{b}{c}u_{*}&-u_{*}\\\
v_{*}&0\end{matrix}\right)\,,$

and from the latter, one can compute the coefficients $c_{1,i}$,
$i=0,1,2,3,4$, and $c_{2,i}$, $i=0,1,2$, and eventually the polynomials
$S_{1}$ and $S_{2}$. In the present case the resulting expressions are quite
involved and not very explicative, however they can be used to numerically
determine the instability region in the complex plane.

##  Appendix E Analysis of the Volterra model

Let us study a set of $n$ coupled Volterra model McKane and Newman (2005) that
describes the interactions of preys and predators in an ecological setting
where they can interact via long-range connections modelled by

$\begin{dcases}\dfrac{dx_{i}}{dt}=-dx_{i}+ac_{1}y_{i}\frac{1}{k^{(in)}_{i}}\sum_{j}A_{ij}x_{j}+(1-a)c_{1}x_{i}\frac{1}{k^{(in)}_{i}}\sum_{j}A_{ij}y_{j}\\\
\dfrac{dy_{i}}{dt}=ry_{i}-sy_{i}^{2}-c_{2}y_{i}\frac{1}{k^{(in)}_{i}}\sum_{j}A_{ij}x_{j}\,.\end{dcases}$

By using the newly introduced Laplace matrix (3) we can rewrite the previous
equations as:

$\begin{dcases}\dfrac{dx_{i}}{dt}=-dx_{i}+c_{1}y_{i}x_{i}+ac_{1}y_{i}\sum_{j}\mathcal{L}_{ij}x_{j}+(1-a)c_{1}x_{i}\sum_{j}\mathcal{L}_{ij}y_{j}\\\
\dfrac{dy_{i}}{dt}=ry_{i}-sy^{2}_{i}-c_{2}y_{i}x_{i}-c_{2}y_{i}\sum_{j}\mathcal{L}_{ij}x_{j}\,,\end{dcases}$
(34)

where one can easily recognise the in-node Volterra model (15) and the
corrections stemming from non-local contributions. The Jacobian matrices are
obtained as:

$\mathbf{J}_{\mathrm{V}}=\left(\begin{matrix}0&c_{1}x^{*}\\\
-c_{2}y^{*}&-sy^{*}\end{matrix}\right)\text{ and
}\mathbf{J}_{2}=\left(\begin{matrix}ac_{1}y^{*}&(1-a)c_{1}x^{*}\\\
-c_{2}y^{*}&0\end{matrix}\right)\,,$

where $(x^{*},y^{*})$ is the homogeneous equilibrium of the Volterra system.

By inserting the given expressions for $\mathbf{J}_{\mathrm{V}}$ and
$\mathbf{J}_{2}$, in the general formulas (B) and (B) we obtain for the
coefficients of $S_{2}(\xi)$

$\displaystyle c_{2,2}$ $\displaystyle=$
$\displaystyle-(a-1)(c_{1}r-ds)\left(4\frac{(a-1)(c_{1}r-ds)}{c_{1}}+a^{2}d\right)$
$\displaystyle c_{2,1}$ $\displaystyle=$
$\displaystyle(c_{1}r(a-2)+2ds)\left(4\frac{(a-1)(c_{1}r-ds)}{c_{1}}+a^{2}d\right)$
$\displaystyle c_{2,0}$ $\displaystyle=$
$\displaystyle-\frac{(c_{1}r(a-2)+2ds)^{2}}{c_{1}}+\left(\frac{c_{1}r(a-2)+2ds}{c_{1}}\right){ads}+(a-1)(c_{1}r-ds)a^{2}d\,,$

and $S_{1}(\xi)$

$\displaystyle c_{1,4}$ $\displaystyle=$
$\displaystyle-(a-1)(c_{1}r-ds)a^{2}d$ $\displaystyle c_{1,3}$
$\displaystyle=$ $\displaystyle
ad\left(-(c_{1}r(a-2)+2ds)a+2\frac{s(a-1)(c_{1}r-ds)}{c_{1}}\right)$
$\displaystyle c_{1,2}$ $\displaystyle=$
$\displaystyle(c_{1}r-ds)a^{2}d+2\frac{c_{1}r(a-2)+2ds}{c_{1}}ads-\frac{(a-1)(c_{1}r-ds)}{c_{1}^{2}}ds^{2}$
$\displaystyle c_{1,1}$ $\displaystyle=$ $\displaystyle-
ds\left(2\frac{c_{1}r-ds}{c_{1}}a+\frac{c_{1}r(a-2)+2ds}{c_{1}}\frac{s}{c_{1}}\right)$
$\displaystyle c_{1,0}$ $\displaystyle=$
$\displaystyle\frac{c_{1}r-ds}{c_{1}^{2}}ds^{2}\,.$

Given such polynomials one can determine the (in)-stability region as shown in
Fig. 4 (main text) or Fig. 7 and thus conclude about the onset of the
instability according to the position of the complex eigenvalues of the
Laplace matrix $\mathcal{L}$. In Fig. 7 we report the region of instability
(grey) for a set of parameters values allowing for the emergence of patterns
for both the reciprocal and non-reciprocal long-range interactions. We can
observe that, contrary to the case shown in the main text, now the instability
region has a non empty intersection with the real axis where it lies the
spectrum of the Laplace operator for reciprocal interactions (black dots). We
can thus determine a web of symmetrical long-range interactions for which the
Volterra model (15) (main text) exhibits an instability and eventually evolves
toward a patchy solution (panel b) of Fig. 7). A similar result holds true
using non-reciprocal long-range interactions, indeed the complex spectrum of
the associated Laplace matrix (white dots in panel a) of Fig. 7) also lies in
the instability region and thus the system converges to a spatially
heterogeneous solution (panel c) of Fig. 7). The underlying long-range
coupling is a directed Erdős-Rényi network with $n=50$ nodes and a probability
for a direct link to exist between two nodes is $p=0.5$.

Figure 7: Instability region and patterns for the Volterra model. We report
the region of the complex plane $(\Re\Lambda,\Im\Lambda)$ for which the
instability condition is satisfied (grey), the instability is at play if at
least one eigenvalue of $\mathcal{L}$ belongs to the region. The model
parameters have been fixed to the values $c_{1}=2$, $c_{2}=13$, $r=1$, $s=1$,
$d=0.02$ and $a=0.05$, and we can observe that the instability region
intersects the real axis and thus both non-reciprocal (white dots) and
reciprocal (black dots) interactions can exhibit eigenvalues entering into the
instability region. This in turns implies the existence of an heterogeneous
solution for both the reciprocal (see panel b) where we report the density of
preys vs time) and non-reciprocal (see panel c) where we report the density of
preys vs time) long-range interactions assumption. In both panels the
horizontal black line denotes the homogeneous equilibrium $y^{*}$.

In Fig. 8 we provide a more global view of the parameters range associated to
bifurcation diagram showing the parameters values, $d$ and $c_{1}$, for which
the instability emerges in the case of both reciprocal interactions and non-
reciprocal ones (black A region) and in the case of only non-reciprocal ones
(white B region), once the remaining parameters have been fixed to some
generic values. One can clearly appreciate how large is the latter compared to
the former, and thus how more often one can find patterns due to non-
reciprocal interactions instead of reciprocal ones.

Figure 8: Bifurcation diagram for the Volterra model, reciprocal interactions.
We fix the parameters values $c_{2}=13$, $r=1$, $s=1$ and $a=0.05$, and we
show the bifurcation region as a function of the remaining parameters $d$ and
$c_{1}$. The C region (grey) corresponds to an unstable homogeneous
equilibrium that also remain unstable once the coupling is present, being the
latter reciprocal or non-reciprocal one; patterns can develop but they are not
due to the interactions. For parameters values in the A region (black), the
stable homogeneous equilibrium is destabilised by the introduction of
reciprocal coupling (as well by a non-reciprocal one); the patterns shown in
the panel b) of Fig. 7 associated to the values $c_{1}=2$ and $d=0.02$ fall in
this class (yellow dot). Finally in the large B region (white) the homogeneous
equilibrium is stable for any reciprocal coupling, and thus no pattern can
develop in this case, however the use of non-reciprocal interactions can drive
the instability and thus the emergence of patterns. The patterns shown in the
panel b) of Fig. 4 (main text) correspond to the values $c_{1}=2$ and $d=0.4$
clearly belonging to the B region (blue star).

##  Appendix F Analysis of the FitzHugh-Nagumo model

Let us consider the FitzHugh-Nagumo model, $\mathbf{f}(\mathbf{x}^{(i)})=(\mu
u_{i}-u_{i}^{3}-v_{i},\gamma(u_{i}-v_{i}))$, with
$\mathbf{x}^{(i)}=(u_{i},v_{i})$, and the non-local coupling

$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})=(\mu
u_{i}-u_{i}^{3}-v_{j},\gamma(u_{j}-v_{i}))\,,$

with $\mathbf{x}^{(j)}=(u_{j},v_{j})$. From Eq. (2) we obtain

$\begin{dcases}\frac{du_{i}}{dt}&=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}\left(\mu
u_{i}-u_{i}^{3}-v_{j}\right)\\\
\frac{dv_{i}}{dt}&=\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}\gamma\left(u_{j}-v_{i}\right)\end{dcases}\quad\forall
i=1,\dots,n\,.$

By using again the definition of $k_{i}^{(in)}=\sum_{j}A_{ij}$ we get

$\begin{dcases}\frac{du_{i}}{dt}&=\mu
u_{i}-u_{i}^{3}-\frac{1}{k_{i}^{(in)}}\sum_{j}A_{ij}v_{j}\\\
\frac{dv_{i}}{dt}&=-\gamma
v_{i}+\frac{\gamma}{k_{i}^{(in)}}\sum_{j}A_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,,$

hence

$\begin{dcases}\frac{du_{i}}{dt}&=\mu
u_{i}-u_{i}^{3}-v_{i}+\sum_{j}\mathcal{L}_{ij}v_{j}\\\
\frac{dv_{i}}{dt}&=\gamma(u_{i}-v_{i})+\gamma\sum_{j}\mathcal{L}_{ij}u_{j}\end{dcases}\quad\forall
i=1,\dots,n\,.$

The Jacobian matrices can be obtained as follows :

$\mathbf{J}_{\mathrm{FHN}}=\left(\begin{matrix}\mu&-1\\\
\gamma&-\gamma\end{matrix}\right)\text{ and
}\mathbf{J}_{2}=\left(\begin{matrix}0&1\\\ \gamma&0\end{matrix}\right)\,.$

By using the formulas (B) and (B), we obtain the polynomials $S_{1}$ and
$S_{2}$:

$S_{1}(\xi)=\gamma(\mu-\gamma)^{2}\left[(1+\xi)^{2}-\mu\right]\text{ and
}S_{2}(\xi)=-4\gamma^{2}(\xi+1)^{2}-\gamma(\mu-\gamma)^{2}\,.$

The instability is thus realised if the complex eigenvalues
$\Lambda^{(\alpha)}$ of $\mathcal{L}$ do satisfy the constraint

$(\Im\Lambda^{(\alpha)})^{2}\geq\frac{(\mu-\gamma)^{2}}{4\gamma}\frac{(1+\Re\Lambda^{(\alpha)})^{2}-\mu}{(1+\Re\Lambda^{(\alpha)})^{2}+1}\,.$

##  Appendix G Analysis of the Stuart-Landau model

We turn now our attention to the study of the paradigmatic model of nonlinear
oscillators given by the Stuart-Landau system (SL) Stuart and DiPrima (1978);
Kuramoto (1984)

$\frac{dw}{dt}=\sigma w-\beta w|w|^{2}\,,$

where $\sigma=\sigma_{\Re}+i\sigma_{\Im}$ and $\beta=\beta_{\Re}+i\beta_{\Im}$
are complex model parameters. One can straightforwardly prove that
$w_{LC}(t)=\sqrt{\sigma_{\Re}/\beta_{\Re}}e^{i\omega t}$,
$\omega=\sigma_{\Im}-\beta_{\Im}\sigma_{\Re}/\beta_{\Re}$, is a limit cycle
solution of the SL system and it is stable provided $\sigma_{\Re}>0$ and
$\beta_{\Re}>0$.

We then assume to have $n$ identical copies of the SL system coupled through
non-reciprocal interactions, $A_{i\ell}\neq A_{\ell i}$, and the hypothesis of
mean field Eq. (20) (main text), hereby reported

$\frac{dw_{j}}{dt}=\frac{\sigma}{k^{(in)}_{j}}\sum_{\ell}A_{j\ell}w_{\ell}-\beta
w_{j}|w_{j}|^{2}=\sigma w_{j}-\beta
w_{j}|w_{j}|^{2}+\sigma\sum_{\ell}\mathcal{L}_{j\ell}w_{\ell}\,.$

Because of the structure of the coupling, $w_{LC}(t)$ is also a solution of
the latter equation. To inquire about its stability we consider the
perturbation given by
$w_{j}(t)=w_{LC}(t)\left(1+u_{j}(t)\right)e^{iv_{j}(t)}$, where $u_{j}(t)$ and
$v_{j}(t)$ are real and small functions nodes dependent. We then insert the
latter into (20) (main text) and we expand by retaining only the first order
terms by obtaining

$\displaystyle\frac{d}{dt}\binom{u_{j}}{v_{j}}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}-2\sigma_{\Re}&0\\\
-2\beta_{\Im}\frac{\sigma_{\Re}}{\beta_{\Re}}&0\end{matrix}\right)\binom{u_{j}}{v_{j}}+\sum_{\ell}\mathcal{L}_{j\ell}\left(\begin{matrix}\sigma_{\Re}&-\sigma_{\Im}\\\
\sigma_{\Im}&\sigma_{\Re}\end{matrix}\right)\binom{u_{\ell}}{v_{\ell}}\,.$
(35)

We invoke once again the existence of an eigenbasis of the Laplace matrix,
$\phi^{(\alpha)}$, $\Lambda^{(\alpha)}$, to decompose the perturbation $u_{j}$
and $v_{j}$, and eventually get

$\displaystyle\frac{d}{dt}\binom{u^{\alpha}}{v^{\alpha}}$ $\displaystyle=$
$\displaystyle\left[\left(\begin{matrix}-2\sigma_{\Re}&0\\\
-2\beta_{\Im}\frac{\sigma_{\Re}}{\beta_{\Re}}&0\end{matrix}\right)+\Lambda^{(\alpha)}\left(\begin{matrix}\sigma_{\Re}&-\sigma_{\Im}\\\
\sigma_{\Im}&\sigma_{\Re}\end{matrix}\right)\right]\binom{u^{\alpha}}{v^{\alpha}}$
(36) $\displaystyle=$
$\displaystyle\left(\mathbf{J}+\Lambda^{(\alpha)}\mathbf{J}_{2}\right)\binom{u^{\alpha}}{v^{\alpha}}=:\mathbf{J}^{(\alpha)}\binom{u^{\alpha}}{v^{\alpha}}\,.$

By inserting the given expressions for $\mathbf{J}$ and $\mathbf{J}_{2}$, in
the general formulas we obtain for the coefficients of $S_{2}(\xi)$

$\displaystyle c_{2,2}$ $\displaystyle=$
$\displaystyle-\sigma_{\Im}^{2}\left(\sigma_{\Re}^{2}+\sigma_{\Im}^{2}\right)$
$\displaystyle c_{2,1}$ $\displaystyle=$ $\displaystyle
2\sigma_{\Im}^{2}\sigma_{\Re}\left(\beta_{\Re}\sigma_{\Re}+\beta_{\Im}\sigma_{\Im}\right)\frac{1}{\beta_{\Re}}$
$\displaystyle c_{2,0}$ $\displaystyle=$
$\displaystyle-\sigma_{\Im}^{2}\sigma_{\Re}^{2}\left(\beta_{\Re}^{2}+\beta_{\Im}^{2}\right)\frac{1}{\beta_{\Re}^{2}}\,,$
(37)

while for $S_{1}(\xi)$

$\displaystyle c_{1,4}$ $\displaystyle=$
$\displaystyle\sigma_{\Re}^{2}\left(\sigma_{\Re}^{2}+\sigma_{\Im}^{2}\right)$
$\displaystyle c_{1,3}$ $\displaystyle=$
$\displaystyle-2\sigma_{\Re}^{2}\left(2\beta_{\Re}\sigma_{\Re}^{2}+\beta_{\Im}\sigma_{\Im}\sigma_{\Re}+\beta_{\Re}\sigma_{\Im}^{2}\right)\frac{1}{\beta_{\Re}}$
$\displaystyle c_{1,2}$ $\displaystyle=$
$\displaystyle\sigma_{\Re}^{2}\left(5\beta_{\Re}\sigma_{\Re}^{2}+4\beta_{\Im}\sigma_{\Im}\sigma_{\Re}+\beta_{\Re}\sigma_{\Im}^{2}\right)\frac{1}{\beta_{\Re}}$
$\displaystyle c_{1,1}$ $\displaystyle=$
$\displaystyle-2\sigma_{\Re}^{3}\left(\beta_{\Re}\sigma_{\Re}+\beta_{\Im}\sigma_{\Im}\right)\frac{1}{\beta_{\Re}}$
$\displaystyle c_{1,0}$ $\displaystyle=$ $\displaystyle 0\,,$ (38)

As previously done in the case of the Volterra model, the explicit knowledge
of the polynomial $S_{1}$ and $S_{2}$ allows to compute the (in)-stability
region as shown in Fig. 6 (main text) or Fig. 9 in a setting where the
instability condition can be realised for both a reciprocal and non-reciprocal
coupling. The instability region (grey) is shown in the complex plane
$(\Re\Lambda,\Im\Lambda)$ together with the spectrum of a reciprocal web of
long-range interactions (black dots) as well with a non-reciprocal one (white
dots); because in both cases there are eigenvalues belonging to the
instability region, the instability is possible and thus a spatio-temporal
pattern emerges (see panel b) in the case of reciprocal can and panel c) for
the non-reciprocal one). The numerical simulations have been perfumed using a
$4$-th order Runge-Kutta method starting from initial conditions
$\delta$-close to the homogeneous limit cycle solution
$w_{LC}(t)=\sqrt{\sigma_{\Re}/\beta\Re}e^{i\omega t}$. In both cases the
maximum of the dispersion relation $\rho_{\alpha}$ is of order of the unity
and thus a relatively small integration span is sufficient to reveal the wavy
solution. The underlying coupling is a directed Erdős-Rényi network with
$n=40$ nodes and a probability for a direct link to exist between two nodes is
$p=0.08$.

Figure 9: Instability region and patterns for the Stuart-Landau model. In
panel a) we report the region of the complex plane $(\Re\Lambda,\Im\Lambda)$
for which the instability condition is satisfied (grey), for the parameters
values $\sigma_{\Re}=1$, $\sigma_{\Im}=4.3$, $\beta_{\Re}=1$ and
$\beta_{\Im}=-2$. We can observe that the instability region intersects the
real axis and thus any kind of coupling, being reciprocal (black dots) or non-
reciprocal one (white dots), can exhibit eigenvalues entering the instability
region and thus allowing for an instability to set on, followed by a spatio-
temporal patterns as shown in panel b) where we report the real part of the
complex state variable $w_{i}$ in the case of a reciprocal coupling and panel
c) for a non-reciprocal one).

In Fig. 10 we report the bifurcation diagram in the plane $\sigma_{\Im}$ and
$\beta_{\Im}$, for $\sigma_{\Re}=\beta_{\Re}=1$. Two regions can be observed;
in region A (black) the instability can be initiated by both a reciprocal and
non-reciprocal web of long-range interactions while in region B (white) only
non-reciprocal interactions can determine an instability and the ensuing wavy
solution.

Figure 10: Bifurcation diagram for the Stuart-Landau model. For the set of
parameters values $\sigma_{\Re}=1.0$ and $\beta_{\Re}=1.0$ we report the range
of the remaining parameters $\sigma_{\Im}$ and $\beta_{\Im}$ for which the
instability can emerge. In the region A (black) the homogeneous equilibrium is
stable and it can be destabilised by a reciprocal coupling as well by a non-
reciprocal one. On the other hand, in the B region (white), the homogeneous
equilibrium is always stable, no pattern can develop using a symmetric
coupling. On the contrary one can found non-reciprocal coupling capable to
destabilise the homogeneous equilibrium and thus the system to develop a wavy
heterogeneous solution. The patterns shown in the panel b) of Fig. 6 (main
text) correspond to the values $\sigma_{\Im}=4$ and $\beta_{\Im}=1$ and
clearly belong to the B region (blue star). The patterns presented in Fig. 9
for parameters values $\sigma_{\Im}=4$ and $\beta_{\Im}=-2$ correspond to the
region A (yellow dot).

##  Appendix H Instability onset in the general $d$-dimensional case

Let us consider $n$ copies of a generic $d$-dimensional system coupled through
non-reciprocal interactions as in Eq. (4), hereby recalled

$\frac{d\mathbf{x}^{(i)}}{dt}=\mathbf{f}(\mathbf{x}^{(i)})+\sum_{j}\mathcal{L}_{ij}\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(j)})\quad\forall
i=1,\dots,n\,.$

Assume there exists an homogeneous solution $\mathbf{x}^{(i)}=\mathbf{s}$ for
the decoupled system, that clearly results also a solution for the coupled
one. To study its stability in the coupled case, let us linearise the system
about this reference solution to obtain (5), that can be further analysed by
projecting it onto the Laplace eigenbasis to eventually get (7).

The stability of the reference solution can be determined by studying the sign
of the real part of the roots of the characteristic polynomial (21)

$p_{d}(\lambda):=(-1)^{d}\left[a_{0}\lambda^{d}+a_{1}\lambda^{d-1}+a_{2}\lambda^{d-2}+\dots+a_{d}\right]\,.$

A straightforward computation, ultimately based on the expansion of the
determinant defining the characteristic polynomial, allows to determine

$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle-\mathrm{tr}\left(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2}\right)-i\Im\Lambda^{(\alpha)}\mathrm{tr}\mathbf{J}_{2}$
(39) $\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{d}\sum_{j=\ell+1}^{d}\Big{\\{}(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{\ell\ell}(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{jj}-(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{\ell
j}(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{j\ell}+$ (40)
$\displaystyle+$
$\displaystyle(\Im\Lambda^{(\alpha)})^{2}\left[(\mathbf{J}_{2})_{j\ell}(\mathbf{J}_{2})_{\ell
j}-(\mathbf{J}_{2})_{\ell\ell}(\mathbf{J}_{2})_{jj}\right]\Big{\\}}+$
$\displaystyle+$ $\displaystyle
i\Im\Lambda^{(\alpha)}\sum_{\ell=1}^{d}\sum_{j=\ell+1}^{d}\Big{\\{}(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{\ell\ell}(\mathbf{J}_{2})_{jj}+(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{jj}(\mathbf{J}_{2})_{\ell\ell}-(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{\ell
j}(\mathbf{J}_{2})_{j\ell}+$ $\displaystyle-$
$\displaystyle(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2})_{j\ell}(\mathbf{J}_{2})_{\ell
j}\Big{\\}}\,.$

Assuming to deal with reciprocal interactions, we can show (see Appendix A)
that the spectrum is real even if the Laplace matrix is not symmetric. The
previous equations simplify by imposing $\Im\Lambda^{(\alpha)}=0$, and in
particular $a_{1}$ and $a_{2}$ are real numbers, that are assumed to be
positive because of the hypothesis of stable equilibrium under symmetric
interactions.

Let us now discuss the role of the auxiliary polynomial $q(\lambda)$. Let
$P(z)$ by a generic polynomial with complex coefficients. Assume $z_{1}$ to be
a simple complex root and to simplify let us assume to be able to factorise
$P(z)=(z-z_{1})P_{1}(z)$ where $P_{1}(z)$ has real coefficients. Let us define
$Q(z)=(z-\overline{z_{1}})P(z)$, then

$Q(z)=(z-\overline{z_{1}})P(z)=(z-\overline{z_{1}})(z-z_{1})P_{1}(z)=(z^{2}-2z\Re
z_{1}+|z_{1}|^{2})P_{1}(z)\,,$

from which it follows that $Q(z)$ has real coefficients and its roots have the
same real part of the roots of $P(z)$ (computed twice in the case of $z_{1}$).
So the general strategy would be to factorise all the complex roots of $P(z)$
and build an auxiliary polynomial $Q(z)$ according to this recipe. In general
we do not know all such roots and thus we can simply assume to apply this
process to all roots, even for the real ones, hence to the whole polynomial,
and thus to define $Q(z)=\overline{P(\overline{z})}P(z)$. The resulting
polynomial will have a degree equals to the double of the degree of $P(z)$,
all its coefficients will be real by construction and its roots will have the
same real part of the roots of $P(z)$.

By applying such recipe to
$p_{d}(\lambda)=\lambda^{d}+a_{1}\lambda^{d-1}+a_{2}\lambda^{d-2}+\dots+a_{d}$
we obtain the polynomial
$q(\lambda)=\lambda^{2d}+b_{1}\lambda^{2d-1}+b_{2}\lambda^{2d-2}+\dots+b_{2d}$
defined in the main text. A direct computation consisting in equating
coefficients of equal powers of $\lambda$ in $q(z)$ and
$\overline{p_{d}(\overline{z})}p_{d}(z)$ allows to determine

$\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle
a_{1}+\overline{a_{1}}=2\Re a_{1}\,,$ $\displaystyle b_{2}$ $\displaystyle=$
$\displaystyle a_{2}+\overline{a_{2}}+a_{1}\overline{a_{1}}=2\Re
a_{2}+|a_{1}|^{2}\,.$

By using the previously obtained expression (39) for $a_{1}$, we can get

$b_{1}=-2\mathrm{tr}\left(\mathbf{J}_{1}+\Re\Lambda^{(\alpha)}\mathbf{J}_{2}\right)\,,$

and we can observe that it coincides (module a factor $2$) with $a_{1}$
obtained once we impose $\Im\Lambda^{(\alpha)}=0$; hence $b_{1}$ results to be
always positive, as $a_{1}$ does. The non-reciprocal interactions cannot thus
change the sign of $b_{1}$. On the other hand by using the expression for
$a_{2}$ given in (40), we can obtain the conditions presented in the main text
(22) and (23) ensuring $b_{2}<0$. Indeed recalling (40) we can write

$b_{2}=2\Re a_{2}+|a_{1}|^{2}=\left(2\Re
a_{2}+|a_{1}|^{2}\right)_{\Im\Lambda^{(\alpha)}=0}+2(\Im\Lambda^{(\alpha)})^{2}\sum_{\ell=1}^{d}\sum_{j=\ell+1}^{d}\Big{\\{}\left[(\mathbf{J}_{2})_{j\ell}(\mathbf{J}_{2})_{\ell
j}-(\mathbf{J}_{2})_{\ell\ell}(\mathbf{J}_{2})_{jj}\right]\Big{\\}}+(\Im\Lambda^{(\alpha)})^{2}\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2}\,.$

By observing that $\left(2\Re
a_{2}+|a_{1}|^{2}\right)_{\Im\Lambda^{(\alpha)}=0}>0$ because of the
assumption of stability of the system using symmetric interaction, we can
conclude that $b_{2}$ can be negative if $(\Im\Lambda^{(\alpha)})^{2}$ is
large enough and

$2\sum_{\ell=1}^{d}\sum_{j=\ell+1}^{d}\Big{\\{}\left[(\mathbf{J}_{2})_{j\ell}(\mathbf{J}_{2})_{\ell
j}-(\mathbf{J}_{2})_{\ell\ell}(\mathbf{J}_{2})_{jj}\right]\Big{\\}}+\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2}<0\,.$

After some algebraic manipulations, we can rewrite the term involving the
double sum as

$\sum_{\ell=1}^{d}\sum_{j=\ell+1}^{d}\Big{\\{}\left[(\mathbf{J}_{2})_{j\ell}(\mathbf{J}_{2})_{\ell
j}-(\mathbf{J}_{2})_{\ell\ell}(\mathbf{J}_{2})_{jj}\right]\Big{\\}}=\frac{\mathrm{tr}\mathbf{J}^{2}_{2}-\left(\mathrm{tr}\mathbf{J}_{2}\right)^{2}}{2}\,,$

from which the condition (22) follows. This is a sufficient conditions for the
destabilisation of the equilibrium solution; if the latter does not hold true
one could in principle look for other conditions capable to change the sign of
one of the remaining coefficients, $b_{j}$, of the polynomial $q(\lambda)$.