text:
\section{Introduction}

The study of random walks has a long tradition in network science \cite{Masuda2017}. Random walks are at the heart of many algorithms to uncover central nodes or communities of densely connected nodes, and they often serve as a first model to understand how the topology of a network affects diffusive processes.
Random Walks have also been studied mathematically and numerically when the underlying topology is a network enriched with additional features. An important family of models consists of temporal networks  \cite{Holme2012,Guide2016}, where edges are dynamical entities and diffusion can only propagate when they are active. In this case, random walks are affected by the interplay between the network topology and the statistical properties of the edge dynamics \cite{Starnini2012,Hoffmann2013,Delvenne2015}. Another important family is made of multiplex networks \cite{Domenico2016,Kivela2014,Aleta2018}, where different types of connections exist between the nodes, as in social  networks \cite{szell2010multirelational,Magnani2013} or transportation networks \cite{Cardillo2013,Gallotti2015} for instance. Multiplex networks have a layered organisation and are usually represented as tensors or by means of a so-called supra-adjacency matrix.  Random walks have also been studied in this context, to uncover how the presence of multiple layers affects diffusion \cite{Domenico2013,Domenico2014} or to define generalized versions of Pagerank \cite{de2015ranking}.

The main objective of this work is to explore phenomena emerging in the case of  networks that are temporal and multiplex \cite{Starnini2017}. To do so, we extend a generic model of continuous-time random walk on  networks \cite{Hoffmann2013}, where the waiting time distribution for an edge, that is the time before an edge gets active, depends solely on its layer. We study the properties of the resulting stochastic process and show that the presence of temporal heterogeneities across layers results in a competition between them and biases in the random walk dynamics. Here competition means that edges in one layer may have a higher probability to be selected by a walker than edges in another layer, due to the statistical properties of their temporal ordering. In other words, competition between layers emerges due to the temporality of the graph, and not as a model parameter as in previous works \cite{Ding2017,Gomez-Gardenes2015,Kleineberg2016,Jang2015}. In addition, we show and explain some apparently counter-intuitive situations,  such as the emergence of a cyclic, rock-paper-scissors precedence between the layers. Note that the notion of non-transitivity is well-known in statistics, and that it has mostly focused on systems having a finite set of possible states, such as in non-transitive dice \cite{Gardner1970}. Our work can be seen as an extension to continuous variables and a study of its impact on diffusion over multilayer networks.
As a second step, we  explore numerically the impact of the above mechanism on the dynamics of a walker and, specifically, we study the coverage of a walker on temporal multilayer network.

\section{Random walks on multiplex temporal networks}
We consider a diffusion process defined as follows \cite{Masuda2017}. The system is made of a given network of nodes and edges, and a random walker jumping between neighbouring nodes according to a given protocol as we now discuss. In a standard framework of temporal networks, which is a natural model for contacts in social networks, we assume that when the random walker arrives on a node, it triggers on each incident outgoing edge an associated random waiting time for so-called `activation' of the edge, and the first edge to reach activation is then followed by the random walker. The process is restarted on the new node reached by the random walker. Here activation is therefore seen as an event of infinitesimal duration, allowing the passage of the random walker. The time at which each edge reaches activation is independently drawn from its own waiting time distribution.
The walk is called
 active since waiting times are drawn once the walker arrives on a node, as opposite to the passive walk where contacts are taking place on edges regardless of the presence of a random walker. Note that, when the network is directed and has no short cycles,  a passive walk may be approximated by an active one provided that waiting time distributions are adapted accordingly \cite{Gueuning2017,Masuda2017}.
Because only the first edge to reach activation is taken by the walker, there is a underlying competition between different edges.
 Even though a general master equation can be derived  \cite{Hoffmann2013}, specific results have often been studied with homogeneous waiting time distribution \cite{Karsai2017,lala,Gueuning2015}.
 Here, we consider a
  situation where edges in the same layer of the multiplex network have the same waiting time distribution but distributions may differ for edges in different layers. The system thus exhibits intra-layer homogeneity and inter-layer heterogeneity.
\subsection{Emergence of biased paths}
Let us consider the trajectory of a random walker. Arriving on node $u$, the walker will leave through the first edge that reaches activation among the $k$ edges connected to $u$.
Denote $T_{1},\ldots,T_{k}$ the random variables associated to the waiting times associated to each edge, with distributions $f_1(t),\ldots,f_k(t)$. The transition time $T$ of the walker, defined as the time before its next jump is given by:
\begin{eqnarray}
    T&=&\min{\left(T_{1},T_{2},....,T_{k}\right)} \label{eq:T}.
\end{eqnarray}

In our case, each edge belongs to a layer and it is thus natural to determine which layer will be selected by the walker. In particular, we will be interested in a notion of precedence between the random variables of the layers.
The random variable $A$ is said to precede the random variable $B$, written $A \prec B$, if $P(A<B)>0.5$, that is if the edge associated to $A$ is more likely to reach activation before the edge associated to $B$.\\

The existence of precedence relations translates into biased paths for the random walker traveling on the network since at each step of the walk some edges may be statistically more likely to be selected by the walker. Therefore, understanding the precedence relation between the random variables associated to the dynamics of the network is of paramount interest in order to understand the resulting diffusion. In the following, we will exemplify the somewhat counter-intuitive properties of precedence,
 then we will numerically investigate its impact on the coverage of a random walker on a multilayer network with inter-layer heterogeneity.

\subsection{Basic properties of precedence }
\subsubsection{Rock-Paper-Scissors  }

One may see precedence as a relation of dominance between edges in competition to attract the random walker. As such one may expect transitivity. However it is not the case, and we encounter circular, rock-paper-scissors situations as follows.
We focus on the triangle network of Figure \ref{fig:sim} and the three following distributions with expectation equal to $1$:\\
\begin{eqnarray}\label{distr}
          X &=&\left\{
                \begin{array}{ll}
              U\left[0.65,0.75\right] &\quad \textrm{with probability} \quad 1-\frac{1}{\varphi}\\
                             U\left[0.3\varphi+0.65,0.3\varphi+0.75\right] &\quad \textrm{with probability} \quad \frac{1}{\varphi};

                \end{array}
              \right.\notag\\
              &&\notag\\
                         Y&\sim&\quad ~U\left[0.9,1.1\right] ;\\
                         &&\notag\\
         Z&\sim&\left\{
                     \begin{array}{ll}
                          U~[0.75,0.85] &\quad\quad\quad \textrm{with probability} \quad \frac{1}{\varphi} \\
                          U\left[\frac{\varphi}{5}+0.95,\frac{\varphi}{5}+1.05\right] &\quad\quad\quad \textrm{with probability} \quad  1-\frac{1}{\varphi} ,
                         \end{array}
                       \right. \notag
\end{eqnarray}
where $U[a,b]$ stands for the uniform distribution on the interval $[a,b]$, and $\varphi=\frac{1+\sqrt{5}}{2}\approx 1.618$.

It is straightforward to show that $Y\prec X$, $Z\prec Y$ and $X\prec Z$. As a consequence, a walker jumping on the network will have a tendency to jump clockwise on the triangle, as illustrated numerically in Figure \ref{fig:sim}.
The emergence of a circular flow leads to correlations between the successive edges on the random walk trajectory, even if edges are chosen independently at random at each step.
Competition becomes  more complex when more than two edges interact together.
For instance, in the above example, even if the pairwise precedence is uniform, that is
 \begin{eqnarray}
 P(Y<X)=P(Z<Y)=P(X<Z)=\frac{1}{\varphi},
 \end{eqnarray}
 when considering the competition between
3 different edges, one finds as $P(X<\min(Y,Z))>\frac{1}{3}$, and thus $X$ will tend to be favoured by the walker against the others two in the presence of the three types of edges simultaneously

\begin{figure}[!h]
\centering\includegraphics[scale=1]{trian}
	\centering\includegraphics[scale=0.5]{1d3}
\caption{Numerical simulations of a random walk on a multilayer triangle. Each (colored) pairwise relation belongs to a different layer and distributions activity of the layers are given by the distribution of (\ref{distr}). The walker position is incremented $+1$ at each clockwise jump, and $-1$ at each anti-clockwise jump. On the right hand side, the four dash-dot lines represent four realisations of random walk starting at $0$ while the blue line represent the average over 1000 independent walks.  The precedence relation is not transitive as a random walker jumping on the triangle will have a tendency to jump clockwise.}
\label{fig:sim}
\end{figure}
\newpage
\subsubsection{There is no most preceding edge }

It is clear that comparing either means or variances of two random variables is not sufficient to determine which one precedes the other, as the random variables $X$, $Y$ and $Z$ in the previous example share equal mean but have different variance.
Even more, it is impossible to find a random variable that precedes any random variable with equal mean.
Without loss of generality, we  prove this statement for random variables with a mean of $1$.
First, we observe that for any random variable $X$ following distribution $f(t)$, it is possible to find a random variable $Y_n$ of distribution $g_n(t)$  such that $Y_n$ precedes $X$, by setting \begin{eqnarray}
g_n(t)= &  &\left\{
 \begin{array}{ll}
  \displaystyle\frac{n}{2} (n-1) & \mbox{if } t\in\left[\frac{1}{n}-\frac{1}{n^2}~,~\frac{1}{n}+\frac{1}{n^2}\right] \\
  \displaystyle\frac{1}{2 \sqrt{n}}&\mbox{if } t\in\left[n-1+\frac{1}{n}-\frac{1}{\sqrt{n}}~,~ n-1+\frac{1}{n}+\frac{1}{\sqrt{n}}\right], \\
 0 & \mbox{otherwise}
  \end{array}
 \right.
\end{eqnarray}
 where $n\in\mbox{I\hspace{-.15em}N}_0$ is chosen large enough. One possible choice of $n$ is such that $\frac{1}{n}+\frac{1}{n^2}<\beta$, where $\beta$ is the $10\%$ quantile of $f$.
From the sequence of random variables $\left(Y_n\right)_{n\in \mbox{I\hspace{-.15em}N}}$ with distribution respectively given by $\left(g_n\right)_{n\in \mbox{I\hspace{-.15em}N}}$, it is possible to extract a subsequence of random variables $\left(Z_n\right)_{n\in \mbox{I\hspace{-.15em}N}}$ whose respective distributions have non-overlapping supports. Such a sequence is increasingly precedent, that is, for any $n \in \mbox{I\hspace{-.15em}N}$ we have $Z_{n+1} \prec Z_n$.

As the series $\left(g_n\right)_{n\in \mbox{I\hspace{-.15em}N}}$ does not converge to a probability distribution, there exists no most preceding random variable.
In terms of diffusion on multiplex networks, this results implies that there is no optimal distribution ensuring that a given layer always captures a majority of the random walk flow, independently of the distributions in the other layers.
After the choice of the other layers have been made, a layer may always find a  distribution that will allow it to be the most precedent.

\subsubsection{Layer Precedence and Node Out-degree}

The notion of precedence between random variables naturally extends to precedence between layers when waiting time distributions inside a layer are homogeneous but vary across layers. In this case, the layer $L_1$ precedes the layer $L_2$ at node $u$ if the walker sitting on $u$ is more likely to perform her next jump through an edge of $L_1$, that is if $L_1$ is more likely to capture the flow passing through $u$. Moreover, in the presence of several layers, $L_1$ precedes $L_2$ and $L_3$ jointly if $L_1$ precedes the artificial layer $L_2 L_3$ obtained by the union of the layers $L_2$ and $L_3$.

Precedence of layers at a node depends on its out-degree $k_i$ on each layer $L_i$, as the layer of the edge selected by a walker is determined by comparing the smallest time $m_i$ to reach activation on each layer $L_i$, where $m_i$ is the minimum of $k_i$ random variables with identical distributions associated to the layer $L_i$. In the simplest case of two duplicate layers $L_1$ and $L_2$ with waiting time on edges of the type $X$ and $Y$ respectively, that is $L_1$ and $L_2$ share the same nodes and have the same edges, the out-going edges consist of $k$ edges of each type, and $L_1$ precedes $L_2$ if   \begin{eqnarray}
\displaystyle\min_{j=1,\ldots,k}X_j \prec \displaystyle\min_{j=1  \ldots,k}Y_j ,
\end{eqnarray}
where $X_j$ and $Y_j$ are duplicates of random variables following the same distribution than $X$ and $Y$ respectively.

The layer precedence may therefore vary between nodes depending on their out-degree, in particular for high and low out-degree nodes. However, it is worth to note that there always exists a threshold out-degree $\nu^*$ above which one layer will always precede the other one. Indeed, the distribution $f$ is said to have a larger minimal weight than the distribution $g$ if there exists $\epsilon > 0$ such that $P(X<\epsilon) > P(Y<\epsilon)$ and $\forall ~0< \sigma < \epsilon$, $P(X<\sigma) \geq P(Y<\sigma)$, where $X$ and $Y$ are two random variables following the distribution $f$ and $g$ respectively. Then, there always exists a threshold out-degree $\nu^*$ above which the minimum of at least $\nu^*$ realizations of the distribution with the larger minimal weight will precede the other one.\\

In order to investigate this effect in a real-world setting, we construct a multiplex network as follows.
We use a dataset of private messages sent on an online social network at the University of California Irvine \cite{Pietro2009}.
This typical social network is then duplicated to create a hypothetical two-layer social network, where each layer can be thought of as corresponding to a different medium of communication.
Each layer is thus identical, and is composed of  $1899$ nodes and $20296$ edges.
 The difference between the layers is induced by the choice of two different distributions of $X$ and $Y$ defined in Eq \ref{distr}, where $Y\prec X$. As $X$ has a larger minimal weight than $Y$, there exists a switch in the precedence between the variables of type $X$ and $Y$. In this case, it is straightforward to show that this switch occurs when at least two edges of each type are competing, that is \begin{eqnarray}\displaystyle\min_{i=1,\ldots,\nu}X_i \prec \displaystyle\min_{i=1  \ldots,\nu}Y_i~~ \forall \nu\leq2.\end{eqnarray}
 Figure \ref{fig:switch_precedence} shows the probability of taking an edge of type $X$ instead of $Y$ with respect to the out-degree of a node in the static aggregated network (by construction, twice the out-degree in each layer), where each edge has two activation times associated to the random variable $X$ and $Y$ respectively. The numerical results confirm a switch of the precedence when the out-degree of the nodes increases.\\
\begin{figure}[h!]
	\centering\includegraphics[scale=0.45]{errorbar_precedence}
	\caption{ We display the average probability ($\pm$ the standard deviation) of taking an edge of $L_1$ versus one of $L_2$, with respect to the out-degree of each node on the replicate layers. We observe a switch in the layers precedence: as the nodes of out-degree equal to one tend to favor the layer $L_2$, nodes of out-degree $\geq2$ favor the layer $L_1$. The distribution of the waiting times of the edges of $L_1$ and $L_2$ are the distributions of the random variables $X$ and $Y$ respectively, defined in Eq \ref{distr}. Dashed line corresponds to a probability equal to $0.5$. }
	\label{fig:switch_precedence}
\end{figure}
\section{Impact of Inter-layer Heterogeneity on the Coverage of the Walker} \label{Mixing}

Finally, we  investigate numerically the impact of  competition in the case of diffusion in a multiplex social network with more than two layers. To do so, we use publicly available data introduced in \cite{Magnani2013} where the layers consist in five kinds of social  relationships between employees. In the numerical simulations, we  consider the three types of random variables defined in Eq.\ref{distr} and focus on the layers associated to Facebook, Leisure and Co-authorship relations. The structure of the graph thus corresponds to real-world interactions, while the inter-events time are chosen for the sake of illustration. Indeed, the presence of three layers and the choice of the specific distributions allow the emergence of properties of precedence we have previously shown, such as a rock-paper-scissors situation, and thus enables the investigation of its impact on the random walk.

\begin{figure}[h!]
	\centering
	\subfloat{
	\includegraphics[scale=0.35]{walk_real_net_crop}
	\label{fig:walk_real_crop}
	}
	\subfloat{
	\includegraphics[scale=0.34]{exploratory_deg1_crop}
	\label{fig:deg1}
	}
	\caption{Exploration rates over $100$ paths starting at each node of the multilayer network. Left: real-world network \cite{Magnani2013} with layers corresponding to Facebook, Leisure and Co-authorship relations respectively; Right: artificial network with out-degree of each node on each layer is set ton one. The waiting time distributions of the layers come from Eq.\ref{distr} and are provided in parenthesis. The exploration rate is larger under homogeneous inter-layer activity (plain lines) than under heterogeneous inter-layer activity (dashed lines) [left], except when out-degrees are low, as rock-paper-scissors situations emerge [right].\\
	}
	\label{fig:exploration}
\end{figure}

We  compute numerically the exploration rate of a walker, or coverage, defined as the percentage of nodes visited by the walker  over time \cite{Domenico2014}. We always assume homogeneous waiting time distribution inside a given layer, and consider homogeneous as well as heterogeneous inter-layer distributions.

 A typical simulation is given in Figure \ref{fig:walk_real_crop}, where we observe that the exploration rate of a random walk tends to grow faster under inter-layer homogeneity compared to inter-layer heterogeneity, irrespective of the choice of the unique distribution.
The slow-down induced by the inter-layer heterogeneity is mainly due to the fact that the flow is captured inside a single layer. Indeed, when graph density of the layers is large enough, one layer precedes the others two jointly, the one with edges of type $X$ in our example. Since the random walker has a tendency to stay in this layer, the walk will mainly take place on this layer, that is less connected than the aggregated network. However, when the graph density of the layers is weak, a rock-paper-scissors situation may arise, as the switch in precedence between layers might not occur for lower out-degree nodes. This emergence promotes the switch of the walker between different layers through lower out-degree nodes, resulting in a more efficient exploratory walk across the network, as illustrated in Figure \ref{fig:deg1} where each node has one outgoing edge in each layer, ensuring low density and homogeneous out-degree distribution. In this case, the exploration rate of the walk is similar to the exploration rate of an homogeneous inter-layer network (or monolayer network), which is the most efficient in terms of coverage rate.

 \section{Discussion}

The main purpose of this work is to investigate the competition between different layers of a multiplex network in situations when the network is temporal.
 In our framework, edges activations are modeled as independent renewal processes, each layer being characterised by a different distribution, and  we highlight implications of the concept of precedence on diffusion.  In particular, we show that precedence may lead to biases between the different layers of the network.
 Despite the simplicity of the process, it may lead to counter-intuitive properties, such as non-transitivity,
  out-degree dependence, or rock-paper-scissors situation. Our numerical results also show that precedence may have important quantitative effects on the speed of diffusion on a multilayer network, as the precedence of one layer over others may hinder the number of edges available to the walker, and hence slow down its exploration of the graph.
We also show that high out-degree nodes are more prone to favor one single layer, while low out-degree nodes exhibit a different effect, and may lead to  a cyclic exploration between the layers.
We study the impact of the precedence on an active random walk, however it is worth noticing that its impact on the passive random walk may lead to opposite bias towards layers. Indeed, in the presence of short cycles for the passive walker, an additional competition will emerge induced by the short cycles \cite{Petit2018}. For non-exponential distributions, this effect results in a backtracking bias towards or against the last traveled edges \cite{Gueuning2017}, typically leading to the emergence of short cycle patterns in human-related network \cite{Saramaki2015}.

  These results remain mostly mathematical as we used toy-model distributions instead of ones modelled on real-life data, and an important next step would be to test the resulting ideas in empirical data of multiplex networks where different layers are associated to different time scales, for instance between physical, mobile phone and social media interactions \cite{sekara2016fundamental}.

\section*{Acknowledgment}
This work has been partially supported by the Concerted  Research Action (ARC) supported by the Federation Wallonia-Brussels Contract ARC 14/19-060  ; and  Flagship European Research Area  Network  (FLAG-ERA)  Joint  Transnational  Call  "FuturICT  2.0".

\bibliographystyle{imaiai}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction\label{Introduction}}

It is well known that the toplogy of the network changes the dynamics
dramatically. After the discovery of the absence of epidemic threshold
in a scale-free network\cite{Pastor-Satorras2001}, many researchers have
 focused on the dynamics on complex network, such as synchronization
\cite{Nishikawa2003, Ichinomiya2004}, pattern formation\cite{Nakao2008},
 and other phenomena.

In this paper, we study  the Bouchaud-M\'ezard (BM)
 model on complex network\cite{Bouchaud2000}.
The power-law behavior of wealth distribution, Pareto's
distribution, has for over a century been one of the main
cornerstones of econophysics\cite{Pareto}. One of the simplest
models that explain this law is that proposed by Bouchaud and
M\'ezard, which consists of multiplicative noise
and globally coupled diffusion. After the proposal of the
BM model, several researchers numerically
investigated the generalized BM model in which diffusion occurs
between adjacent nodes in a complex network\cite{Garlaschelli2004,
Souma2001}. While, the research revealed that the network topology
alters the wealth distribution, there has been no quantitative
theory that sufficiently explains these simulation results.

Recently, we proposed a new theory for the BM model on a random
network\cite{Ichinomiya2012}. Using several assumptions that we
describe in the next section, we derived equations that determine
the static probability density function (PDF) of the wealth. The
results of this analysis were compared with those of the numerical
simulations and good agreement was obtained.

However, the wealth distributions were analyzed only on a random
network in the previous paper, and the distributions on other
complex networks was left as an open problem. The aim of this paper
is develop our previous work so that it is applicable to a general
complex network. Using the same techniques applied in the previous
paper, we derive the equations that determine the static PDF for the
BM model on a complex network for a given adjacency matrix. The
results are evaluated by a comparison with the numerical simulation
results, and our method is revealed to perform well for many network
systems.

This paper is organized as follows. In section \ref{Theory}, we
derive the self-consistent equations that determine the PDF of the
BM model on a complex network whose adjacency matrix is given. The
results obtained with this PDF are tested in sec. \ref{Simulation}
by comparing with simulations on several networks such as a random
network, the Watts--Strogatz (WS) network, and a real social
network. Finally, we discuss the results and future problems and
summarize the paper in sec. IV.

\section{Theory\label{Theory}}

We consider the BM model on a undirected complex network consisting
of $N$ nodes, whose adjacency matrix is given by $A=(a_{ij})$. The
dynamics of $x_i$, i.e., the wealth on node $i$, are determined by
the following Ito-type stochastic differential equations.
\begin{equation}
 dx_i = J\sum_{j=1}^N a_{ij} (x_j-x_i) dt + \sqrt{2}\sigma dW_i,\label{083834_8Aug12}
\end{equation}
where $J$, $\sigma$ and $W_i$ represents the diffusion constant of
wealth, strength of noise, and the standard Brownian motion, respectively.
This equation has no static PDF; however, the normalized  wealth
$x/\langle x \rangle$, where $\langle \cdots \rangle$ represents an
average over all nodes, can have a static PDF. For simplicity, $x$
is used in place of $x/\langle x  \rangle$ in the following
discussion.

To obtain a static $x_i$ distribution, we make the ``adiabatic and
independent'' approximation introduced in our previous
paper\cite{Ichinomiya2012}. We first assume that all $x_i$ values
are independent and that the correlation between different nodes is
negligible. Under this assumption, we assume the PDF
can be decomposed as $\rho(x_1,x_2,\cdots,
 x_N)=\rho_1(x_1)\rho_2(x_2)\cdots \rho_N(x_N)$.  Second, we assume that the rate of change of
$\bar x_i =\frac{1}{d_i}\sum_j a_{ij}x_j$, the average $x$ around
node $i$, is much slower than that of $x_i$, where $d_i$ represents
the degree of node $i$. Using these assumptions, we can calculate
the static PDF as follows.

Suppose $\bar{x}_i$ is constant. Then  $\rho_i(x|\bar{x}_i)$,  the
conditional PDF of $x$ on node $i$ for a given $\bar{x}_i$, is a
solution of the following Fokker--Planck equation.

\begin{equation}
 \frac{\partial \rho_i}{\partial t} = -\frac{\partial}{\partial
  x}[d_i(\bar{x}_i-x)\rho_i] +\sigma^2 \frac{\partial^2}{\partial x^2}
  [x^2 \rho_i]\;.
\end{equation}

The static solution of this equation, $\rho_{eq,i}(x|\bar{x}_i)$, is
given by
\begin{equation}
 \rho_{eq,i}(x|\bar x_i)=C_{d_i}(\bar x_i)\exp(-\alpha_{d_i} \bar x_i/x)x^{-2-\alpha_{d_i}}\;,\label{164143_7Aug12}
\end{equation}
where $\alpha_d =Jd/\sigma^2$ and  $C_d(x) = (\alpha_d
 x)^{1+\alpha_d}/\Gamma(1+\alpha_d)$. Using the adiabatic assumption, the static PDF of $x$ at node $i$ is given by
\begin{equation}
 \rho_{eq,i}(x) = \int d\bar{x} P_i(\bar{x})\rho_{eq,i}(x | \bar x)\;,\label{164204_7Aug12}
\end{equation}
where $P_i(\bar x)$ represents the PDF of $\bar x_i$.

To obtain $P_i(\bar x)$, we use the independent assumption again
together with the central-limit theorem. Under this assumption,
$P_i(\bar x)$ can be assumed to be the distribution of the average
of the independent variables $x_j$, where $j$ represents a node
adjacent to node $i$. If the variances of $x_j$ are finite for all
$j$, we can apply the central-limit theorem to approximate $P_i(\bar
x)$ by the Gaussian distribution
\begin{equation}
 P_i(\bar x ) = \frac{1}{\sqrt{2\pi} S_i} \exp\left[-\frac{(\bar x -M_i)^2}{2S_i^2}\right]\;,\label{091257_8Aug12}
\end{equation}
where $M_i$ and $S_i$ represent the average and variance of $\bar
x_i$, respectively. From Lindeberg's theorem, these values are given
by
\begin{equation}
 M_i = \sum_j a_{ij}\mu_j/d_i\;,\label{090008_8Aug12}
\end{equation}
and
\begin{equation}
S_i^2 = \sum_j a_{ij} s_j^2/d_i^2\;,\label{090021_8Aug12}
\end{equation}
where $\mu_j$ and $s_j^2$ are the average and variance of $x_j$,
respectively.

We can determine $M_i$ and $S_i$ by calculating $\mu_i$ and $s_i^2$
for $x$ using Eqs. (\ref{164143_7Aug12}), (\ref{164204_7Aug12}),
 (\ref{090008_8Aug12}), and (\ref{090021_8Aug12}). By using
 $\int_0^{\infty} dx \rho_{eq,i}(x|\tilde x) x = \tilde x$ and
$\int_0^{\infty} dx \rho_{eq,i}(x|\tilde x) x^2  =
\frac{\alpha_{d_i}}{\alpha_{d_i}-1}\tilde x^2$ for $\alpha > 1 $, we
find
\begin{equation}
 \int_0^\infty dx x\rho_i(x) = \int_{0}^{\infty} d\bar{x} P_i(\bar{x})
  \int_0^\infty dx x \rho_i(x|\bar{x}) =\int_{0}^{\infty} d\bar{x} \bar
  x P_i(\bar{x}) \sim M_i\;,\label{170529_22Aug12}
\end{equation}
and
\begin{equation}
  \int dx x^2 \rho_{i}(x) = \int_0^{\infty} d\bar x
  \frac{\alpha_{d_i}}{\alpha_{d_i}-1} \tilde x^2 P_i(\bar x) \sim
  \frac{\alpha_{d_i}}{\alpha_{d_i}-1} \left(M_i+S_i^2\right)\;,\label{170556_22Aug12}
\end{equation}
where we assume $M_i \gg \sqrt{S_i}$. We thus obtain $M_i=\mu_i=1$
for all $i$, and $S_i^2$ is the solution of
\begin{equation}
 s_i^2=\frac{\alpha_{d_i} S_i^2 + 1}{\alpha_{d_i} - 1}\;.\label{133408_24Aug12}
\end{equation}
Using Eq. (\ref{090021_8Aug12}), we can rewrite this equation as
\begin{equation}
 s_i^2=\frac{1}{\alpha_{d_i}-1} \left(1+\frac{\alpha_{d_i}}{d_i^2}\sum_j a_{ij} s_j^2\right)\;.\label{175036_7Aug12}
\end{equation}
For convenience and in anticipation for later analysis, we define a
new matrix $R=(r_{ij})$ as
\begin{equation}
 r_{ij}= \frac{\alpha_{d_i}}{d_i^2(\alpha_{d_i}-1)} a_{ij}\;.\label{defOfR}
\end{equation}
We then finally obtain the following equation for $s_i^2$:
\begin{equation}
 \sum_j (\delta_{ij}-r_{ij})s^2_j = 1/(\alpha_{d_i}-1)\label{eqForS}\;.
\end{equation}
Eqs. (\ref{164143_7Aug12}), (\ref{164204_7Aug12}),
(\ref{091257_8Aug12}), and (\ref{eqForS}) define the PDF of $x$.

Before concluding this section, we note in the following about wealth
condensation, defined as the divergence of $s_i^2$. Because $s_i^2
\ge 0$ for all $i$, Eq. (\ref{eqForS}) must have a solution
$s_j^2>0$ to have any meaningful result. If $J$ is large enough,
this condition is satisfied provided that $\mbox{min}(d_i)\ge 2$.
For a $J$ value that satisfies $\alpha_{d_i}>1$ for all $i$,
$R_{ij}$ is a non-negative matrix. Therefore, the solution of Eq.
(\ref{eqForS}) satisfies $s_i^2
>0$ if $\rho(R) <1$, where $\rho(R)$ represents the spectral radius
of $R$. From Frobenius's theorem, $\rho(R)$ is equal to the largest
eigenvalue of $R$, and thus we conclude that Eq. (\ref{eqForS}) has
a meaningful solution if the largest eigenvalue of $R$ is less than
1. Noting that $d_i = \sum_{j} a_{ij}$,  $\rho(R) \le
\mbox{max}(\sum_j |r_{ij}|)$, and $r_{ij} \rightarrow a_{ij}/d_i^2$
for $J\rightarrow \infty$, $\rho(R)$ is less than 1 for large $J$,
and a non-condensed phase appears, provided that $\mbox{min}(d_i)
\ge 2$.

In contrast, we always have a condensed phase at small $J$ due to
the divergence of the $r_{ij}$ elements. As we decrease $J$ from a
large value, the $r_{ij}$  values increase and diverge when
$\alpha_{d_i} = 1$; $\rho(R)$ also diverges in this case. Therefore,
the wealth is condensed if $J/\sigma^2  \le \mbox{max}(1/d_i)$.

\section{Simulations on various network models\label{Simulation}}

To test our theory, we compare the results with numerical simulation
results on several networks. We begin with the network created by
the Erd\"os-R\'enyi algorithm and then consider several networks
created by the WS model\cite{Watts1998} to test the effect of
clustering. It is shown that our theory fails when the rewiring rate
$p$ is too small, but good agreement is obtained when $p\ge 0.1$.
Finally, we test our theory on a real network. It would be ideal to
apply the theory to a real economic network, but as there is
currently no available data, we use the American college football
network obtained by Girvan and Newman\cite{Girvan} instead. For this
case, the PDF obtained from the simulations shows good agreement
with our theory. For all the simulations, $\sigma^2$ is set to 1.

\subsection{Random network}

The network created by the Erd\"os-R\'enyi algorithm, shown in Fig.
\ref{140327_24Aug12}, has 200 nodes and a mean degree of $\langle k
\rangle = 10$.
\begin{figure}[t]
 \resizebox{!}{!}{\includegraphics{randGraph.eps}}
 \caption{The random network used for the simulations in sec. III.A}
 \label{140327_24Aug12}
\end{figure}
\begin{figure}[t]
 \resizebox{.3\textwidth}{!}{\includegraphics{randj03.eps}}
 \resizebox{.3\textwidth}{!}{\includegraphics{randj05.eps}}
 \resizebox{.3\textwidth}{!}{\includegraphics{randj2.eps}}
 \caption{A log-log plot of the PDF on the random network obtained by numerical simulation
 for $J=0.3$ (left), 0.5 (middle), and 2.0(right). }
 \label{170117_22Aug12}
\end{figure}

In Fig. \ref{170117_22Aug12} we plot the PDFs obtained by numerical
simulation for $J=0.3, 0.5$, and 2.0. The PDFs obtained using our
theory are indicated by the solid line and show good agreement with
the numerical simulation for all $x$ for the $J=2.0$ case, but there
are minor discrepancies at small $x$ values in the cases of $J=0.3$
and 0.5.
The discrepancy is due to the $Mi \gg
 \sqrt{S_i}$ assumption used for deriving Eqs. (\ref{170529_22Aug12}) and (\ref{170556_22Aug12}).
The maximum $S_i^2$ value calculated by Eqs. (\ref{eqForS}) and
(\ref{133408_24Aug12}) is 0.33 for $J=0.3$ and 0.09 for $J=0.5$,
while $M_i = 1$. Thus, the $M_i \gg S_i$ approximation is not valid
as $\sqrt{S_i}$ is larger than 0.3. However, the approximation is
valid for the $J=2.0$ case as the maximum $S_i^2$ value is 0.01.

\subsection{Small-World Network}

Next, we consider three networks created by the WS algorithm. The
networks have $N=100$ nodes and a mean degree of $\langle d \rangle
= 8$ with different rewiring rates of $p=0.01, 0.1$, and 0.2. Fig.
\ref{155309_15Aug12} shows both the networks and PDFs. In the case
of $p=0.01$, our theory differs from the numerical simulation, while
we find good agreement for the $p=0.1$ and $p=0.2$ cases.
\begin{figure}[t]
 \resizebox{.3\textwidth}{!}{\includegraphics{ws001.eps}}
 \resizebox{.3\textwidth}{!}{\includegraphics{ws01.eps}}
 \resizebox{.3\textwidth}{!}{\includegraphics{ws02.eps}}\\
\resizebox{.3\textwidth}{!}{\includegraphics{ws001j2.eps}}
\resizebox{.3\textwidth}{!}{\includegraphics{ws01j2.eps}}
\resizebox{.3\textwidth}{!}{\includegraphics{ws02j2.eps}}
\caption{Upper: The small-world networks created by the
Watts--Strogatz algorithm with rewiring probabilities of
$p=0.01,0.1,$ and 0.2 used to test our theory. Lower: The PDFs
obtained from simulations on each networks with $J=2.0$. The solid
lines show the results of the method proposed here.}
\label{155309_15Aug12}
\end{figure}

To investigate the reason for this discrepancy, we examine the
scatter plots of $x$ and $\bar x$ for $p=0.01$ and $p=0.1$ in Fig.
\ref{scatterPlot}. For $p=0.01$, we find that there is a strong
correlation between $x$ and $\bar x$ that weakens for $p=0.1$. The
differences can be quantified by the correlation, which is 0.74 for
$p=0.01$ and 0.48 for $p=0.1$. This result suggests that the discrepancy
is caused by the breakdown of independent assumption.

It is reasonable to question here when and why the independent
assumption fails. The small-world property does not appear to affect
the independence of wealth because the networks show high-clustering
and small betweenness for both the $p=0.01$ and 0.1 cases. We return
to this point in the discussion section.

\begin{figure}[t]
 \resizebox{.45\textwidth}{!}{\includegraphics{scatterWS001J2.eps}}
 \resizebox{.45\textwidth}{!}{\includegraphics{scatterWS01J2.eps}}
 \caption {Scatter plots of $x$ and $\bar x$ for the BM model on WS
 networks with $p=0.01$ and 0.1 and with $J=2$.}
 \label{scatterPlot}
\end{figure}

\subsection{American college football network}

The American college football network obtained by Girvan and
Newman\cite{Girvan} includes 115 nodes and 613 edges. The simulated
PDFs for $J=0.3$ and 0.5 are shown in Fig. \ref{093508_20Aug12}. In
both cases, our theory reproduces the PDF well.

The agreement between theory and simulation seems excellent. We have not
understand the reason for this good agreement, however, one possible
reason is that this network has community structure.
As Girvan and Newman showed,  this  network has clear communities.
As we will see in the later section, the spatial correlation of wealth
would be small for globally coupled network.
Therefore it is natural to assume that wealth correlation is also small, if
the network is consisted from several densely connected subgraphs.
However, we still need further study to explain this excellent agreement.

\begin{figure}[t]
 \resizebox{.45\textwidth}{!}{\includegraphics{fbj03.eps}}
 \resizebox{.45\textwidth}{!}{\includegraphics{fbj05.eps}}
 \caption{The simulated PDFs of the BM model on the American college football network for $J=0.3$ and 0.5.
 The solid lines represent the results obtained by our theory.}
 \label{093508_20Aug12}
\end{figure}

\section{Discussion and Conclusions}

Here, we have developed a theory for the BM model on complex
networks. By generalizing our previous work, we were able to propose
a theory applicable to a complex network with a given adjacency
matrix. Using the adiabatic and independent assumptions, we derived
the equations that determine the static PDF of the wealth. The
result was compared to numerical simulations, and we found our
theory works well on a Erd\"os-R\'enyi network, WS networks with $p \ge 0.1$,
and the American college football network. We also found that the
theory did not perform well for a WS network with $p= 0.01$ as the
independent assumption is no longer valid.

The theory does not apply to a WS network with a small $p$ value
because of the large spatial correlation in this case. Thus, we may
ask when and why the spatial correlation becomes large. While it is
difficult to give a complete answer to this question, the following
discussion suggests that the properties of the Laplacian are
essential to this problem.

Suppose that $J$ is very large and the fluctuations of $x$ are very
small. Under such circumstances, we can approximate $xdW$ in Eq.
(\ref{083834_8Aug12}) as $dW$ and obtain
\begin{equation}
 dx_i = - J \sum_j L_{ij} x_j  dt +\sqrt 2 \sigma dW\;,\label{150804_20Aug12}
\end{equation}
where $L=(L_{ij})$ is the Laplacian matrix $L_{ij}= \delta_{ij} d_i
- a_{ij}$. To investigate the correlation in this toy model, we
assume for simplicity that the network is undirected and connected.
Introducing the eigenvalues and corresponding normalized
eigenvectors of $L$ as $\lambda_1 = 0 <  \lambda_2 \le \cdots
\lambda_N$ and $\mathbf{v}_1,\mathbf{v}_2,\cdots \mathbf{v}_N$, we
can write $(x_1, x_2, \cdots , x_N)^T = \sum_{i=1}^N q_i
\mathbf{v}_i$. Then $L$ is diagonalized and Eq.
(\ref{150804_20Aug12}) can be written as
\begin{equation}
 dq_i = - J \lambda_i q_i dt +\sqrt 2 \sigma dW\;.
\end{equation}
The static distribution of $q_i (2\le i \le N)$ is
calculated as
\begin{equation}
 P(q_i) = \frac{\sqrt{J \lambda_i}}{\sqrt{2\pi}
  \sigma}\exp\left[-\frac{J \lambda_iq_i^2}{2\sigma^2}\right]\;.\label{133706_28Aug12}
\end{equation}
The $q_1$ term does not have a static distribution because
$\lambda_1 =0$. However, because $\mathbf{v}_1 =
\frac{1}{\sqrt{N}}(1,1,\cdots,1)$, this term gives the uniform
change, $x_i\rightarrow x_i+C$, and does not contribute to the
correlation.

From Eq. (\ref{133706_28Aug12}) and the relation $x_i = \sum_j q_j
(\mathbf{v}_j)_i$, we obtain the covariance
\begin{equation}
 \langle (x_i-\langle x \rangle)( x_j-\langle x \rangle) \rangle =
  \sum_{l=2}^N \frac{\sigma^2}{J\lambda_l}(\mathbf{v}_l)_i(\mathbf{v}_l)_j\;.\label{145452_24Aug12}
\end{equation}
Our independent assumption implies that the off-diagonal elements of
Eq. (\ref{145452_24Aug12}) are smaller than the diagonal elements.
This condition is satisfied if every $\mathbf{v}_i$  is "localized", in other
words, has only one large components.
It has been shown by some researchers that eigenvectors is localized in many complex networks \cite{McGraw2008,Nakao2010}.
Our theory will be applicable to these network models.

We also consider here the case for which the adiabatic assumption is
not valid. If the node degree is large enough, then $\tilde x$
changes at a much slower rate than that of $x$. Therefore, this
assumption will be valid if the degree of each node is large. On the
other hand, if the degree of each node is small, then not only the
adiabatic assumption but also the use of the central-limit theorem
cannot be justified. For example, our theory will not work well for
the ``star''-network in which almost all nodes have a degree of 1.

We note that theory for the BM model in the wealth-condensed phase
is still elusive, as the standard central-limit theory is not
applicable when $\langle x^2 \rangle $ diverges. Thus, to discuss
the static properties, we will need to generalize the central-limit
theorem so that it can be adopted even if the variance diverges.

Finally, we note that our developed technique is applicable to other
dynamical systems subject to noise. Stochastic dynamics on complex
network is attracting the interests of many reserchers recently, however,
dynamics with multiplicative noise has not studied so far.
We believe this method will shed new light on other problems regarding
the dynamics on complex networks.