# Multiple colour interacting urns on complete graphs

Benito Pires  and  Rafael A. Rosales  Departamento de Computação e Matemática,
Faculdade de Filosofia, Ciências e Letras, Universidade de São Paulo,
14040-901, Ribeirão Preto - SP, Brazil<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We present a multiple colour generalisation of the model of graph interacting
urns studied by Benaïm _et. al._ , Random Struct. Alg., 46: 614-634, 2015. We
show that for complete graphs and for a broad class of reinforcement functions
governing the addition of balls in the urns, the process of colour proportions
at each urn converges almost surely to the fixed points of the reinforcement
function.

###### Key words and phrases:

generalized Pólya urns, stochastic approximations, reinforced processes

###### 2020 Mathematics Subject Classification:

Primary: 60K35, 37C10; Secondary: 62L20, 37A50

The first author was partially supported by grant #2019/10269-3 and
#2022/14130-2 São Paulo Research Foundation (FAPESP)

##  1\. Introduction

Interacting Pólya urn models have recently gained considerable attention, see
for instance [1], [5], [8], [9], [10] and [12]. In this article we consider a
multiple colour generalisation of the single-colour model introduced in [3],
in which the interaction depends on the geometry of an underlying graph. To
describe the model in [3], let $\rscal{G}=(V,E)$ be a finite, connected,
undirected graph such that $V=[d]=\\{1,\ldots,d\\}$ and at each vertex $u\in
V$ an urn with $B_{u}(0)\geq 1$ balls is placed at time $n=0$. At each time
$n\geq 1$ and to each edge $\\{u,v\\}\in E$, one ball is added either to $u$
or to $v$, and the probability that the ball is added to $u$ is

(1) $\mathbb{P}\big{(}u\,\,\textrm{is chosen among}\,\,\\{u,v\\}\,\,\textrm{at
time}\,\,n+1\big{)}=\dfrac{B_{u}(n)^{\alpha}}{B_{u}(n)^{\alpha}+B_{v}(n)^{\alpha}},$

where $B_{u}(n)$ denotes the number of balls in the urn $u$ at time $n\geq 0$,
and $\alpha>0$ is a parameter that reinforces the effect of the number of
balls at each vertex. The almost sure asymptotic behaviour of the proportions

(2) $X_{u}(n)=\dfrac{B_{u}(n)}{\sum_{v\in V}B_{v}(n)}$

of balls in each urn at each time is studied in [3]. In such a study, the
geometry of $\rscal{G}$ and the value of $\alpha$ play a fundamental role. A
complete description for finite, connected graphs in the linear case, i.e.
when $\alpha=1$, is given in [4] and [11]. Notice that the transition
probability in (1) expressed in terms of (2) takes the form
$\psi\big{(}X_{v}(n)/X_{u}(n)\big{)}$, where
$\psi(t)=1/\big{(}1+t^{\alpha}\big{)}$.

To introduce our model, suppose that $\rscal{G}=(V,E)$ is a finite,
undirected, complete graph with at least $d\geq 2$ vertices. Assume that an
urn with $B^{i}_{u}(0)\geq 1$ balls of colour $i\in C=\\{1,2,\ldots,c\\}$,
$c\geq 2$, is placed at each vertex $u\in V$. To simplify matters, we assume
that there exists $b_{0}\geq 1$ such that $\sum_{u\in V}B_{u}^{i}(0)=b_{0}$
for all $i\in C$, that is, we start the process with the same total number of
balls of each colour. The number of balls of colour $i$ in the urn $u\in V$ at
time $n\geq 0$ is denoted by $B_{u}^{i}(n)$. Let

(3) $X_{u}^{i}(n)=\frac{B_{u}^{i}(n)}{\sum_{v\in V}B_{v}^{i}(n)}$

be the proportion of balls of colour $i$ at time $n$ in the urn $u$. At each
time $n$ and to each edge $\\{u,v\\}\in E$, a ball of each colour is added
either to $u$ or to $v$, and the probability that the ball of colour $i$ is
added to $u$ is

(4) $\mathbb{P}\bigg{(}\textrm{\parbox{170.71652pt}{\centering the ball of
colour $i$ at edge $\\{u,v\\}$\\\ goes to vertex $u$ at step $n+1$
\@add@centering} }\\!\\!\bigg{)}=\frac{e^{\sum_{k\in
C{\setminus}\\{i\\}}\alpha X_{u}^{k}(n)}}{e^{\sum_{k\in
C{\setminus}\\{i\\}}\alpha X_{u}^{k}(n)}+e^{\sum_{k\in
C{\setminus}\\{i\\}}\alpha X_{v}^{k}(n)}}.$

When $\alpha<0$, the probability that a ball of colour $i$ is put into the urn
$u$ decreases exponentially with the proportion of balls of all other colours
in that urn. This model can therefore be used to describe a competition among
the different colours for the urns. The case $\alpha>0$, on the other hand,
defines a model for cooperative dynamics. The reinforcement parameter $\alpha$
strengthens the competition or the cooperation among different colours. The
model may be used to describe the dynamics of residential segregation of
individuals from $C$ groups into $d$ regions identified with the vertices of
$\rscal{G}$. It may also be conceived as repeated population game among agents
of $C$ types or the competition of $d$ companies for $C$ products; providing a
generalisation of the motivating examples of the single colour model in [3].

Setting $t=\sum_{k\in C{\setminus}\\{i\\}}\big{(}\alpha X_{u}^{k}(n)-\alpha
X_{v}^{k}(n)\big{)}$, shows that the transition probability in (4) is given by
a reinforcement function of the form $\psi(t)=1/\big{(}1+e^{-t}\big{)}$. This
is merely an example of a broader family of functions considered throughout.
The results described in Section 2 are valid for a broader family of smooth
functions $\varphi:\mathbb{R}\to(0,1)$ satisfying the following hypotheses:

  1. ($h_{1}$)

there exists $M>0$ such that $0<\varphi^{\prime}(t)<M$ for all
$t\in\mathbb{R}$;

  2. ($h_{2}$)

$\varphi(t)+\varphi(-t)=1$ for all $t\in\mathbb{R}$.

The class of functions satisfying these hypotheses includes among others the
following functions:

$\varphi(t)=\dfrac{1}{1+e^{-t}},\quad\varphi(t)=\dfrac{1}{2}(1+\tanh
t),\quad\varphi(t)=\dfrac{1}{2}\left(1+\dfrac{t}{\sqrt{t^{2}+1}}\right).$

In these terms, the nonlinear generalised Pólya urn considered throughout is
set by defining

(5) $\mathbb{P}{\bigg{(}\,\textrm{\parbox{170.71652pt}{\centering the ball of
colour $i$ at edge $\\{u,v\\}$\\\ goes to vertex $u$ at step $n+1$
\@add@centering}\\! }\bigg{)}}=\varphi\bigg{(}\sum\nolimits_{k\in
C{\setminus}\\{i\\}}\left(\alpha X_{u}^{k}(n)-\alpha
X_{v}^{k}(n)\right)\\!\bigg{)}.$

The main results of this article are described in Section 2. Their proofs are
presented in Sections 3, 4 and 5.

##  2\. Statement of main results

By using (3), define the process of colour proportions as

(6) $X=\\{X(n)\\}_{n\geq 0},$

where

$X(n)=\big{(}X^{1}_{1}(n),\ldots,X_{d}^{1}(n)\ldots,X_{1}^{c}(n),\ldots,X_{d}^{c}(n)\big{)}.$

Observe that for each $n$, $X(n)$ takes values on

$\Delta^{c}=\bigg{\\{}(x_{u}^{i})\in\mathbb{R}^{dc}:x_{u}^{i}\geq 0,\forall
u,i,\,\,\textrm{and}\,\,{\sum_{u=1}^{d}}x_{u}^{i}=1,\forall i\bigg{\\}}.$

We assume hereafter that the process $X$ is defined on an appropriate
probability space $(\Omega,\mathcal{F},\mathbb{P})$.

We are principally concerned with a strong law for $X$. To state our main
result, let $\pi={(\pi_{u}^{i}):}\Delta^{\\!c}\to\Delta^{\\!c}$ be the map
whose coordinate-function at ${x}=\big{(}x_{u}^{i}\big{)}\in\Delta^{\\!c}$
takes the value

(7) $\pi_{u}^{i}({x})=\dfrac{1}{|E|}\sum_{w\in
V{\setminus}\\{u\\}}\varphi\left(\sum_{k\in C{\setminus}\\{i\\}}\left(\alpha
x_{u}^{k}-\alpha x_{w}^{k}\right)\right).$

It follows from the hypotheses ($h_{1}$) and ($h_{2}$) that $\pi$ takes
$\Delta^{\\!c}$ into the interior of $\Delta^{c}$, as $\pi_{u}^{i}>0$. Denote
by Fix$(\pi)$ the set of fixed points of $\pi$, i.e.,
Fix$(\pi)=\\{x\in\mathbb{R}^{dc}:\pi(x)=x\\}$.

######  Theorem 1.

Let $\rscal{G}=(V,E)$ be a finite complete graph, $X=\\{X(n)\\}_{n\geq 0}$ be
the process of colour proportions at each urn given by (6), and $\pi$ be the
map defined in (7). If the set of fixed points of $\pi$ is finite, then
$X=\\{X(n)\\}_{n\geq 0}$ converges almost surely to the fixed points of $\pi$,
that is,

$\sum_{x\in\
{\textup{Fix}}(\pi)}\mathbb{P}\left(\lim_{n\to\infty}X(n)=x\right)=1.$

It is important to remark that Theorem 1 is a global result, demanding strong
hypotheses. More precisely, the graph needs to be complete and the
reinforcement function needs to be biologically or physically inspired so that
there exists a strict Lyapunov function defined on the entire phase space.
This is exactly the case of this article since the function $\varphi$ in the
right-hand-side of (5) satisfies the hypotheses $(h1)$ and $(h2)$, which
appear in a biologically inspired neural network (Hopfield network) that
models associative memories in live beings.

As the next result clarifies, the hypothesis that $\pi$ has finitely many
fixed points in Theorem 1 is automatic when the parameter $\alpha$ in (7) has
small absolute value. In this case, there exists a unique asymptotic state,
i.e., a unique point $x\in\Delta^{c}$ such that $\lim_{n\to\infty}X(n)=x$. On
the other hand, it is an open question whether for all $\alpha$ satisfying
$|\alpha|\geq
d/\big{(}4(c-1)\sup_{t\in\mathbb{R}}|\varphi^{\prime}(t)|\big{)}$, the number
of fixed points of $\pi$ is finite.

######  Corollary 1.

Let $\rscal{G}=(V,E)$ be a finite complete graph, $X=\\{X(n)\\}_{n\geq 0}$ be
the process of colour proportions at each urn given by (6), and $\pi$ be the
map defined in (7). If
$|\alpha|<d\big{/}\big{(}4(c-1)\sup_{t\in\mathbb{R}}|\varphi^{\prime}(t)|\big{)}$,
then $X=\\{X(n)\\}_{n\geq 0}$ converges almost surely to the unique fixed
point $x^{*}=\left(\frac{1}{d},\ldots,\frac{1}{d}\right)\in\mathbb{R}^{dc}$ of
$\pi$.

In what follows, we denote by
$\rho(A)=\max\,\\{|\lambda|:\lambda\,\,\textrm{is eigenvalue of $A$}\\}$ the
spectral radius of a square matrix $A$. Given $x\in\Delta^{c}$, we let
$J\pi(x)$ denote the Jacobian matrix of the map $\pi$ at $x$. We say that a
fixed-point $x$ of $\pi$ is stable if $\rho\big{(}J\pi(x)\big{)}<1$. If
$\alpha$ has large absolute value, then $\pi$ may have multiple stable fixed
points, as the following example shows.

######  Example 1.

Let $V=\\{1,2,3\\}$, $C=\\{1,2\\}$, $\alpha=\frac{615}{182}\ln 9$,
$\varphi(t)=1/(1+e^{-t})$, $t\in\mathbb{R}$, and $\pi$ be the map defined in
(7). Then $\textrm{Fix}\,(\pi)$ contains the stable fixed points:

(8)
$\displaystyle\textstyle\left(\frac{23}{615},\frac{1}{3},\frac{129}{205},\frac{23}{615},\frac{1}{3},\frac{129}{205}\right),$
$\displaystyle\textstyle\left(\frac{23}{615},\frac{129}{205},\frac{1}{3},\frac{23}{615},\frac{129}{205},\frac{1}{3}\right),$
$\displaystyle\textstyle\left(\frac{1}{3},\frac{23}{615},\frac{129}{205},\frac{1}{3},\frac{23}{615},\frac{129}{205}\right)$
$\displaystyle\textstyle\left(\frac{129}{205},\frac{1}{3},\frac{23}{615},\frac{129}{205},\frac{1}{3},\frac{23}{615}\right),$
$\displaystyle\textstyle\left(\frac{129}{205},\frac{23}{615},\frac{1}{3},\frac{129}{205},\frac{23}{615},\frac{1}{3}\right),$
$\displaystyle\textstyle\left(\frac{1}{3},\frac{129}{205},\frac{23}{615},\frac{1}{3},\frac{129}{205},\frac{23}{615}\right)$

Figure 1 presents six sample paths for a single colour of the process
$X=\\{X(n)\\}_{n\geq 0}$ started at
$X(0)=\left(\frac{1}{3},\frac{1}{3},\frac{1}{3},\frac{1}{3},\frac{1}{3},\frac{1}{3}\right)$.
The figure includes simulations for two different reinforcement strengths,
$\alpha=\frac{3}{4}$ and $\alpha=({615}/{182})\ln 9$. Open circles are located
at the fixed points of $\pi$. In the later case, the six non-central circles
are located at the points with coordinates determined by (8). Continuous light
curves represent the flow of the ODE in (9). These simulations are consistent
with the general result stated in Theorem 1 and also with Corollary 1.

|
---|---
$\alpha=\frac{3}{4}$ | $\alpha=\frac{615}{182}\ln 9$
Figure 1. Sample paths for a single colour and different reinforcement
parameters $\alpha$.

##  3\. Proof of Theorem 1

The almost sure convergence of $X$ is established by using the ODE method for
stochastic approximations (see Appendix A). This technique shows that the
limit set of the process $X$ is almost surely connected and contained in the
chain-recurrent set (see Definition 2) of the ordinary differential equation

(9) $\dot{x}=F(x):=-x+\pi(x),\quad x\in\Delta^{\\!c}.$

If $\pi$ has finitely many fixed points, then the chain-recurrent set of (9)
equals the set of fixed points of $\pi$ (see Corollary 2). To achieve this
result, we construct a strict Lyapunov function for (9) in terms of energy
functions of Hopfield’s Networks.

We will need some notions from the theory of ordinary differential equations.
The function $F$ in (9) is defined on the entire $\mathbb{R}^{dc}$. By
($h_{1}$) and ($h_{2}$), $F$ is a smooth map with bounded partial derivatives.
Under such conditions, there exists a semi-flow associated to (9) and denoted
by $\phi_{t}:\mathbb{R}^{dc}\to\mathbb{R}^{dc}$, $t\geq 0$. A set
$S\subset\mathbb{R}^{dc}$ is called invariant by (9) if
$\phi_{t}\big{(}{x}\big{)}\in S$ for all $t\geq 0$ and ${x}\in S$. A point
${x}\in\mathbb{R}^{dc}$ is an equilibrium point if $\phi_{t}({x})={x}$ for all
$t\geq 0$. We denote by $\boldsymbol{E}$ the set of equilibrium points of (9).

######  Definition 1 (strict Lyapunov function).

A continuous function $L:S\to\mathbb{R}$ is a _strict Lyapunov function_ for
(9) on an invariant set $S\subset\mathbb{R}^{dc}$ if the function
$t\in[0,\infty)\mapsto L\big{(}\phi_{t}({x})\big{)}$ is strictly decreasing
for all ${x}\in S{\setminus}\boldsymbol{E}$.

######  Theorem 2.

The following statements concerning the differential equation (9) hold:

  * $(i)$

$\Delta^{\\!c}\subset\mathbb{R}^{dc}$ is an invariant set;

  * $(\mathit{ii})$

There exists a strict Lyapunov function $L:\Delta^{\\!c}\to\mathbb{R}$ on
$\Delta^{\\!c}$.

###### Proof of Theorem 2.$(i)$.

Let ${x}=(x_{u}^{i})\in\Delta^{\\!c}$, that is, assume that $x_{u}^{i}\geq 0$
for all $u\in V$, $i\in C$, and $\sum_{u\in V}x_{u}^{i}=1$ for all $i\in C$.
By (7) and ($h_{2}$), since $\rscal{G}=(V,E)$ is an undirected complete graph,
we have that

$\displaystyle\sum_{u\in
V}\pi_{u}^{i}({x})=\dfrac{1}{{|E|}}\sum_{\begin{subarray}{c}u,w\in V\\\
u<w\end{subarray}}\Bigg{[}$ $\displaystyle\varphi\Bigg{(}\sum_{k\in
C{\setminus}\\{i\\}}\Big{(}\alpha x_{u}^{k}-\alpha x_{w}^{k}\Big{)}\Bigg{)}$
$\displaystyle+\ \ $ $\displaystyle\varphi\Bigg{(}\sum_{k\in
C{\setminus}\\{i\\}}\Big{(}\alpha x_{w}^{k}-\alpha
x_{u}^{k}\Big{)}\Bigg{)}\Bigg{]}={\dfrac{1}{{|E|}}\,{|E|}=1}.$

In this way, $\pi({x})$ takes $\Delta^{\\!c}$ into its interior because
$\pi_{u}^{i}({x})>0$ for all $u,i$, and $\sum_{u\in V}\pi_{u}^{i}({x})=1$ for
all $i\in C$. In particular, if ${x}$ is a point in the boundary of
$\Delta^{\\!c}$, we have that $F({x})=-{x}+\pi({x})$ is the displacement
vector from ${x}\in\partial{\Delta^{\\!c}}$ to the point
$\pi\big{(}{{x}}\big{)}$ in the interior of $\Delta^{\\!c}$. Hence, by the
convexity of $\Delta^{\\!c}$, we have that $F({x})$ points towards the
interior of $\Delta^{\\!c}$. As a result, $\Delta^{\\!c}$ is invariant by (9).
∎

###### Proof of Theorem 2.$(\mathit{ii})$.

Introducing the new variables $x_{uv}^{i}=\alpha x_{u}^{i}-\alpha x_{v}^{i}$,
$i\in C,\,u,v\in V$, $u\neq v$, in (9) leads to the system of $c\cdot{|E|}$
differential equations

(10) $\dot{x}_{uv}^{i}=-x_{uv}^{i}+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}\varphi\left(\sum_{{k\in
C{\setminus}\\{i\\}}}x_{uw}^{k}\right)-\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}\varphi\left(\sum_{k\in
C{\setminus}\\{i\\}}x_{vw}^{k}\right).$

By using ($h_{2}$), we may rewrite (10) in the form

(11) $\displaystyle\dot{x}_{uv}^{i}=$ $\displaystyle-
x_{uv}^{i}+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}\varphi\left(\sum_{{k\in
C{\setminus}\\{i\\}}}x_{uw}^{k}\right)+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}\varphi\left(\sum_{k\in
C{\setminus}\\{i\\}}x_{wv}^{k}\right)$
$\displaystyle-\dfrac{\alpha}{{|E|}}\sum_{w\in V{\setminus}\\{v\\}}1.$

Setting ${x}_{uv}=\left[x_{uv}^{1},x_{uv}^{2},\ldots,x_{uv}^{c}\right]^{t}$,
$u\neq v$, and rewriting (11) in vectorial form leads to

(12) $\displaystyle\dot{{x}}_{uv}=$
$\displaystyle-{x}_{uv}+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}{\Phi}(B{x}_{uw})+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}{\Phi}(B{x}_{wv})$
$\displaystyle-\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}\mathbf{1},\qquad u\neq v,$

where $\mathbf{1}=[1,1,\ldots,1]^{t}$, $B=(b_{ij})$ is the non-singular
$c\times c$ matrix with entries $b_{ij}=1-\delta_{ij}$, where $\delta_{ij}=1$
if $i=j$ and $\delta_{ij}=0$ otherwise, and
${\Phi}:\mathbb{R}^{c}\to\mathbb{R}^{c}$ is the map defined by
$\Phi([p_{1},p_{2},\ldots,p_{c}]^{t})=[\varphi(p_{1}),\varphi(p_{2}),\ldots,\varphi(p_{c})]^{t}$.
Introducing the new variables
${y}_{uv}=\left[y_{uv}^{1},y_{uv}^{2},\ldots,y_{uv}^{c}\right]^{t}=B{x}_{uv}$,
$u\neq v$, in (12) yields

(13) $\dot{{y}}_{uv}=-{y}_{uv}+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}B{\Phi}({y}_{uw})+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}B{\Phi}({y}_{wv})-\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}B\mathbf{1}.$

After rewriting (13) in the scalar form, we obtain

(14) $\displaystyle\dot{y}_{uv}^{i}=$ $\displaystyle-
y_{uv}^{i}+\dfrac{\alpha}{{|E|}}\sum_{w\in V{\setminus}\\{u\\}}\sum_{k\in
C{\setminus}\\{i\\}}\varphi\big{(}y_{uw}^{k}\big{)}+\dfrac{\alpha}{{|E|}}\sum_{w\in
V{\setminus}\\{v\\}}\sum_{k\in
C{\setminus}\\{i\\}}\varphi\big{(}y_{wv}^{k}\big{)}$
$\displaystyle-\dfrac{\alpha}{{|E|}}(d-1)(c-1).$

Again, let us introduce the new variables

(15) $z_{iuv}=y_{uv}^{i},\quad i\in C,\quad u,v\in V,\quad u\neq v.$

Then we may rewrite (14) in the form

(16)
$\dot{z}_{iuv}=-z_{iuv}+\sum_{(j,r,s)\in\Lambda}w_{iuv}^{jrs}\,\varphi(z_{jrs})+m$

where

(17) $\begin{gathered}\Lambda=\\{(i,r,s):i\in C,\,r,s\in V,\,r\neq
s\\},\\\\[7.5pt] m=-\dfrac{\alpha}{{|E|}}(d-1)(c-1),\quad
w_{iuv}^{jrs}=\begin{cases}\dfrac{2\alpha}{{|E|}}&\textrm{if}\quad j\neq
i,\,r=u,\,s=v\\\\[7.22743pt] \dfrac{\alpha}{{|E|}}&\textrm{if}\quad j\neq
i,\,r=u,\,s\neq v\\\\[7.22743pt] \dfrac{\alpha}{{|E|}}&\textrm{if}\quad j\neq
i,\,r\neq u,\,s=v\\\\[7.22743pt] \phantom{1}0&\textrm{if}\quad j\neq i,\,r\neq
u,\,s\neq v\\\\[7.22743pt] \phantom{1}0&\textrm{if}\quad
j=i\end{cases}\end{gathered}$

It is clear that the multidimensional matrix $W=\big{(}w_{iuv}^{jrs}\big{)}$
is symmetric in the sense that $w_{iuv}^{jrs}=w_{jrs}^{iuv}$ for all
$(i,u,v)\in\Lambda$, $(j,r,s)\in\Lambda$.

The system of equations (16) is analogous to the one governing the dynamics of
a Hopfield network. Provided that $W=\big{(}w_{iuv}^{jrs}\big{)}$ is
symmetric, [7, p. 2090, Equation (7)] proved (see Appendix B) that the map
$L:\mathbb{R}^{c|E|}\to\mathbb{R}$ defined by

$\displaystyle L({z})=-\frac{1}{2}\sum_{(i,u,v)\in\Lambda}$
$\displaystyle\sum_{(j,r,s)\in\Lambda}w_{iuv}^{jrs}\,\varphi\big{(}z_{iuv}\big{)}\varphi\big{(}z_{jrs}\big{)}$
$\displaystyle+\sum_{(i,u,v)\in\Lambda}\int^{\varphi(z_{iuv})}\varphi^{-1}(\eta)\mathrm{d}\eta+\sum_{(i,u,v)\in\Lambda}m\varphi(z_{iuv}),$

is strict Lyapunov function for (16). In this way, if $z:t\in[0,\infty)\mapsto
z(t)\in\mathbb{R}^{c|E|}$ is a non-constant solution of (16), then the
function $t\in[0,\infty)\mapsto L\big{(}z(t)\big{)}$ is strictly decreasing.
Now we will explain how to translate such a result about the solutions of (16)
into a result about the solutions of (9). Let
$\Psi:\mathbb{R}^{dc}\to\mathbb{R}^{c|E|}$ be the map that at
$x=(x_{u}^{i})\in\mathbb{R}^{dc}$ takes the value
$z=(z_{iuv})\in\mathbb{R}^{c|E|}$, where

(18) $z=\Psi(x),\quad z_{iuv}=\sum_{k\in C{\setminus}\\{i\\}}\big{(}\alpha
x_{u}^{k}-\alpha x_{v}^{k}\big{)},\quad i\in C,\quad u,v\in V,\quad u\neq v.$

It follows from all the change of variables made within this section that if
$x:t\in[0,\infty)\mapsto x(t)\in\mathbb{R}^{dc}$ is a solution of (16) then
$z:[0,\infty)\to\mathbb{R}^{c|E|}$ defined by
$z(t)=\Psi\big{(}x(t)\big{)},t\in\mathbb{R}$, is a solution of (16). To show
that $L\circ\Psi$ is a Lyapunov function for (9) on $\Delta^{c}$, we will show
that if $x:t\in[0,\infty)\mapsto x(t)\in\Delta^{c}$ is a non-constant solution
of (9) on $\Delta^{c}$, then $t\in[0,\infty)\mapsto
L\Big{(}\Psi\big{(}x(t)\big{)}\Big{)}$ is strictly decreasing. To proceed with
that reasoning, let $x:t\in[0,\infty)\mapsto x(t)\in\Delta^{c}$ be a non-
constant solution of $\eqref{eqn:THE_field}$ and let $z:t\in[0,\infty)\mapsto
z(t)=\Psi\big{(}x(t)\big{)}$ be the corresponding solution of (16). The proof
will be finished if we show that $z$ is non-constant because in that case, as
already remarked, the function $t\in[0,\infty)\mapsto
L\big{(}z(t)\big{)}=L\Big{(}\Psi\big{(}x(t)\big{)}\Big{)}$ is strictly
decreasing. By way of contradiction, assume that $z$ is constant. By (18) and
(7), there exist $b_{u}^{i}\in\mathbb{R}$, $i\in C$, $u\in V$, such that

(19) $b_{u}^{i}=\pi_{u}^{i}\big{(}x(t)\big{)},\quad t\geq 0,\quad i\in C,\quad
u\in V.$

In this way, $b:=(b_{u}^{i})=\pi\big{(}x(t)\big{)}\in\Delta^{c}$. In
particular, we have that

(20) $\sum_{u\in V}b_{u}^{i}=1,\quad\forall i\in C.$

By (9) and (19), we obtain that $x:t\in[0,\infty)\mapsto
x(t)\in\mathbb{R}^{dc}$ satisfies the uncoupled system of differential
equations:

$\dot{x}_{u}^{i}(t)=-x_{u}^{i}(t)+b_{u}^{i},\quad t\geq 0,\quad i\in C,\quad
u\in V.$

Applying integrating factor to the previous ode, we conclude that there exist
$a_{u}^{i}\in\mathbb{R}$, $i\in C$, $u\in V$ such that

(21) $x_{u}^{i}(t)=a_{u}^{i}e^{-t}+b_{u}^{i},\quad t\geq 0,\quad i\in C,\quad
u\in V.$

Now given, $u,v\in V$ with $u\neq v$, let
${x}_{uv}(t)=\left[x_{uv}^{1}(t),x_{uv}^{2}(t),\ldots,x_{uv}^{c}(t)\right]^{t}$,
where $x_{uv}^{i}(t)=\alpha x_{u}^{i}(t)-x_{v}^{i}(t)$. Since, by hypothesis,
$z$ is constant, we have that $\dot{z}(t)=0$ for all $t\geq 0$. Hence, by
(18), we have that

$0=\dot{z}_{iuv}(t)=\sum_{k\in
C{\setminus}\\{i\\}}\big{(}\alpha\dot{x}_{u}^{k}(t)-\alpha\dot{x}_{v}^{k}(t)\big{)}=\left[B\dot{x}_{uv}(t)\right]_{i},\quad
t\geq 0,\quad i\in C.$

That is, $B\dot{x}_{uv}(t)=0$ for all $t\geq 0$. Since $B$ is non-singular, we
have that $\dot{x}_{uv}(t)=0$ for all $t\geq 0$ or, equivalently,

(22) $\dot{x}_{u}^{i}(t)=\dot{x}_{v}^{i}(t),\quad t\geq 0,\quad i\in C,\quad
u,v\in V,\quad u\neq v.$

Derivating (21) and using (22) imply that that there exist $a^{i}$, $i\in C$,
such that

$x_{u}^{i}(t)=a^{i}e^{-t}+b_{u}^{i},\quad t\geq 0,\quad i\in C,\quad u\in V.$

Summing over $u\in V$ and using (20) yield

$1=a^{i}e^{-t}+1,\quad\forall t\geq 0,\quad\forall i\in C,$

which is a contradiction. In fact, since $x$ is a non-constant solution, we
have that $a^{i}\neq 0$ for at least one $i\in C$. ∎

######  Definition 2 (chain-recurrent set).

Let $\delta>0$ and $T>0$. A $(\delta,T)$-pseudo orbit of (9) from
$x\in\Delta^{\\!c}$ to $y\in\Delta^{\\!c}$ is a finite sequence of partial
orbits $\\{\phi_{t}(y_{i}):0\leq t\leq t_{i}\\}$; $i=0,\ldots,k-1$; $t_{i}\geq
T$ of (9) such that

$\|y_{0}-x\|<\delta,\quad\|\phi_{t_{i}}(y_{i})-y_{i+1}\|<\delta,\,\,i=0,\ldots,k-1,\quad\textrm{and}\quad
y_{k}=y.$

A point $x\in\Delta^{\\!c}$ is chain-recurrent if for every $\delta>0$ and
$T>0$ there is a $(\delta,T)$-pseudo orbit from $x$ to itself. The chain-
recurrent set of (9), denoted by $\mathcal{R}$, is the set of all chain-
recurrent points.

######  Corollary 2.

The chain-recurrent set $\mathcal{R}$ of (9) equals the set of fixed points of
the map $\pi$.

###### Proof.

The set of fixed points of $\pi$, $\textup{Fix}(\pi)$, equals the set
$\boldsymbol{E}$ of equilibrium points of (9). It is an elementary fact that
$\boldsymbol{E}\subset\mathcal{R}$. By hypothesis, $\textup{Fix}(\pi)$, and
therefore $\boldsymbol{E}$, is finite. Hence, if $L$ is the strict Lyapunov
function in Theorem 2, then $L(\boldsymbol{E})$ is finite. In this case, by
Proposition 3.2 in [2], it follows that $\mathcal{R}\subset\boldsymbol{E}$. We
have proved that $\mathcal{R}=\boldsymbol{E}=\textrm{Fix}\,(\pi)$. ∎

######  Definition 3 (limit set).

Let $X=\\{X(n)\\}_{n\geq 0}$ be the process given by (6), which is defined on
the probability space $(\Omega,\mathcal{F},\mathbb{P})$. Given
$\omega\in\Omega$, we say that $x\in\Delta^{\\!c}$ is a limit point of
$\\{X(n,\omega)\\}_{n\geq 0}$ if there exists a strictly increasing sequence
of integers $1\leq n_{1}<n_{2}<\cdots$ such that
$\lim_{k\to\infty}X(n_{k},\omega)=x$.

###### Proof of Theorem 1.

By the ODE method for stochastic approximation (see Appendix A), there exists
a full measure set $\Omega^{\prime}\subset\Omega$ such that for each
$\omega\in\Omega^{\prime}$, the limit set of
$\left\\{X(n,\omega)\right\\}_{n\geq 0}$ is a connected subset of the set
$\mathcal{R}$ of chain-recurrent points of (9). By Corollary 2, the connected
components of $\mathcal{R}$ consists of the fixed points of $\pi$, which are
finitely many. Hence, the limit set of $\left\\{X(n,\omega)\right\\}_{n\geq
0}$ is a set consisting of a fixed point of $\pi$. Since the limit set is a
unique point, we conclude that $\lim_{n\to\infty}X(n,\omega)$ exists and
belongs to $\textrm{Fix}\,(\pi)$. ∎

##  4\. Proof of Corollary 1

To prove Corollary 1, we will show that if $|\alpha|$ is enough small, then
the map $\pi$ defined in (7) is a contraction on $\Delta^{\\!c}$ with respect
to the norm $\|x\|_{1}=\sum_{(u,i)\in V\times C}|x_{u}^{i}|$, that is, $\pi$
satisfies the following conditions:

  * (i)

$\pi$ takes the compact, convex set $\Delta^{\\!c}$ into itself;

  * (ii)

there exists $0<\lambda<1$ such that
$\|\pi(x)-\pi(y)\|_{1}\leq\lambda\|x-y\|_{1}$, $\forall x,y\in\Delta^{\\!c}$.

Condition (i) was already verified, see the proof of Theorem 2. To show that
Condition (ii) is true, let $J\pi(x)=\big{(}(\partial\pi_{u}^{i}/\partial
x_{v}^{j})(x)\big{)}$ denote the Jacobian matrix of $\pi$ at
$x\in\Delta^{\\!c}$. Since $\Delta^{\\!c}$ is compact and convex, by the Mean
Value Theorem, $\pi$ will be a contraction on $\Delta^{\\!c}$ provided
$\sup_{x\in\Delta^{\\!c}}\left\|J\pi(x)\right\|_{1}<1$ or equivalently, if the
following condition is true

  * (iii)

$\displaystyle\sum_{(u,i)\in V\times
C}\left|\dfrac{\partial\pi_{u}^{i}}{\partial x_{v}^{j}}(x)\right|<1$,
$\quad\forall(v,j)\in V\times C$, $\forall x\in\Delta^{\\!c}$.

In fact, for all $(v,j)\in V\times C$ and $x\in\Delta^{\\!c}$, we have that

$\dfrac{\partial\pi_{u}^{i}}{\partial
x_{v}^{j}}(x)=\begin{cases}\phantom{-}0&\textrm{if}\quad i=j\\\
\phantom{-}\sum_{w\in
V{\setminus}\\{u\\}}\dfrac{\alpha}{|E|}\varphi^{\prime}\left(\sum_{k\in
C{\setminus}\\{i\\}}\left(\alpha x_{u}^{k}-\alpha
x_{w}^{k}\right)\right)&\textrm{if}\quad i\neq j,\,\,u=v\\\
-\dfrac{\alpha}{|E|}\varphi^{\prime}\left(\sum_{k\in
C{\setminus}\\{i\\}}\left(\alpha x_{u}^{k}-\alpha
x_{v}^{k}\right)\right)&\textrm{if}\quad i\neq j,\,\,u\neq v\end{cases}.$

In this way, for all $(v,j)\in V\times C$ and $x\in\Delta^{\\!c}$, since
${|E|}=d(d-1)/2$, we have that

$\displaystyle\sum_{(u,i)\in V\times
C}\left|\dfrac{\partial\pi_{u}^{i}}{\partial x_{v}^{j}}(x)\right|$
$\displaystyle\leq
2(c-1)(d-1)\dfrac{|\alpha|}{|E|}\sup_{t\in\mathbb{R}}\left|\varphi^{\prime}(t)\right|$
$\displaystyle=4(c-1)\dfrac{|\alpha|}{d}\sup_{t\in\mathbb{R}}\left|\varphi^{\prime}(t)\right|<1.$

Hence, Condition (iii) holds and $\pi$ is a contraction on $\Delta^{\\!c}$.
Now since $\Delta^{\\!c}$ is a compact (therefore complete) subset of
$\mathbb{R}^{dc}$, by the Banach Fixed Point Theorem, we have that $\pi$ has a
unique fixed point. We claim that such a fixed point is
$x^{*}=\left(\frac{1}{d},\ldots,\frac{1}{d}\right)$. In fact, by $(h2)$, we
have that $\varphi(0)=\frac{1}{2}$. By (7), it follows that

$\pi_{u}^{i}(x^{*})=\dfrac{1}{|E|}\sum_{w\in
V{\setminus}\\{u\\}}\varphi\left(\sum_{k\in
C{\setminus}\\{i\\}}\left(\dfrac{\alpha}{d}-\dfrac{\alpha}{d}\right)\right)=\dfrac{2}{d(d-1)}(d-1)\varphi(0)=\dfrac{1}{d}=(x^{*})_{u}^{i}.$

To conclude the proof, apply Theorem 1.

##  5\. Proof of Example 1

First we provide conditions under which a point
$x=(x_{1}^{1},x_{2}^{1},x_{3}^{1},x_{1}^{2},x_{2}^{2},x_{3}^{2})\in\Delta^{2}$
is a fixed-point of $\pi$. We write $x=(y,z)$, where

$y=(y_{1},y_{2},y_{3})=(x_{1}^{1},x_{2}^{1},x_{3}^{1})\quad\textrm{and}\quad
z=(z_{1},z_{2},z_{3})=(x_{1}^{2},x_{2}^{2},x_{3}^{2}).$

By (7), we have that $x=(y,z)$ is a fixed-point of $\pi$ if and only if the
following system of equations is satisfied

(23)
$\begin{cases}y_{u}=\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
z_{u}-\alpha z_{w}\big{)}\\\\[7.22743pt]
z_{u}=\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
y_{u}-\alpha y_{w}\big{)}\end{cases},\quad u\in\\{1,2,3\\}.$

######  Lemma 1.

Let $x=(y,z)\in\Delta^{2}$ be such that $z=y$ and

(24)
$y_{u}=\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
z_{u}-\alpha z_{w}\big{)},\quad u\in\\{1,2,3\\}.$

Then $x$ is a fixed-point of $\pi$.

###### Proof.

By hypothesis the first three equations in (23) are satisfied. To prove that
the fourth one is true, let $u=1$. Then, we have that

$\displaystyle z_{1}=y_{1}$ $\displaystyle=$
$\displaystyle\dfrac{1}{3}\varphi\big{(}\alpha z_{1}-\alpha
z_{2}\big{)}+\dfrac{1}{3}\varphi\big{(}\alpha z_{1}-\alpha
z_{3}\big{)}=\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{1\\}}\varphi\big{(}\alpha
y_{1}-\alpha y_{w}\big{)},$

which shows that the fourth equation in (23) holds. Proceeding in the same way
for $u=2$ and $u=3$, we obtain that the fifth and the sixth equations in (23)
are true, implying that $x=(y,z)=(y,y)$ is a fixed-point of $\pi$. ∎

In what follows, given a permutation $\sigma$ on the set $V=\\{1,2,3\\}$ and a
vector $w=(w_{1},w_{2},w_{3})$, we use the notation
$\omega^{\sigma}=\big{(}\omega_{\sigma(1)},\omega_{\sigma(2)},\omega_{\sigma(3)}\big{)}$.

######  Lemma 2.

If $x=(y,z)\in\Delta^{2}$ satisfies (24), then for any permutation $\sigma$ on
the set $V=\\{1,2,3\\}$, we have that $x^{*}=(y^{\sigma},z^{\sigma})$ also
satisfies (24).

###### Proof.

Suppose that $x=(y,z)$ satisfies $\eqref{yu}$ and let $\sigma$ be a
permutation on the set $V=\\{1,2,3\\}$. Let $x^{*}=(y^{*},z^{*})$, where
$y^{*}=y^{\sigma}$ e $z^{*}=z^{\sigma}$. Then, we have that

$\displaystyle y_{u}^{*}=y_{\sigma(u)}$ $\displaystyle=$
$\displaystyle\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{\sigma(u)\\}}\varphi\big{(}\alpha
z_{\sigma(u)}-\alpha z_{w}\big{)}$ $\displaystyle=$
$\displaystyle\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
z_{\sigma(u)}-\alpha z_{\sigma(w)}\big{)}$ $\displaystyle=$
$\displaystyle\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
z_{u}^{*}-\alpha z_{w}^{*}\big{)}$

In this way, $x^{*}=(y^{*},z^{*})=(y^{\sigma},z^{\sigma})$ satisfies (24). ∎

###### Proof of Example 1.

First we consider the point

$x=(y,z)=\left(\dfrac{23}{615},\dfrac{1}{3},\dfrac{129}{205},\dfrac{23}{615},\dfrac{1}{3},\dfrac{129}{205}\right).$

To show that $x$ is a fixed-point of $\pi$, we will verify the hypotheses of
Lemma 1. First notice that $z=y$. Moreover, it is elementary to show that if
$\alpha=\frac{615}{182}\ln\,9$, then the system of equations in (1) is
verified. Hence, by Lemma 1, $x$ is a fixed-point of $\pi$. Now let $x^{*}$ be
any of the other $5$ points in (8). As
$x^{*}=(y^{*},z^{*})=(y^{\sigma},z^{\sigma})$, from Lemmas 1 and 2 it follows
that $x^{*}$ is a fixed-point of $\pi$.

Now it remains to prove that at each of the $6$ fixed-points, each eigenvalue
of $J\pi(x)$ has absolute value less than $1$. By (7), we have that at each
any point $x=(y,z)$ of $\Delta^{2}$,

$\pi_{u}^{i}(x)=\pi_{u}^{i}(y,z)=\begin{cases}\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
z_{u}-\alpha z_{w}\big{)}&\textrm{if}\quad i=1\\\\[7.22743pt]
\dfrac{1}{3}\sum_{w\in\\{1,2,3\\}{\setminus}\\{u\\}}\varphi\big{(}\alpha
y_{u}-\alpha y_{w}\big{)}&\textrm{if}\quad i=2\end{cases},\quad
u\in\\{1,2,3\\}.$

In particular, the Jacobian matrix of $\pi$ at $x=(y,y)$ is given by

(25) $J\pi(x)=J\pi(y,y)=\begin{pmatrix}0&M(y)\\\ M(y)&0\end{pmatrix},$

where $0$ denotes the null matrix of order $3$,

$M(y)=\begin{pmatrix}\frac{\alpha}{3}\varphi^{\prime}\big{(}y_{12}\big{)}+\frac{\alpha}{3}\varphi^{\prime}\big{(}y_{13}\big{)}&-\frac{\alpha}{3}\varphi^{\prime}(y_{12})&-\frac{\alpha}{3}\varphi^{\prime}(y_{13})\\\\[5.0pt]
-\frac{\alpha}{3}\varphi^{\prime}(y_{21})&\frac{\alpha}{3}\varphi^{\prime}\big{(}y_{21}\big{)}+\frac{\alpha}{3}\varphi^{\prime}\big{(}y_{23}\big{)}&-\frac{\alpha}{3}\varphi^{\prime}(y_{23})\\\\[5.0pt]
-\frac{\alpha}{3}\varphi^{\prime}(y_{31})&-\frac{\alpha}{3}\varphi^{\prime}(y_{32})&\frac{\alpha}{3}\varphi^{\prime}(y_{31})+\frac{\alpha}{3}\varphi^{\prime}(y_{32})\end{pmatrix},$

and $y_{uw}=\alpha y_{u}-\alpha y_{w}$ for all $u,w\in\\{1,2,3\\}$ with $u\neq
w$.

Since $\varphi(t)=1/(1+e^{-t})$, by $(h2)$,
$\varphi^{\prime}(t)=(1-\varphi(t))\varphi(t)=\varphi^{\prime}(-t)$.
Therefore, since $y_{uw}=-y_{wu}$, we have that $M(y)$ is a symmetric matrix.
Now we will do the calculation for the point

$x=(y,y)=\left(\dfrac{23}{615},\dfrac{1}{3},\dfrac{129}{205},\dfrac{23}{615},\dfrac{1}{3},\dfrac{129}{205}\right).$

Since $\alpha=\frac{615}{182}\ln 9$, we have that

$y_{12}=-\ln 9,\quad y_{13}=-2\ln 9,\quad y_{23}=-\ln 9.$

Therefore, since $\varphi^{\prime}(t)=(1-\varphi(t))\varphi(t)$, we have that

(26)
$M(y)=\frac{\alpha}{3}\begin{pmatrix}\phantom{-}\frac{8577}{84050}&-\frac{9}{100}&-\frac{81}{6724}\\\\[5.0pt]
-\frac{9}{100}&\phantom{-}\frac{9}{50}&-\frac{9}{100}\\\\[5.0pt]
-\frac{81}{6724}&-\frac{9}{100}&\phantom{-}\frac{8577}{84050}\end{pmatrix}$

We will need the following lemma.

######  Lemma 3.

If $\lambda$ is a real eigenvalue of $J\pi(x)$, then $\lambda^{2}$ is a real
eigenvalue of $\big{(}M(y)\big{)}^{2}$.

###### Proof.

Let $\lambda$ be an eigenvalue of $J\pi(x)$ associated to an eigenvector

$v=\begin{pmatrix}p\\\ q\end{pmatrix},\quad\textrm{where}\quad
p=\begin{pmatrix}p_{1}\\\ p_{2}\\\ p_{3}\end{pmatrix}\quad\textrm{and}\quad
q=\begin{pmatrix}q_{1}\\\ q_{2}\\\ q_{3}\end{pmatrix}.$

Then, by (25), we have that

$\begin{pmatrix}\lambda p\\\ \lambda q\end{pmatrix}=\lambda
v=J\pi(x)v=\begin{pmatrix}0&M(y)\\\ M(y)&0\end{pmatrix}\begin{pmatrix}p\\\
q\end{pmatrix}=\begin{pmatrix}M(y)q\\\ M(y)p\end{pmatrix}.$

In particular, we have that

$\big{(}M(y)\big{)}^{2}p=M(y)M(y)p=M(y)\lambda q=\lambda M(y)q=\lambda^{2}p,$

showing that $\lambda^{2}$ is eigenvalue of $\big{(}M(y)\big{)}^{2}$. ∎

By (26), the eigenvalues of $M(y)$ are approximately $0$, $0.28$ and $0.66$.
Therefore, the eigenvalues of $\big{(}M(y)\big{)}^{2}$ are real and less than
$1$. By (25), $J\pi(x)$ is symmetric, therefore the eigenvalues of $J\pi(x)$
are real. By Lemma 3, the eigenvalues of $J\pi(x)$ are real numbers with
absolute value less than $1$. More precisely, if $\lambda$ is eigenvalue of
$J\pi(x)$, then $\lambda$ is real and $\lambda^{2}\in[0,1)$. This implies that
$\lambda\in(-1,1)$. Hence, the spectral radius of $J\pi(x)$ is less than $1$.
The same procedure can be applied to the other $5$ fixed-points of $\pi$. ∎

##  Appendix A Stochastic approximations

In what follows, denote by $(\Omega,\mathcal{F},\mathbb{P})$ the probability
space where the process $X=\\{X(n)\\}_{n\geq 0}$ given by (6) is defined. Set
$\mathcal{F}_{0}=\\{\emptyset,\Omega\\}$ and for each $n\geq 1$, let
$\mathcal{F}_{n}$ be the $\sigma$-algebra generated by the process $X$ from
time 1 up to time $n$. We denote by $F:\mathbb{R}^{dc}\to\mathbb{R}^{dc}$ the
vector field defined in (9).

######  Lemma 4.

The process $X=\\{X(n)\\}_{n\geq 0}$ defined by (6) satisfies the recursion

(27) $X(n+1)-X(n)={\gamma_{n}}\big{(}F(X(n))+{U(n)}\big{)},\quad n\geq 0,$

where

  1. $(i)$

$\\{\gamma_{n}\\}_{n\geq 0}$ is a deterministic sequence of positive numbers;

  2. $(ii)$

$\sum_{n\geq 0}\gamma_{n}=\infty$ and $\sum_{n\geq 0}\gamma_{n}^{2}<\infty$;

  3. $(iii)$

$U(n):\Omega\to\mathbb{R}^{dc}$ is a $\mathcal{F}_{n}$-measurable map for each
$n\geq 0$;

  4. $(iv)$

$\mathbb{E}[{U(n)}\mid\mathcal{F}_{n}]=0$ for all $n\geq 0$;

  5. $(v)$

For each $T>0$, almost surely it holds that

$\lim_{n\to\infty}\Bigg{(}\sup_{\\{\,r\,:\,0\,\leq\,\tau_{r}-\tau_{n}\,\leq\,T\,\\}}\Bigg{\|}\sum_{k=n}^{r-1}\gamma_{k}U(k)\Bigg{\|}\Bigg{)}=0,$

where $\tau_{0}=0$ and $\tau_{n}=\sum_{k=0}^{n-1}\gamma_{k}$.

###### Proof of $(i)$-$(iv)$.

For any $\\{u,v\\}\in E$, $i\in C$, and $n\geq 0$, let $B_{v\to u}^{i}(n+1)$
be the event

$B_{v\to u}^{i}(n+1)=\Big{\\{}\,\text{the ball of colour $i$ at edge
$\\{u,v\\}$ goes to vertex $u$ at $n+1$}\,\Big{\\}}.$

Define

(28) $N^{i}_{u}(n+1)=\sum_{w\in V{\setminus}\\{u\\}}\mathbf{1}_{B^{i}_{w\to
u}(n+1)}$

that is, $N^{i}_{u}(n+1)$ is the number of balls of colour $i$ put into the
urn at vertex $u$ at time $n+1$. We recall that the number of balls of colour
$i$ at urn $u$ at time $n$ is denoted by $B_{u}^{i}(n)$. Moreover, as an
hypothesis of the model, we assumed that there exists $b_{0}\geq 1$ such that
$\sum_{u\in V}B_{u}^{i}(0)=b_{0}$ for all $i\in C$. In this way, for each
$i\in C$ and $n\geq 0$,

$\sum_{u\in V}B_{u}^{i}(n)=b_{0}+n{|E|}.$

By (3), the increment in the proportion of balls of colour $i$ at vertex $u$
at time $n+1$ equals

(29) $\displaystyle X^{i}_{u}(n+1)-X^{i}_{u}(n)$
$\displaystyle=\frac{B^{i}_{u}(n)+N^{i}_{u}(n+1)}{b_{0}+(n+1){|E|}}-\frac{B^{i}_{u}(n)}{b_{0}+n{|E|}}$
$\displaystyle=\frac{-{|E|}X^{i}_{u}(n)+N^{i}_{u}(n+1)}{b_{0}+(n+1){|E|}}$
$\displaystyle=\frac{1}{b_{0}/{|E|}+n+1}\Big{(}-X^{i}_{u}(n)+\frac{1}{{|E|}}N^{i}_{u}(n+1)\Big{)}.$

Set

(30)
${\gamma_{n}}=1/(b_{0}/{|E|}+n+1)\quad\textrm{and}\quad\xi^{i}_{u}({n})=\frac{1}{{|E|}}N^{i}_{u}(n+1).$

For each ${n\geq 1}$, define $\xi(n)=\big{(}\xi^{1}_{1}(n)$, $\ldots$,
$\xi^{1}_{d}(n)$, $\ldots$ $\xi^{c}_{1}(n)$, $\ldots$,
$\xi^{c}_{d}(n)\big{)}$. By considering the expression for the increments of
each component of $X$ given by (29), it follows that for any $n\geq 0$,

(31) $X(n+1)-X(n)=\gamma_{{n}}\big{(}-X(n)+\xi({n})\big{)}.$

According to (30) and (28), for any $i\in C$,

(32)
$\displaystyle\begin{split}\mathbb{E}\big{[}\xi^{i}_{u}({n})\mid\mathcal{F}_{n}\big{]}&=\frac{1}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}\mathbb{E}\big{[}\mathbf{1}_{B^{i}_{w\to
u}(n+1)}\mid\mathcal{F}_{n}\big{]}\\\ &=\frac{1}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}\mathbb{P}\big{(}B^{i}_{w\to
u}(n+1)\mid\mathcal{F}_{n}\big{)}.\end{split}$

Markov’s property implies that $\mathbb{P}\big{(}B^{i}_{w\to
u}(n+1)\mid\mathcal{F}_{n}\big{)}=\mathbb{P}\big{(}B^{i}_{w\to
u}(n+1)\mid{X(n)}\big{)}$. Hence, by (32), (5), (7), we obtain

(33)
$\mathbb{E}\big{[}\xi^{i}_{u}({n})\mid\mathcal{F}_{n}\big{]}=\frac{1}{{|E|}}\sum_{w\in
V{\setminus}\\{u\\}}\mathbb{P}\big{(}B^{i}_{w\to
u}(n+1)\mid{X(n)}\big{)}=\pi_{u}^{i}\big{(}X(n)\big{)}.$

By (33) and (9),

(34)
$-X(n)=F\big{(}X(n)\big{)}-\pi\big{(}X(n)\big{)}=F\big{(}X(n)\big{)}-\mathbb{E}\big{[}\xi({n+1})\mid\mathcal{F}_{n}\big{]}.$

Replacing $-X(n)$ in the right hand side of (31) by (34) yields

$X(n+1)-X(n)=\gamma_{{n}}\left(F\big{(}X(n)\big{)}+U({n})\right),$

where $U(n)=\xi(n)-\mathbb{E}[\xi(n)\mid\mathcal{F}_{n}]$.

We have proved that the increments of $X$ satisfy (27). Moreover, it follows
from the expressions of $\gamma_{n}$ and $U(n)$ that $(i)$-$(iv)$ hold.

Proof of $(v)$: Let $M_{n}=\sum_{k=0}^{n}\gamma_{k}U(k)$. An immediate
verification shows that $\\{M_{n};n\geq 0\\}$ is a martingale adapted to the
filtration $\\{\mathcal{F}_{n+1}:n\geq 0\\}$, that is
$\mathbb{E}[M_{n+1}\mid\mathcal{F}_{n+1}]=M_{n}$. It then follows that almost
surely

$\displaystyle{\sum_{k=0}^{n}}\mathbb{E}\big{[}\|M_{k+1}-M_{k}\|^{2}\big{|}\
{\mathcal{F}_{k+1}}\big{]}$
$\displaystyle={\sum_{k=0}^{n}}\gamma_{k+1}^{2}\mathbb{E}\big{[}\|{U(k+1)}\|^{2}\big{|}\
{\mathcal{F}_{k+1}}\big{]}$
$\displaystyle\leq{\sum_{k=0}^{n}}\gamma_{k+1}^{2}\bigg{(}\sum_{u\in
V}\sum_{i\in C}{\xi^{i}_{u}(k+1)}\bigg{)}^{2}$
$\displaystyle\leq\bigg{(}cd\frac{(d-1)}{{|E|}}\bigg{)}^{2}{\sum_{k=0}^{n}}\gamma^{2}_{k+1}$
$\displaystyle\leq 4c^{2}\sum_{k\geq 0}\gamma^{2}_{k+1}<\infty.$

This implies that $M_{n}$ converges almost surely to a finite limit, see
Theorem 5.4.9 in [6], and hence that $\\{M_{n}:n\geq 0\\}$ is a Cauchy
sequence. This concludes the proof of assertion $(v)$. ∎

In what follows, we will need Definitions 2 and 3.

######  Theorem 3.

Let $X=\\{X(n)\\}_{n\geq 0}$ be the process given by (6), which is defined on
the probability space $(\Omega,\mathcal{F},\mathbb{P})$. There exists a full
measure set $\Omega^{\prime}\subset\Omega$ such that for each
$\omega\in\Omega^{\prime}$, the limit set of $\\{X(n,\omega)\\}_{n\geq 0}$ is
a connected set contained in the chain-recurrent set of (9).

###### Proof.

This is a well established result for any process satisfying the list of
conditions in the statement of Lemma 4, see for instance [2]. ∎

##  Appendix B Hopfield’s energy function

######  Theorem 4 (Hopfield).

Let $\Lambda$ be a finite set and
$G=(G_{\mu})_{\mu\in\Lambda}:\mathbb{R}^{|\Lambda|}\to\mathbb{R}^{|\Lambda|}$
be the vector field whose coordinate-function at $z=(z_{\mu})$ takes the value

(35)
$G_{\mu}(z)=-z_{\mu}+\sum_{\nu\in\Lambda}w_{\mu\nu}\varphi(z_{\nu})+I_{\mu},\quad\mu\in\Lambda,$

where $W=(w_{\mu\nu})$ is a symmetric real matrix, $I_{\mu},\mu\in\Lambda$,
are constants and $\varphi:\mathbb{R}\to(0,1)$ is any function satisfying the
hypothesis ($h_{1}$). Then the function
$L:\mathbb{R}^{|\Lambda|}\to\mathbb{R}$ defined by

$L(z)=-\frac{1}{2}\sum_{\mu,\nu\in\Lambda}w_{\mu\nu}\,\varphi(z_{\mu})\varphi(z_{\nu})+\sum_{\mu\in\Lambda}\int^{\varphi(z_{\mu})}\varphi^{-1}(\eta)\mathrm{d}\eta+\sum_{\mu\in\Lambda}I_{\mu}\varphi(z_{\mu}),$

is a strict Lyapunov function for the ordinary differential equation

(36) $\dot{z}=G(z)$

on $\mathbb{R}^{|\Lambda|}$ and it is called Hopfield’s energy function.

###### Proof.

The proof we present is adapted from [7], where it appears in a different
context and with another notation. By direct derivation, since
$W=(w_{\mu\nu})$ is symmetric, we have that for all $t\geq 0$,

(37)
$\begin{split}\dfrac{d}{dt}L\big{(}z(t)\big{)}=\sum_{\mu\in\Lambda}\Bigg{(}-\sum_{\nu\in\Lambda}w_{\mu\nu}\varphi(z_{\nu}(t))&\varphi^{\prime}(z_{\mu}(t))+z_{\mu}(t)\varphi^{\prime}(z_{\mu}(t))\\\
&+I_{\mu}\varphi^{\prime}(z_{\mu}(t))\Bigg{)}G_{\mu}\big{(}z(t)\big{)},\end{split}$

where $z(t)$ denotes the solution of (36) that satisfies
$z(0)=z^{(0)}\in\mathbb{R}^{|\Lambda|}$. By using (35), we may rewrite (37) in
the form

$\dfrac{d}{dt}L\big{(}z(t)\big{)}=-\sum_{\mu\in\Lambda}\varphi^{\prime}\big{(}z_{\mu}(t)\big{)}\left[G_{\mu}\big{(}z(t)\big{)}\right]^{2}.$

If $z^{(0)}$ is not an equilibrium point of (36), then
$\left[G_{\mu}\big{(}z(t)\big{)}\right]^{2}>0$ for all $t\geq 0$. Moreover, by
the hypothesis ($h_{1}$), we have $\varphi^{\prime}\big{(}z_{\mu}(t)\big{)}>0$
for all $t\geq 0$. This shows that $L$ is a strict Lyapunov function for (36)
on $\mathbb{R}^{|\lambda|}$. The other condition in hypothesis ($h_{1}$)
regarding the boundedness of $\varphi^{\prime}$ is important to ensure that
the solutions of (36) are defined for all $t\geq 0$. ∎

To use Theorem 4 in our article, consider the notation used in (16) and set

$\displaystyle\Lambda=\\{(i,u,v):i\in C,u,v\in V,u\neq v\\},$
$\displaystyle\mu=(i,u,v)\in\Lambda,\quad\nu=(j,r,s)\in\Lambda,\quad
I_{\mu}=m,\quad w_{\mu\nu}=w_{iuv}^{jrs}.$

Acknowledgements. We acknowledge Paulo José Rodriguez for writing the code of
the simulations shown in Figure 1.

## References

  * [1] Giacomo Aletti and Andrea Ghiglietti, _Interacting generalized friedman’s urn systems_ , Stochastic Processes and their Applications 127 (2017), no. 8, 2650–2678.
  * [2] Michel Benaïm, _A dynamical system approach to stochastic approximations_ , SIAM J. Control Optim. 34 (1996), no. 2, 437–472. MR 1377706
  * [3] Michel Benaïm, Itai Benjamini, Jun Chen, and Yuri Lima, _A generalized Pólya’s urn with graph based interactions_ , Random Structures Algorithms 46 (2015), no. 4, 614–634. MR 3346459
  * [4] Jun Chen and Cyrille Lucas, _A generalized Pólya’s urn with graph based interactions: convergence at linearity_ , Electron. Commun. Probab. 19 (2014), no. 67, 13.
  * [5] Marcelo Costa and Jonathan Jordan, _Phase transitions in non-linear urns with interacting types_ , Bernoulli 28 (2022), no. 4, 2546–2562. MR 4474553
  * [6] Rick Durrett, _Probability: theory and examples_ , fourth ed., Cambridge Series in Statistical and Probabilistic Mathematics, vol. 31, Cambridge University Press, Cambridge, Cambridge, UK, 2010. MR 2722836
  * [7] John J. Hopfield, _Neurons with graded response have collective computational properties like those of 2-state neuros_ , Proceedings of the National Academy of Sciences of the United States of America 81 (1984), no. 10, 3088–3092.
  * [8] Gursharn Kaur and Neeraja Sahasrabudhe, _Interacting urns on a finite directed graph_ , Journal of Applied Probability (2022), 1–23.
  * [9] Sophie Laruelle and Gilles Pagès, _Nonlinear randomized urn models: a stochastic approximation viewpoint_ , Electron. J. Probab. 24 (2019), Paper No. 98, 47. MR 4017116
  * [10] Mickaël Launay and Vlada Limic, _Generalized interacting urn models_ , 2012, https://arxiv.org/abs/1207.5635.
  * [11] Yuri Lima, _Graph-based Pólya’s urn: completion of the linear case_ , Stoch. Dyn. 16 (2016), no. 2, 1660007, 13. MR 3470556
  * [12] Wioletta M. Ruszel and Debleena Thacker, _Positive reinforced generalized time-dependent pólya urns via stochastic approximation_ , 2022, https://arxiv.org/abs/2201.12603.


# Ballistic magnetic thermal transport coupled to phonons

P. Zavitsanos1 and X. Zotos1,2,3 1Department of Physics, University of Crete,
70013 Heraklion, Greece  2Foundation for Research and Technology - Hellas,
71110 Heraklion, Greece  3Leibniz Institute for Solid State and Materials
Research IFW Dresden, 01171 Dresden, Germany

###### Abstract

Motivated by thermal conductivity experiments in spin chain compounds, we
propose a phenomenological model to account for a ballistic magnetic transport
coupled to a diffusive phononic one, along the line of the seminal two-
temperature diffusive transport Sanders-Walton model. Although the expression
for the effective thermal conductivity is identical to that of Sanders-Walton,
the interpretation is entirely different, as the ”magnetic conductivity” is
replaced by an ”effective transfer conductivity” between the magnetic and
phononic component. This model also reveals the fascinating possibility of
visualizing the ballistic character of magnetic transport, for appropriately
chosen material parameters, as a two peak counter-propagating feature in the
phononic temperature. It is also appropriate for the analysis of any thermal
transport experiment involving a diffusive component coupled to a ballistic
one.

The thermal transport by magnetic excitations has been extensively studied
over the last few years hess  and it was established as a novel thermal
conduction mechanism besides the well known electronic and phononic ones. In
particular, it was studied in one dimensional quantum magnets ballistic  where
it was shown that, in compounds accurately described by the Heisenberg
spin-1/2 chain model, there is a ballistic magnetic component of thermal
transport znp  interacting with the phononic one.

Besides steady state studies of the thermal conductivity, the flash
methodparker  provides information on the dynamic (in time) propagation of
heat and in particular on the interaction between the magnetic and phononic
components flash . These studies could provide a playground for confronting
experiments to recent theoretical developments on the far-out of equilibrium
dynamics in integrable spin Hamiltonians ghd .

However, the magnetic excitations in a quantum magnet are interacting with the
phonons and disentangling their contribution in the total thermal conductivity
is an important issue. A minimal phenomenological framework of thermal
conduction in a diffusive magnon plus phonon system was proposed by Sanders
and Walton (SW) sw  and it has become the standard model for analyzing
magnetic thermal transport experiments ott . In the SW model the magnetic
subsystem is assumed diffusive so it is interesting to revisit this model in
the case of a ballistic magnetic component more appropriate for the spin-1/2
Heisenberg chain compounds and not only. In the following we will try to
highlight the differences in the expected effective thermal conductivity and
thermal pulse propagation between the two models, namely the two-temperature
SW diffusion model and the present advection-diffusion one.

In SW the starting relation is the equilibration in time of the phonon
($T_{p}$) and magnetic ($T_{m}$) temperature difference,

$\frac{\partial\Delta T}{\partial t}=-\frac{\Delta T}{\tau},~{}~{}~{}\Delta
T=T_{p}-T_{m}$ (1)

where $\tau$ is a characteristic phenomenological relaxation time. This basic
relation is also satisfied by a system of individual contributions,

$\displaystyle\frac{\partial T_{p}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{c_{m}}{C}\frac{T_{m}-T_{p}}{\tau}$
$\displaystyle\frac{\partial T_{+}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{c_{p}}{C}\frac{T_{p}-T_{+}}{\tau}$
$\displaystyle\frac{\partial T_{-}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{c_{p}}{C}\frac{T_{p}-T_{-}}{\tau},$ (2)

from right (left) moving magnetic carriers with temperature $T_{\pm}$
($T_{m}=(T_{+}+T_{-})/2$) and phonons at temperarure $T_{p}$. Here, $c_{\pm}$
are the corresponding magnetic ($c_{m}=c_{+}+c_{-}$) and $c_{p}$ phonon
specific heats, $C=c_{p}+c_{m}$ the total specific heat.

These relations are for space independent temperature profiles. We can extend
them to a space dependent energy diffusion equation for the phonon subsystem
and two advection equations for the ballistic magnetic system,

$\displaystyle\frac{\partial\epsilon_{p}}{\partial
t}=D\frac{\partial^{2}\epsilon_{p}}{\partial
x^{2}}+\frac{c_{p}c_{m}}{C}\frac{T_{m}-T_{p}}{\tau}$
$\displaystyle\frac{\partial\epsilon_{\pm}}{\partial t}\pm
v\frac{\partial\epsilon_{\pm}}{\partial
x}=\frac{c_{p}c_{\pm}}{C}\frac{T_{p}-T_{\pm}}{\tau}.$ (3)

$D$ is the phonon diffusion constant, $v$ the characteristic velocity of
magnetic excitations, $\epsilon_{p}$ the phonon energy density,
$\epsilon_{\pm}$ the magnetic ones and $\delta\epsilon_{p,\pm}=c_{p,\pm}\delta
T_{p,\pm}$. To have a concrete model in mind, for the low energy spinon gas in
the 1D spin-1/2 Heisenberg chain with energy dispersion $\epsilon_{\rm
spinon}\simeq v|p|$, $\epsilon_{\pm}=\frac{\pi}{12}\frac{T_{\pm}^{2}}{v}$,
$c_{\pm}=\frac{\pi}{6}\frac{T_{\pm}}{v}$ and $v$ the spinon velocity. Here and
in the following, we consider small deviations from thermal equilibrium and
take the specific heats independent of temperature.

Furthermore, the total energy current density $Q=Q_{p}+Q_{m}$ is given by the
phonon $Q_{p}$ and magnetic $Q_{m}$ energy currents,

$\displaystyle Q_{p}$ $\displaystyle=$ $\displaystyle-\kappa_{p}\frac{\partial
T_{p}}{\partial x}$ $\displaystyle Q_{m}$ $\displaystyle=$ $\displaystyle
v\epsilon_{+}-v\epsilon_{-},$ (4)

with $\kappa_{p}=c_{p}D$ the phonon thermal conductivity. Reverting back to
temperature dependent equations,

$\displaystyle\frac{\partial T_{p}}{\partial
t}=D\frac{\partial^{2}T_{p}}{\partial
x^{2}}+\frac{c_{m}}{C}\frac{T_{m}-T_{p}}{\tau}$ $\displaystyle\frac{\partial
T_{\pm}}{\partial t}\pm v\frac{\partial T_{\pm}}{\partial
x}=\frac{c_{p}}{C}\frac{T_{p}-T_{\pm}}{\tau}$ $\displaystyle
Q=-\kappa_{p}\frac{\partial T_{p}}{\partial x}+v\frac{c_{m}}{2}(T_{+}-T_{-}).$
(5)

We will first consider the effective thermal conductivity of a system
$-L/2<x<L/2$ in a steady state with only a phonon energy currrent $Q|_{x=\pm
L/2}=Q_{p}$ at its borders and no magnetic current $(T_{+}-T_{-})|_{x=\pm
L/2}=0$. The steady state equations become,

$\displaystyle v\frac{\partial T_{m}}{\partial x}=-\frac{c_{p}}{C}\frac{\Delta
T_{m}}{\tau},~{}~{}~{}\Delta T_{m}=\frac{T_{+}-T_{-}}{2}$ $\displaystyle
D\frac{\partial^{2}T_{p}}{\partial
x^{2}}+\frac{c_{m}}{C}\frac{T_{m}-T_{p}}{\tau}=0,$ (6)
$Q=-\kappa_{p}\frac{\partial T_{p}}{\partial x}+vc_{m}\Delta T_{m}.$ (7)

Solving (7) for $\Delta T_{m}$ we obtain,

$\displaystyle\Delta T_{m}$ $\displaystyle=$
$\displaystyle\frac{1}{vc_{m}}(Q+\kappa_{p}\frac{\partial T_{p}}{\partial x})$
$\displaystyle\frac{\partial T_{m}}{\partial x}$ $\displaystyle=$
$\displaystyle-\frac{1}{\tilde{\kappa}_{m}}(Q+\kappa_{p}\frac{\partial
T_{p}}{\partial x}),$ (8)

where ${\tilde{\kappa}}_{m}=(c_{m}/c_{p})Cv^{2}\tau$ is an effective transfer
thermal conductivity. We solve this equation by, (i) assuming that
$c_{p},c_{m}$ are temperature independent (which strictly speaking is not the
case) and (ii) taking the boundary condition $T_{p}(x=0)=T_{m}(x=0)=T_{0}$. We
find,

$T_{m}=T_{0}-\frac{\kappa_{p}}{\tilde{\kappa}_{m}}(T_{p}-T_{0})-\frac{1}{\tilde{\kappa}_{m}}Qx$
(9)

and by substituting in (6),

$\displaystyle\frac{\partial^{2}(T_{p}-T_{0})}{\partial
x^{2}}-A^{2}(T_{p}-T_{0})-\frac{A^{2}}{\kappa_{t}}Qx=0$ $\displaystyle
A^{2}=\frac{c_{p}c_{m}}{C\tau}\cdot\frac{\kappa_{t}}{\kappa_{p}{\tilde{\kappa}}_{m}},~{}~{}~{}\kappa_{t}=\kappa_{p}+{\tilde{\kappa}}_{m}.$
(10)

The solution of (10) with the boundary condition
$\partial(T_{p}-T_{0})/\partial x=-Q/\kappa_{p}$ gives the phonon temparature
profile,

$T_{p}=T_{0}-\frac{x}{\kappa_{t}}Q-\frac{\tilde{\kappa}_{m}}{\kappa_{t}\kappa_{p}}\frac{\sinh
Ax}{A\cosh AL/2}Q,$ (11)

and (9) the magnetic temperature one.

The effective thermal conductivity obtained from $\kappa_{eff}=-QL/\Delta
T_{p}$ sw  is given by,

$\displaystyle\kappa_{eff}$ $\displaystyle=$
$\displaystyle\kappa_{t}\Big{(}1+\frac{\tilde{\kappa}_{m}}{\kappa_{p}}\frac{\tanh(AL/2)}{AL/2}\Big{)}^{-1},$
$\displaystyle\kappa_{eff}$ $\displaystyle\sim$
$\displaystyle\kappa_{p},~{}~{}~{}AL\rightarrow 0$ $\displaystyle\kappa_{eff}$
$\displaystyle\sim$
$\displaystyle\kappa_{p}+\tilde{\kappa}_{m}=\kappa_{t},~{}~{}~{}AL\rightarrow\infty.$
(12)

The above relations are identical to those of the SW two-temperature model
with the replacement of the magnetic conductivity by the effective magnetic
transfer one $\tilde{\kappa}_{m}=(c_{m}/c_{p})Cv^{2}\tau$.

Next, we will discuss the time dependent evolution of phonon and magnetic
temperature profiles (5) that can be probed by the flash method parker ; flash
. Considering an open system, $0<x<L$, we set as zero energy current boundary
conditions, $\partial T_{p}/\partial x|_{x=0,L}=0$ and $T_{+}=T_{-}|_{x=0,L}$.
We seek solutions of the form,

$\displaystyle T_{p}$ $\displaystyle=$
$\displaystyle\frac{a_{0}}{2}+\sum_{n=1}^{+\infty}a_{n}\cos
q_{n}x,~{}~{}~{}q_{n}=\frac{\pi n}{L}$ $\displaystyle T_{\pm}$
$\displaystyle=$ $\displaystyle\frac{b_{0}}{2}+\sum_{n=1}^{+\infty}b_{n}\cos
q_{n}x\pm c_{n}\sin q_{n}x.$ (13)

By substituting (13) in (5), we obtain the time dependence of
$a_{n},b_{n},c_{n}$,

$\displaystyle{\dot{a}}_{0}+\frac{c_{m}}{C\tau}(a_{0}-b_{0})=0$
$\displaystyle{\dot{b}}_{0}+\frac{c_{p}}{C\tau}(b_{0}-a_{0})=0$ (14)

with solution,

$\displaystyle
a_{0}(t)=a_{0}(0)-\frac{c_{m}}{C}(a_{0}-b_{0})|_{t=0}(1-e^{-t/\tau})$
$\displaystyle
b_{0}(t)=b_{0}(0)-\frac{c_{p}}{C}(b_{0}-a_{0})|_{t=0}(1-e^{-t/\tau}).$ (15)

For finite wavevector $q_{n}$,

$\displaystyle{\dot{a}}_{n}+(Dq_{n}^{2})a_{n}+\bar{c}_{m}(a_{n}-b_{n})=0$
$\displaystyle{\dot{b}}_{n}+(vq_{n})c_{n}+\bar{c}_{p}(b_{n}-a_{n})=0$
$\displaystyle{\dot{c}}_{n}-(vq_{n})b_{n}+\bar{c}_{p}c_{n}=0,$ (16)

where $\bar{c}_{p}=c_{p}/(C\tau),~{}~{}~{}\bar{c}_{m}=c_{m}/(C\tau)$ are
$O(1/\tau)$. Solutions of the form $e^{\lambda t}$ are obtained by solving the
characteristic 3rd order polynomial equation,

$\displaystyle{\tilde{\lambda}}^{3}+(Dq_{n}^{2}+\bar{c}_{m}-\bar{c}_{p}){\tilde{\lambda}}^{2}+((vq_{n})^{2}-\bar{c}_{p}\bar{c}_{m}){\tilde{\lambda}}+$
$\displaystyle(Dq_{n}^{2}+\bar{c}_{m}-\bar{c}_{p})(vq_{n})^{2}=0$
$\displaystyle\lambda={\tilde{\lambda}}-\bar{c}_{p}.$ (17)

The roots of this polynomial, although known cardano , are not physically
transparent. To proceed, we will make the physical assumption that the
relaxation time $\tau$ is the shortest scale in the problem. Thus we expect
two roots of the order $1/\tau$ presenting the relaxation of the magnetic
excitations and one root of the order of the diffusion constant. Once the
roots determined, the constants $\alpha_{i},\beta_{i},\gamma_{i}$ of the time
evolutions,
$a_{n}(t)=\sum_{i=1}^{3}\alpha_{i}e^{\lambda_{i}t},~{}~{}b_{n}(t)=\sum_{i=1}^{3}\beta_{i}e^{\lambda_{i}t},~{}~{}c_{n}(t)=\sum_{i=1}^{3}\gamma_{i}e^{\lambda_{i}t}$
are evaluated from (17) and the initial conditions
$a_{n}(t=0),b_{n}(t=0),c_{n}(t=0)$.

As an example, in the recent experiment flash  the relevant quantities for
SrCuO2 are, $\kappa\simeq 50$ W/Km, $\kappa_{p}\simeq 8$ W/Km, $c_{p}\simeq
2.8\cdot 10^{6}$ J/Km3, $c_{m}\simeq 3\cdot 10^{4}$ J/Km3 ($c_{p}\sim
100c_{m}$), $D=\kappa_{p}/c_{p}\sim 3\cdot 10^{-6}$ m2/s, $v\simeq 2\cdot
10^{4}$ m/s , $\tau\sim O(10^{-12}$ s), which gives ${\tilde{\kappa}}_{m}\sim
O(10$ W/Km). The largest uncertainty in these parameters is in the relaxation
time $\tau$. For a typical sample of length $L\sim 1$ mm these imply three
well separated characteristic time scales, $D/L^{2}\sim 3$ s-1, $(v/L)\sim
2\cdot 10^{7}$ s-1 and $1/\tau\sim 10^{12}$ s-1. For these parameters, $AL\sim
O(10^{5})$ so that $\kappa_{eff}\simeq\kappa_{p}+\tilde{\kappa}_{m}$.

Furthermore, for these experimental values, keeping the dominant terms in (17)
(e.g. dropping the 1st term, taking $Dq_{n}^{2}\tau\rightarrow 0$ and
$c_{m}<<c_{p}$) we obtain ,

$\displaystyle\lambda_{2}$ $\displaystyle\simeq$
$\displaystyle-\frac{1}{\tau}+\frac{(vq_{n})^{2}}{{\bar{c}}_{m}}$
$\displaystyle\lambda_{3}$ $\displaystyle\simeq$
$\displaystyle-\frac{c_{p}}{C\tau}-\frac{(vq_{n})^{2}}{{\bar{c}}_{m}}.$ (18)

Next, taking $\lambda\sim O(\epsilon)$, $\epsilon\tau<<1$, substituting in
(17) $\tilde{\lambda}=\epsilon+\bar{c}_{p}$ and keeping 1st order terms in
$\epsilon$, we find,

$\lambda_{1}\simeq-(D+\frac{\bar{c}_{m}}{\bar{c}_{p}}\cdot\frac{v^{2}}{\bar{c}_{p}})q_{n}^{2}.$
(19)

This relation implies a total diffusion constant composed of a phononic
component $D\simeq 3\cdot 10^{-6}$ m2/s enhanced by the ballistic magnetic
component
$\frac{\bar{c}_{m}}{\bar{c}_{p}}\cdot\frac{v^{2}}{\bar{c}_{p}}=\frac{c_{m}}{c_{p}}\cdot\frac{v^{2}C\tau}{c_{p}}\sim
O(10^{-6}$ m2/s). Consistently, multiplying (19) by $c_{p}$ we recover (12)
with the second term corresponding to the effective magnetic transfer
${\tilde{\kappa}}_{m}$. In this limit,
$\alpha_{1}\simeq\beta_{1}\simeq-\beta_{2}\simeq a_{n}(t=0)$ and
$\alpha_{3}\simeq\beta_{3}\simeq\gamma_{i}\simeq 0$. Thus, a flash method
experiment that probes the long time behavior of the temperature profile, when
analyzed in terms of a diffusion equation, gives an effective diffusion
constant with a phononic and magnetic contibution.

Whether we have three real roots or one real and two complex conjugate ones,
indicating oscillatory behavior, depends on the sign of the discriminant in
the roots $\lambda_{2,3}$,

$\Delta=(\bar{c}_{m})^{2}-4(vq_{n})^{2}.$ (20)

Assuming the experimental values quoted above and tuning the relaxation time,
we find oscillatory behavior (complex roots) of the magnetic component
relaxation for a window of $\tau$ less than about $\sim 10^{-9}$ sec giving,

$\lambda_{2,3}\sim-\frac{1}{\tau}\pm i(vq_{n}).$ (21)

Thus, the typical scenario emerging is that, within a time $\tau$ there is
equilibration of the magnetic temperature profile to the phononic one,
eventually with an oscillatory behavior, followed by diffusive propagation of
the combined system with an effective diffusion constant.

Note that, in a flash experiment the heat is deposited on the phonon system
and then it relaxes to the magnetic system, propagating coupled thereafter.
And, at any instant, it is the phonon temperature that it is probed. It would
be particularly interesting if for some material parameters we can visualize
the ballistic propagation of the magnetic temperature profile. It is clear
that at any time the magnetic profile separates in two wavefronts propagating
left-right while relaxing to the phonon bath. The question is whether this
two-bump profile can be manifested on the phonon temperature profile. This
would be a telltale sign of ballistic propagation. To realize such a scenario
within the diffusion time, it is favorable to have a magnetic system with
large specific heat that will rapidly propagate, thus large velocity and
relaxation time comparable to the diffusion time.

Figure 1: Temperature profiles at $t=0.01s$ for
$L=1mm,~{}~{}D_{p}/L^{2}=1s^{-1},~{}~{}c_{p}/C=0.1,~{}~{}c_{m}/C=0.9,~{}~{}v=20m/s,~{}~{}\tau=10^{-3}s$.
Initial temperature profiles,
$T_{p}=(40/\sqrt{\pi})e^{-(40(x-1/2))^{2}},~{}~{}T_{m}(t=0)=0$.

It is clear that the parameter space of magnetic materials to search for such
a behavior is very extended. Here we will discuss just an arbitrary example
with model parameters indicated in Fig.1, where we observe such an evolution
by numerical integration of (5). We first take as initial phonon temperature
profile a narrow Gaussian at the center of the sample. Note, that for these
parameters, the discriminant (20) is large and negative, indicative of a
strong oscillatory behavior.

Figure 2: Phonon temperature profiles for the same parameters as in Fig.1.
Initial temperature profiles
$T_{p}(t=0)=(40/\sqrt{\pi})e^{-(40x)^{2}},~{}~{}T_{m}(t=0)=0$. Figure 3:
Temperature profiles at $x/L=1$ for the same parameters as in Fig.1. Initial
temperature profiles
$T_{p}=(40/\sqrt{\pi})e^{-(40x))^{2}},~{}~{}T_{m}(t=0)=0$.

Next, for an initial heat pulse at the left-side of the sample deposited at
the phonon subsystem as in a flash experiment, we observe in Fig.2 the
propagation of a heat wave, very unlike to a diffusive one, due to the
hybridization of the phonon and magnetic components. It results to a
characteristic non-monotonic behavior at the right side of the sample as shown
in Fig.3. Certainly, further study is required to establish the range of
parameters where such a singular temperature evolution can be observed.

In conclusion, to appropriatelly describe the ballistic thermal transport of
magnetic excitations coupled to phonons we have introduced a phenomenological
advection-diffusion model. This model shows qualitatively different behavior
to the SW two-temperature diffusion model. In particular, it implies a greatly
enhanced effective diffusion constant for the material parameters of a recent
dynamic heat experiment due to the coupled propagation of magnetic and
phononic excitations.

Most important, this model predicts a specific two-bump counter-propagating
temperature profiles for certain material parameter range. It would by
particularly interesting to find materials with appropriate parameters in
order to observe such a behavior in a dynamic heat propagation experiment,
e.g. direct evidence of ballistic transport by spinons in spin-1/2 Heisenberg
chain compounds. Last but not least, althought the proposed model was
motivated by the ballistic magnetic transport coupled to the diffusive phonon
one in spin chains, it can be applied to any two-component system with coupled
advection-diffusion transport.

##  I Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft through Grant
HE3439/13. X.Z. acknowledges fruitful discussions with Profs. P. van
Loosdrecht and C. Hess on existing and future dynamic heat experiments.

## References

  * (1) C. Hess, Phys. Rep. 811, 1 (2019).
  * (2) N. Hlubek, P. Ribeiro, R. Saint-Martin, A. Revcolevschi, G. Roth, G. Behr, B. Büchner, C. Hess, Phys. Rev. B81, 020405 (2010).
  * (3) X. Zotos, F. Naef and P. Prelovsek, Physical Review B 55, 11029 (1997).
  * (4) W. Parker, R. J. Jenkins, C. P. Butler, and G. L. Abbott, J. Appl. Phys. 32, 1679 (1961).
  * (5) M. Montagnese, M. Otter, X. Zotos, D.A. Fishman, N. Hlubek, O. Mityashkin, C. Hess, R. Saint-Martin, S. Singh, A. Revcolevschi, P.H.M. van Loosdrecht, Phys. Rev. Lett. 110, 147206 (2013).
  * (6) O.A. Castro-Alvaredo, B. Doyon and T. Yoshimura Phys. Rev. X6, 041065 (2016); B. Bertini, M. Collura, J. De Nardis and M. Fagotti, Phys. Rev. Lett. 117, 207201 (2016).
  * (7) D.J. Sanders and D. Walton, Phys. Rev. B 15, 1489 (1977).
  * (8) A. V. Sologubenko, K. Giannò, H. R. Ott, A. Vietkine and A. Revcolevschi, Phys. Rev. B64, 054412 (2001).
  * (9) Abramowitz, Milton; Stegun, Irene A., eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover (1965), chap. 22 p. 773.