Synchronization is optimal in non-diagonalizable networks
Takashi Nishikawa
We consider the problem of maximizing the synchronizability of oscillator networks by assigning weights and directions to the links of a given interaction topology. We first extend the well-known master stability formalism to the case of non-diagonalizable networks. We then show that, unless some oscillator is connected to all the others, networks of maximum synchronizability are necessarily non-diagonalizable and can always be obtained by imposing unidirectional information flow with normalized input strengths. The extension makes the formalism applicable to all possible network structures, while the maximization results provide insights into hierarchical structures observed in complex networks in which synchronization plays a significant role.
arxiv.tex
% January 21, 2006
\documentclass[aps,pre,twocolumn,superscriptaddress,showpacs]{revtex4}
\usepackage{amsmath,epsfig,amsfonts,amssymb,latexsym,times}
\baselineskip 25pt

\begin{document}

{\small Phys. Rev. E \textbf{73}, 065106 (2006)}
  
\title{Synchronization is optimal in non-diagonalizable networks}

\author{Takashi Nishikawa}
\affiliation{Department of Mathematics, Southern Methodist University, Dallas, TX 75275}

\author{Adilson E. Motter}
\affiliation{Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208}
\affiliation{CNLS and Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545}

\date{June 30, 2006}

\begin{abstract}
We consider the problem of maximizing the synchronizability of oscillator networks
by assigning weights and directions to the links of a given interaction topology.
We first extend the well-known master stability formalism to the case of non-diagonalizable 
networks.
We then show that, unless some oscillator is connected to all the others,
networks of maximum synchronizability are necessarily non-diagonalizable and can always be
obtained by imposing unidirectional information flow with normalized input strengths.
The extension makes the formalism applicable to all possible network structures,
while the maximization results provide insights into hierarchical structures
observed in complex networks in which synchronization plays a significant role.
\end{abstract}

\pacs{05.45.Xt, 89.75.-k, 87.18.Sn}

\maketitle

Under extensive study in recent years is how the collective dynamics of a complex
network is influenced by the structural properties of the network \cite{recent},
such as clustering coefficient~\cite{wattsstrogatz1998}, average network distance~\cite{watts1999},
%degree 
connectivity
distribution~\cite{barabasi1999},  assortativity~\cite{Newman2002}, and weight distribution
\cite{Motter2005b,Motter2005,Zhou2006}. The effects of these properties on synchronization has particularly
attracted the attention of researchers, partly because of the elegant analysis due to
Pecora and Carroll~\cite{pecora1998} which allows us to isolate the 
%effect 
contribution
of
the network structure in terms of the 
%Laplacian eigenvalues. 
eigenvalues of the coupling matrix.
%%Studying such effects not only provides us with insights into the dynamics of real-world complex
%%networks, but also guides us in designing large artificial networks.  

Synchronizability of complex networks of oscillators generally has been shown to improve as the
average network distance decreases, with one notable exception:
%In scale-free networks~\cite{barabasi1999}, 
in random scale-free networks,
which are characterized by a
strong heterogeneity of the 
%degree 
connectivity
distribution \cite{barabasi1999}, 
synchronization was shown to become
more difficult as the heterogeneity increases~\cite{Nishikawa2003}, even though the average
network distance decreases at the same time. Motivated by this counter-intuitive effect,
%of the degree heterogeneity, 
%several 
researchers have pursued
ways to enhance the synchronizability of scale-free networks by introducing directionality and
weight to each link in the network~\cite{Motter2005b,Motter2005,Chavez2005}. 
A natural question arising in this context is: \textit{Given a network of oscillators with a fixed
topology of interactions, which assignment of weights and directions maximizes
its synchronizability?}  By maximization, we mean that the synchronized states are stable for the widest
possible range of the parameter representing the overall coupling strength.


The study of such a question not only provides us with insights into the dynamics of 
real-world complex networks but also guides us in designing large artificial networks.
Metabolic networks---the system of hundreds of interconnected biochemical reactions
responsible for the biomass and energy production in a cell---is a prototypic
example where the weights and directions of feasible links (metabolic fluxes) are adjusted 
to optimize fitness, which is likely to account for robustness of synchronized behavior
against environmental changes \cite{fishcer:2005}. Other 
%relevant 
examples 
%can be found ranging 
range
from the enhancement of neuronal synchronization for a given topology of synaptic
connections in the brain, to the design of interaction schemes that optimize the performance
of computational tasks based on the synchronization of processes in computer networks 
\cite{korniss:2003}. The adjustment of flows in power grids and communication patterns
in social organizations are additional examples where directional and weighted patterns
can be favored because they can better facilitate the synchronized or coordinated behavior
on which the functioning of these networks is based.


%Such a problem arises naturally in complex networks in the real world.  
%For example, natural organisms have networks of hundreds of metabolites connected by potential biochemical reactions, the rates (weight of links) of which might have been %adjusted under evolution to optimize their fitness, which may be related to the robustness of synchronized behavior against environmental changes.  
%Another example is the network of communication traffic constrained within the network of 
%communication channels, such as neural networks in the brain or the social networks of people, 
%in which certain traffic patterns may 
%be favored because they can better facilitate synchronized or coordinated behaviors on which some functions of the network are based.

Here we show that the answer to the question of maximum synchronizability 
%typically 
falls outside the framework of the Pecora-Carroll analysis, which is built
on the assumption that the 
%Laplacian 
network dynamics can be linearly decomposed into eigenmodes, i.e., the
coupling matrix of the network is diagonalizable. Indeed, we show that maximally
synchronizable networks are always \emph{non-diagonalizable}
(except for the 
%unrealistic 
extreme configurations where a node is connected to all the others)
%unless some node 
%has the full degree.
%is directly connected to all the others.
and can be constructed for any given interaction
topology by imposing that the network: ({\it i}) embeds an {\em oriented spanning tree}, 
({\it ii}) has {\em no directed loops}, and ({\it iii}) has {\em normalized} input strength in each node.
%We also show that these optimal networks are obtained whenever there are {\em no directed loops}, 
%the total input strength of each node is {\em normalized}, and the 
%network embeds an {\em oriented spanning tree}.  
The fact that the 
%Laplacian 
%coupling matrix 
networks
are not necessarily
diagonalizable has been largely overlooked in the literature, apparently because 
most previous works have focused 
on networks of symmetrically coupled oscillators,
which are guaranteed to be diagonalizable. However, 
the same does not hold true in general when the network is directed,
as required in the realistic modeling of many complex systems.
%as observed in realistic systems.
%In fact, the non-diagonalizable networks that are optimal are quite realistic, because many networks in the real world are naturally directed.
Here we develop a new theory that extends the Pecora-Carroll
analysis to the case of non-diagonalizable networks.  
We show that in this case the
synchronizability  is still determined 
%in the same way 
by the eigenvalues of the 
%Laplacian 
coupling matrix, but the speed at which the system converges toward the synchronized
state may be significantly slower. 
This theory is a first example of going beyond the traditional framework for studying complex systems based on either decomposition into eigenmodes or some sort of superposition principle.
%{\bf The main focus of complex systems research was either on systems that can be decomposed into eigenmodes or
%systems in which some sort of superposition principle applies.}
%This theory is a first example of a general framework for treating complex systems that do not have these properties.
%{\bf The study of complex system has traditionally been focused mainly on either systems that can be decomposed into eigenmodes or
%systems in which some sort of superposition principle applies.  The theory presented below is a first general framework for treating complex systems that do not have these properties.}

Consider $n$ identical oscillators whose individual dynamics without coupling is 
governed by $\dot{\mathbf{x}} = \mathbf{F}(\mathbf{x})$, $\mathbf{x}\in I\!\! R^m$. 
Now consider the network of these oscillators coupled via an output signal function 
$\mathbf{H}: I\!\! R^m \to I\!\! R^m$ along a network with
a symmetric adjacency matrix $A=(A_{ij})$ defined by $A_{ij} = 1$ if oscillators
$i$ and $j$ ($\neq i$) are connected and $A_{ij} = 0$ otherwise.  
Let $W_{ij}\ge 0$ denote the strength of the coupling that oscillator $i$
receives from $j$.  Thus, $A$ represents the topology of interactions
and $W=(W_{ij})$ represents the assignment of weights and directions.  The system of
equations governing the dynamics of the oscillator network  
can then be written as $\dot{\mathbf{x}}_i = \mathbf{F}(\mathbf{x}_i) + \sigma \sum_{j=1}^n A_{ij}W_{ij}[\mathbf{H}(\mathbf{x}_j) - \mathbf{H}(\mathbf{x}_i)]$
or, equivalently,
%\begin{eqnarray}
%\dot{\mathbf{x}}_i &=& \mathbf{F}(\mathbf{x}_i) + \sigma \sum_{j=1}^n A_{ij}W_{ij}[\mathbf{H}(\mathbf{x}_j) - \mathbf{H}(\mathbf{x}_i)]\nonumber\\
% &=& \mathbf{F}(\mathbf{x}_i) - \sigma \sum_{j=1}^n L_{ij}\mathbf{H}(\mathbf{x}_j),\quad i=1,\ldots ,n, \label{eqn:main}
%\end{eqnarray}
\begin{equation}
\dot{\mathbf{x}}_i =\mathbf{F}(\mathbf{x}_i) - \sigma \sum_{j=1}^n L_{ij}\mathbf{H}(\mathbf{x}_j),\quad i=1,\ldots ,n,
\label{eqn:main}
\end{equation}
where $\sigma$ is the parameter controlling the overall coupling strength
and $L = (L_{ij})$ is the 
%Laplacian 
coupling
matrix of the directed weighted network, defined
by $L_{ij} = -A_{ij}W_{ij}$ if $i \neq j$ and $L_{ii} = -\sum_{j\neq i} L_{ij}$. 
Note that $L$ is not necessarily symmetric because the network is not constrained to
be undirected.

The maximization problem considered in this paper can
%then 
be formulated as follows. For a given topology of interactions 
between oscillators (represented by $A$), we want to find the 
assignment of weights and directions (represented by $W$) that
maximizes the synchronizability of the network. In order to address
this question, 
we need a condition for the network to synchronize.
For any solution $\mathbf{x}=\mathbf{s}(t)$ of the individual dynamics
$\dot{\mathbf{x}} = \mathbf{F}(\mathbf{x})$, the completely synchronous
state $\mathbf{x}_i = \mathbf{s}(t)$, $i=1,\dotsc,n$ is automatically
a solution of the entire system~\eqref{eqn:main}.
%because $\sum_j L_{ij} = 0$ for all $i$.
The question then is  to determine when
this solution is stable against small perturbations.  This synchronization 
condition can be derived by extending the linear stability analysis of 
Pecora and Carroll~\cite{pecora1998} to the case where $L$ is not
necessarily diagonalizable, as follows.

 
The starting point of our analysis is the observation that, for any $n \times n$
matrix $L$, there exists an invertible matrix $P$ of generalized eigenvectors of $L$
which transforms $L$ into Jordan canonical form as $P^{-1}LP = J$, where
\begin{equation}
\label{eqn:jordan}
J = \begin{pmatrix}
0 & & &\\
& B_1 & &\\
& & \ddots & \\
& & & B_l
\end{pmatrix},\;\;
B_j = \begin{pmatrix}
\lambda & & & \\
1 & \lambda & & \\
& \ddots & \ddots & \\
& & 1 & \lambda \\
\end{pmatrix},
\end{equation}
and $\lambda$ is one of the (possibly complex) eigenvalues of $L$. 
The stability of the synchronous solution of Eq.~\eqref{eqn:main}
is determined by the variational equation
$ \dot{\xi} = D\mathbf{F}(\mathbf{s})\xi - \sigma D\mathbf{H}(\mathbf{s})\xi L^T$,
where $\xi = (\boldsymbol{\xi}_1, \ldots ,\boldsymbol{\xi}_n)$ and $\boldsymbol{\xi}_i$
is the perturbation to the $i$th oscillator.
By applying the change of variable
$\eta = \xi P^{-T}$, we get
\begin{equation}
\label{eqn:variational2}
\dot{{\eta}} = D\mathbf{F}(\mathbf{s}){\eta} - \sigma D\mathbf{H}(\mathbf{s}){\eta} J^T.
\end{equation}
Each block of the Jordan canonical form corresponds to a subset of equations in
\eqref{eqn:variational2}. For example, if block $B_j$ is $k \times k$, then it
takes the form
\begin{align}
\dot{\boldsymbol{\eta}}_{1} &= [D\mathbf{F}(\mathbf{s})
- \alpha D\mathbf{H}(\mathbf{s})]\boldsymbol{\eta}_{1} \label{eqn:eta1}\\
\dot{\boldsymbol{\eta}}_{2} &= [D\mathbf{F}(\mathbf{s}) 
- \alpha D\mathbf{H}(\mathbf{s})]\boldsymbol{\eta}_{2} - \sigma D\mathbf{H}(\mathbf{s})\boldsymbol{\eta}_{1} \label{eqn:eta2}\\
&\cdots\nonumber\\
\dot{\boldsymbol{\eta}}_{k} &= [D\mathbf{F}(\mathbf{s}) 
- \alpha D\mathbf{H}(\mathbf{s})]\boldsymbol{\eta}_{k} - \sigma D\mathbf{H}(\mathbf{s})\boldsymbol{\eta}_{k-1}, \label{eqn:eta3}
\end{align}
where $\alpha=\sigma\lambda$ and $\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\ldots,\boldsymbol{\eta}_{k}$ 
%are the modes
%of perturbations along bases of the generalized eigenspace corresponding to the eigenvalue $\lambda$. 
are perturbation modes in the generalized eigenspace of eigenvalue $\lambda$.  

For $\alpha$ regarded as a complex parameter, Eq.~\eqref{eqn:eta1} is 
%called the 
a master stability equation
and its largest Lyapunov exponent $\Lambda(\alpha)$, called 
%the 
master stability function~\cite{pecora1998}, 
determines the stability of Eq.~\eqref{eqn:eta1}: %  against perturbation $\boldsymbol{\eta}_{1}$:
%of the zero solution of Eq.~\eqref{eqn:eta1}:  
it is linearly stable 
%if and only if 
iff $\Lambda(\sigma\lambda) < 0$. 
The condition for Eq.~\eqref{eqn:eta2} to be stable is apparently more involved but
can be 
%argued 
formulated
as follows. The linear stability of Eq.~\eqref{eqn:eta1}
implies that $\boldsymbol{\eta}_1$ converges to zero exponentially as $t \to \infty$.
Assuming that the norm of $D\mathbf{H}(\mathbf{s})$ is bounded,
we have that the second term in Eq.~\eqref{eqn:eta2} is exponentially small.
Then, the same condition $\Lambda(\sigma\lambda) < 0$, now applied to Eq.~\eqref{eqn:eta2}, 
guarantees 
%that 
the stabilizing effect of 
%the first term overcomes the effect of the second
both the first and second terms, 
resulting in exponential convergence of $\boldsymbol{\eta}_2$ to zero as $t \to \infty$.
The same argument applied repeatedly shows that $\boldsymbol{\eta}_3,\ldots,\boldsymbol{\eta}_k$ must also converge to zero if
$\Lambda(\sigma\lambda) < 0$.  This shows
that $\Lambda(\sigma\lambda) < 0$ is a necessary and sufficient condition for the linear
stability of the equations corresponding to each full block $B_j$.
%Jordan block
%%Matrix $L$ always has a zero eigenvalue corresponding to perturbations within the synchronization manifold.
%%Thus, taking all blocks into account, we have that the synchronous solution is stable
%%if and only if $\Lambda(\sigma\lambda) < 0$ for each of the other eigenvalues $\lambda$.
This condition is valid not only in diagonalizable~\cite{pecora1998} but also in
non-diagonalizable networks.



However, it is worthwhile noting a crucial difference between the diagonalizable and
non-diagonalizable cases.  If $L$ is diagonalizable, 
then all Jordan blocks are $1 \times 1$, so there would be no equations like
\eqref{eqn:eta2} or \eqref{eqn:eta3}, and each mode of perturbation 
is decoupled from the others. Thus, the exponential convergence occurs
independently and simultaneously.  
On the other hand, if $L$ is not diagonalizable, some modes
of perturbation may suffer from a long transient.  For instance,
%
if we have a network of linearly coupled phase oscillators, $\dot\theta_i=\omega-\sigma\sum_j L_{ij} \theta_j$, $\theta_i \in S^1$, 
then we can explicitly solve Eqs.~\eqref{eqn:eta1}--\eqref{eqn:eta3}  for the solution $s(t) = \omega t$ to obtain the last perturbation mode $\eta_k= e^{-\alpha t}\sum_{i=0}^{k-1} c_i t^i$, 
where the constants $c_i$ depend on the initial condition.
%if we have a network of phase oscillators, $\dot\theta_i=-\sigma\sum_j L_{ij} \sin\theta_j$, 
%then we can explicitly solve Eqs.~\eqref{eqn:eta1}--\eqref{eqn:eta3}  for the steady state $s(t)\equiv 0$ to obtain the last perturbation mode $\eta_k= e^{-\alpha t}\sum_{i=0}^{k-1} c_i t^i$, 
%where the constants $c_i$ depend on the initial condition.  
Therefore, the larger the size $k$ of the Jordan block, the longer the transient.


Turning our attention back to the maximization problem, we first note that
%We now turn our attention back to the maximization problem. 
%The synchronizability of the network can be expressed solely in terms of 
the eigenvalues $\lambda_1, \ldots ,\lambda_n$ of matrix $L$ can be
%hereafter 
ordered such that $0=\lambda_1 \le\text{Re}\,\lambda_2 \le \ldots \le \text{Re}\,\lambda_n$,
where one eigenvalue is always zero because $L$ has zero row sum and all the others
are guaranteed to have nonnegative real parts because of the Gerschgorin Circle Theorem.
Thus, taking all the Jordan blocks into account,
it follows from our stability analysis that the synchronous solution is stable if and
only if 
%$\Lambda(\sigma\lambda_i) < 0$ for $i=2,\ldots , n$
\begin{equation}
\mbox{$\Lambda(\sigma\lambda_i) < 0$ for $i=2,\ldots , n$}.
\end{equation} 
%($\lambda_1$ corresponds to perturbations parallel to the synchronization manifold).
$\Lambda(\sigma\lambda_1)=\Lambda(0)\ge 0$ is the largest Lyapunov exponent of the
individual oscillators and corresponds to the stability along the synchronization manifold.
We next
%first 
note that Re$\,\lambda_2 > 0$ if and only if the network embeds an
oriented spanning tree, i.e., there is a node from which all other nodes can be reached
by following
% the 
directed links. This condition follows from the recent Ref.~\cite{Wu2005}
and generalizes the notion of connecteness to directed networks. We assume this condition 
here to ensure that the network is compatible with 
synchronization.
%We now turn our attention back to the maximization problem. 
%Note first that the eigenvalues $\lambda_1,\ldots,\lambda_n$ of $L$ are in general complex with nonnegative real part (by the Gerschgorin Circle Theorem) and there is always %a zero eigenvalue (which we denote here by $\lambda_1=0$) because $L$ has zero row sum.  Moreover, the zero eigenvalue is simple if and only if there is a node in the %network from which all
%other nodes can be reached 
%%following the directions of the link~\cite{connected}.  [12/31/2005]
%by following the links~\cite{connected}.  
%We assume this here to ensure that the network is compatible with synchronization.  

In most of the previously studied cases, 
the master stability function $\Lambda(\alpha)$, determined by $\mathbf{F}$, $\mathbf{H}$, and $\mathbf{s}$, 
has been found to be negative in a single convex bounded
region of the complex plane~\cite{conjugate}.
This implies the existence of a single interval 
($\sigma_{\min}$, $\sigma_{\max}$)
of the overall coupling strength $\sigma$ for which synchronization is stable. Thus, the synchronizability
of the network can be measured in terms of the relative interval $\sigma_{\max}/\sigma_{\min}$:
%of the coupling strength:
the network becomes {\it more
synchronizable} as 
%the ratio 
$\sigma_{\max}/\sigma_{\min}$ becomes larger.
%, where $\sigma_{\max}$ and $\sigma_{\min}$ are respectively the maximum and minimum
%values of $\sigma$ for which the synchronous solution is stable.  
In the special case of undirected networks,
the eigenvalues of $L$ are real, and this measure of synchronizability is proportional to the ratio 
$\lambda_2/ \lambda_n$
%: the larger this ratio, the more synchronizable the network~
\cite{barahona2002}.
%In the more general case of possibly directed networks, the network becomes more
%synchronizable as the ratio $\sigma_{\max}/\sigma_{\min}$ becomes larger.

A critical observation is that in order for the ratio $\sigma_{\max} / \sigma_{\min}$ 
to achieve absolute maximum 
%for all choices of $\Lambda(\alpha)$ with a convex stability region, [12/31/2005]
for any given $\Lambda(\alpha)$ with a convex stability region, 
%%we must have all nonzero eigenvalues to be  real and equal to each other.
all nonzero eigenvalues must be  real and equal to each other.
%$\text{Re}\,\lambda_2=\cdots=\text{Re}\,\lambda_n>0$ and 
%$\text{Im}\,\lambda_2=\cdots=\text{Im}\,\lambda_n=0$.  
The condition that the eigenvalues must be real follows from
the convexity of the stability region and the fact that
complex eigenvalues appear in conjugate pairs, 
while the condition that they must be equal follows
from the fact that, for real eigenvalues, the  ratio $\sigma_{\max} / \sigma_{\min}$ 
is proportional to $\lambda_2/\lambda_n$.
%This is so because having one eigenvalue different from others 
%can only make the ratio smaller but not larger.
%The condition that the eigenvalues must be real follows from
%the convexity of the stability region and the fact that complex eigenvalues
%appear in conjugate pairs. 
%For real eigenvalues, we have
%$\sigma_{\max} / \sigma_{\min} \propto \lambda_2/\lambda_n$.
Thus, a network 
%for which
with 
\begin{equation}
0=\lambda_1 < \lambda_2=\cdots=\lambda_n
\label{eq.8}
\end{equation} 
has the {\it widest possible range of coupling strength} in which synchronization is stable,
independently of the individual node dynamics $\mathbf{F}$, output function $\mathbf{H}$, and
synchronous state $\mathbf{s}$, as long as the stability region is convex \cite{non_convex}.


%We show that, for any connected topology of interactions, there
%are assignments of directions and weights for which the resulting network is
%maximally synchronizable and non-diagonalizable. Moreover, u
Under the mild
assumption that the interaction topology
allows no oscillator to interact with all the other oscillators, {\it any maximally 
synchronizable network is necessarily non-diagonalizable}. 
%The latter observation
This
comes from the fact that if $L$ is diagonalizable and satisfies the optimality
condition (\ref{eq.8}) with nonzero eigenvalues equal to $\lambda>0$, then all
the rows of the characteristic matrix $L-\lambda I$ must be equal. In terms of
the network structure, this means that each node must either have uniform output
to all the other nodes
%, and at least one of them must do so,
(at least one of them must do so)
or have no output at all.
These exceptional cases include globally connected networks and directed star
configurations.  
However, it is uncommon in a large complex
network that an oscillator can communicate with all the other oscillators.
%However, it is usually costly to ensure that an oscillator can 
%communicate with all others and this is indeed uncommon in large complex
%networks. 
Therefore, our extension of the master stability analysis to 
non-diagonalizable networks was indeed necessary to properly address the
optimization problem.


%Suppose that a given topology of interaction allows no oscillator to interact with all other oscillators
%, i.e., no node has the full degree.  
%This is a mild assumption that a typical complex network satisfies~\cite{exceptions}.  
%The assumption is also realistic from the design viewpoint, since it may be costly to 
%ensure that an oscillator can communicate with all others.  Under this assumption, 
%one can show~\cite{detail} that for a maximally synchronizable network, $L$ must be non-diagonalizable, 
%i.e., \textit{synchronization is optimal in non-diagonalizable networks}. 
%This comes from the fact that if $L$ is diagonalizable and satisfy the optimality condition above,
%then all the rows of matrix $L-\lambda I$ must be equal. In terms of the network topology,
%this implyes that each node must either have uniform output to all other nodes or have no output at all.  
%The non-diagonalizability of optimal networks means that the extension of the master stability 
%analysis was indeed necessary to properly address the optimization problem.


%In the case of global coupling, where each oscillator can interact with all other oscillators,
%the uniform and undirected assignment of weights leads to this optimum situation.  For such
%assignment we have $L_{ij} = -1$ and $L_{ii} = n-1$ for all $i,j\ne i$, and the eigenvalues
%are $\lambda_1=0$ and $\lambda_2=\cdots=\lambda_n =n$. However, for other topologies of
%interaction, this assignment scheme in general does not lead to identical eigenvalues, 
%indicating that synchronizability is not maximal.

Having observed that optimal networks are rarely diagonalizable,
we now show that, for {\it any connected topology of interactions}, 
there are assignments of directions and weights for which {\it the
resulting network is non-diagonalizable and maximally synchronizable}.
We first note that maximum synchronizability can always be achieved 
by imposing that the network ({\it i}) embeds an oriented spanning tree,
({\it ii}) has no directed loops, and ({\it iii}) has normalized input strengths in
each node, i.e., the total input is the same for all nodes that receive input.
%with input connections.
Condition ({\it i}) guarantees that  Re$\,\lambda_2 > 0$,  condition ({\it ii}) 
guarantees that the eigenvalues are real, and condition ({\it iii}) then
implies the identity (\ref{eq.8}) among the nonzero eigenvalues.
In such optimal networks, we can always rank the nodes so that each 
node receives inputs only from nodes that are higher in the ranking 
(see Fig.~\ref{fig:ost}(a) for an example).  
In this hierarchical 
structure, information flows only from top to bottom 
of the ranking, without feedback.
The optimality can be formally confirmed by noting that indexing 
nodes  according to the ranking makes $L$ a lower triangular matrix 
with  $0,\lambda,\ldots,\lambda$ on the diagonal, which means that 
$\lambda_2=\cdots=\lambda_n=\lambda$, where $\lambda>0$ is the total
input strength in $n-1$ of the nodes. An important class of such 
maximally synchronizable networks consists of the oriented spanning 
trees themselves, where the normalization condition leads to 
uniform weights for all links of the tree (see Fig.~\ref{fig:ost}(b)
for an example).  This example shows that any interaction topology
admits at least $n-1$, but usually many more, optimal non-diagonalizable 
networks. Indeed, from the Matrix-Tree Theorem it follows that the number 
of all oriented spanning trees is $\Pi_{i=2}^{n}\mu_i$, where $\mu_2,\ldots,\mu_n$ 
are the nonzero 
Laplacian eigenvalues of the underlying {\it undirected} network defined 
by matrix $A$. For a globally connected network, for example, the number 
is $n^{n-1}$, which is huge even for
relatively small networks. 
%With the exception of the star configuration, all these oriented spanning trees are non-diagonalizable.
All these oriented spanning trees are non-diagonalizable, except for the star configuration. 
Oriented spanning trees can be explicitly constructed by 
the well-known procedure called the breadth-first search, which spans
all nodes starting from an arbitrary root node.

%FIGURE 2
\begin{figure}[t]
\begin{center}
\epsfig{figure=ost.eps,width=8.5cm}
\caption{(Color online) (a) Example of optimal assignment of weights and directions within a given interaction topology.  
The total input strength in each node is normalized to $\lambda$, where thick, medium, and thin arrows indicate weight $\lambda$,  $2\lambda/3$, and $\lambda/3$, respectively, and dashed lines have zero weight.   
Nodes are numbered and colored to show the hierarchical structure in which connections are only from a higher level to a lower level, with no feedback loops.
%Nodes are ranked so that connections are only from a higher number to a lower number, with no feedback loops.
(b) Example of oriented spanning tree within the same 
interaction topology as in (a), constructed by the breadth-first search.  
%The numbers indicate the order of discovery (which is a possible ranking) and colors represent the levels in the hierarchy.
}
\label{fig:ost}
\end{center}
\end{figure}

Physically, the optimality conditions ({\it i})-({\it iii}) can be understood as follows.  
The top node in the ranking receives no input and acts as a {\it master oscillator} that 
dominates  the network dynamics. If the coupling strength $\sigma$ is chosen so that 
$\Lambda(\sigma\lambda) < 0$, 
then the oscillators that are
immediately lower in the hierarchy and receive input from
the master will synchronize themselves with the master.   Any
oscillator receiving input only from these oscillators and the master
must also synchronize, since normalization of the total input
strength makes the equation effectively look as if it were receiving
input from a single oscillator that is synchronized with the master.  Repeating the same argument 
for the rest of the network, we see that under conditions ({\it i})-({\it iii})
all oscillators must eventually synchronize and they do so for
the entire range of $\sigma$ 
%such that
where $\Lambda(\sigma\lambda) < 0$.  


Interestingly, {\it undirected} tree 
%topology has
networks have
been found to be among the most difficult to 
synchronize~\cite{Yook2005}, in striking contrast to our result that {\em directed} spanning trees 
lead to the most synchronizable
% networks.
configurations.  
  This highlights the significance 
of directionality of the interactions in determining the synchronizability of the networks~\cite{Zheng2000}.  On the other hand, the choice of the master oscillator
% that controls the dynamics of
in a maximally synchronizable network is completely arbitrary, 
despite the intuition that the nodes with largest 
connectivity would be the most natural choice.
%degree (the number of links connected to a node) is the most natural choice for the master.
Moreover, 
%the assignment of direction for each link 
the directions of the links
in such a network 
%does not necessarily correlate with node ordering based on their properties, such as the degree or age, in spite of such suggestion as in~\cite{Hwang2005} that it may be related to the age.
are not necessarily related to the properties of the nodes they connect,
%, such as the age or the degree
even though there has been a suggestion that it would be related to the age of the nodes~\cite{Hwang2005}.
In contrast, under the stricter constraint 
that all feasible input connections
%are required to 
have the same strength in each node,
it was found~\cite{Motter2005b,Motter2005} that maximum synchronizability 
is achieved when 
the individual
input strength is inversely proportional to the 
%degree,
connectivity of the node, 
which is consistent with our result that normalization is key to ensuring optimality.
%In contrast, different constraints were used in \cite{Motter2005b,Motter2005} to conclude that higher synchronizability is achieved when input strength is normalized using the degrees, which is consistent with our results.

The optimality conditions ({\it i})-({\it iii}) suggest that in designing a network for which synchronization is desired, it is generally advantageous
to avoid feedback loops and to normalize input strength.
%In designing a network for which synchronization is desired,
%the optimality conditions ({\it i})-({\it iii}) indicate that it is generally advantageous
%to avoid feedback loops and to normalize input strength.
%into nodes 
Because these conditions typically lead to assigning nonzero weights only to a subset of all 
possible links,
%Because most of such assignments explore only some of the many potential links, 
%our result 
this interesting result
can be interpreted as a synchronization version of the paradox of Braess for traffic flow~\cite{braess},
in which removing links leads counter-intuitively to improved performance of the network.
Furthermore, such assignment of weights not only maximize the synchronizability, but also {\it minimize} the coupling cost.
%it also does it at a fairly low cost.
%assignments satisfying the conditions ({\it i})-({\it iii}) allow a network with
%arbitrary topology of interaction to achieve the same synchronization efficiency as that 
%for globally connected networks. 
%Avoiding directed loops and normalizing total input strength do not only maximize 
%the synchronizability, but it also does it at a fairly low cost.  
%To quantify the efficiency, the total cost can be defined by the sum of the total input strength
The coupling cost can be defined as the sum of the input strengths
of all nodes at the synchronization threshold~\cite{Motter2005b}.  
If $\Lambda(\alpha) < 0 $ in $(\alpha_1, \alpha_2)$, then the coupling cost for any network can only be as small as $\alpha_1(n-1)$, which can be achieved by networks with global uniform coupling.  A surprising fact, however, is that this minimum can also be achieved by the maximally synchronizable networks as well.
In other words, our optimality conditions allow a network constrained
by an arbitrary topology to synchronize with the best possible efficiency.
It is interesting to point out that efficiency optimization of traffic flow on a transportation network model leads to a hierarchical structure similar to that possessed by our maximally synchronizable networks~\cite{Barthelemy2006}.

Our characterization of the maximally synchronizable networks can be used to test the 
widely assumed hypothesis that synchronizability plays an important role in the evolution of many real-world complex networks. The loop structure of
the metabolic network of {\em E. coli} suggests that having fewer loops may have been beneficial for the cell (the details will be published elsewhere), while recent experimental findings~\cite{neuron} suggest the significance of hierarchical structures in neuronal networks.  Exploring more real data to systematically test this hypothesis is of critical fundamental for a better understanding of complex networks.

A. E. M. was supported by DOE under Contract No. W-7405-ENG-36.


\begin{thebibliography}{99}

\bibitem{recent}
Interesting recent works include:
%Self-Sustained Activity in a Small-World Network of Excitable Neurons
A. Roxin, H. Riecke, and S.~A. Solla, Phys. Rev. Lett. {\bf 92}, 198101 (2004);
%Communication and Synchronization in Disconnected Networks with Dynamic Topology: Moving Neighborhood Networks
J.~D. Skufca and E.~M. Bollt,  Math. Biosci. Eng. {\bf 1}, 347 (2004);
%Coherence in scale-free networks of chaotic maps
P.~G. Lind, J.~A.~C. Gallas, and H.~J. Herrmann, Phys. Rev. E {\bf 70}, 056207 (2004);
%Synchronization in power-law networks
L.~Kocarev and P.~Amato, Chaos {\bf 15}, 024101 (2005); 
%Clustering and the synchronization of oscillator networks
P.~N. McGraw and M. Menzinger, Phys. Rev. E {\bf 72}, 015101(R) (2005);
%Onset of synchronization in large networks of coupled oscillators 
J.~G. Restrepo, E. Ott, and B. R. Hunt,  Phys. Rev. E {\bf 71}, 036151 (2005);
%Frequency clustering of coupled phase oscillators on small-world networks,
L.~G. Morelli, H. Cerdeira, and D.~H. Zanette, Eur. Phys. J. B {\bf 43}, 243 (2005);
%
L.~Donetti, P.~Hurtado, and M.~A.~Mu\~noz, Phys. Rev. Lett. {\bf 95}, 188701 (2005).

\bibitem{wattsstrogatz1998}
D.~J. Watts and S.~H. Strogatz, Nature (London) {\bf 393}, 440 (1998).

\bibitem{watts1999}
D.~J. Watts, {\it Small Worlds} (Princeton University  Press, Princeton, 1999).

\bibitem{barabasi1999}
A.-L. Barab\'asi and  R. Albert, Science {\bf 286}, 509 (1999).

\bibitem{Newman2002}
M.~E.~J. Newman, Phys. Rev. Lett. {\bf 89}, 208701 (2002).

\bibitem{Motter2005b}
A.~E. Motter, C.~S. Zhou, and J. Kurths, Europhys. Lett. {\bf 69}, 334 (2005); Phys. Rev. E {\bf 71}, 016116 (2005). 

\bibitem{Motter2005}
A.~E. Motter, C.~S. Zhou, and J. Kurths, AIP Conference Proceedings {\bf 776}, 201 (2005).

\bibitem{Zhou2006}
C.~S. Zhou, A.~E. Motter, and J.~Kurths, Phys. Rev. Lett. {\bf 96}, 034101 (2006).

\bibitem{pecora1998}
L.~M. Pecora and T.~L. Carroll, Phys. Rev. Lett. {\bf 80}, 2109 (1998).

\bibitem{Nishikawa2003}
%T. Nishikawa, A.~E. Motter, Y.-C. Lai, and F.~C. Hoppensteadt, 
T. Nishikawa {\it et al.},
Phys. Rev. Lett. {\bf 91}, 014101 (2003).

\bibitem {Chavez2005}
%M. Chavez, D.-U. Hwang, A. Amann, H.~G.~E. Hentschel, and S. Boccaletti,
M. Chavez {\it et al.}, 
Phys. Rev. Lett. {\bf 94}, 218701 (2005).

\bibitem{fishcer:2005}
E. Fischer and U. Sauer, Nature Genetics {\bf 37}, 636 (2005).

\bibitem{korniss:2003}
G. Korniss {\it et al.},  Science {\bf 299},  677 (2003).

%\bibitem{connected}
%This follows directly from a result in C.~W.  Wu, Linear Algebra and Its Applications {\bf 402},  207 (2005).  
%In the case of undirected networks, the condition is equivalent to the usual notion of connectedness.

\bibitem{Wu2005}
C.~W.  Wu, Linear Algebra and Its Applications {\bf 402},  207 (2005).

\bibitem{conjugate}
In addition, it follows from the complex conjugate of Eq.~\eqref{eqn:eta1}
that $\Lambda(\alpha)$ is always symmetric about the real axis.

\bibitem{barahona2002}
M. Barahona and L.~M. Pecora, Phys. Rev. Lett. {\bf 89}, 054101 (2002).

\bibitem{non_convex}
If the number $n$ of oscillators is even, the same would hold true even
for non-convex regions.

%\bibitem{exceptions}
%Exceptional cases include networks with the global coupling (complete graph) and the star configuration.
 
%\bibitem{detail}
%The detailed proofs will be published elsewhere.

%\bibitem{zero}
%$\Lambda(\sigma\lambda_1)=\Lambda(0)$ is the largest Lyapunov exponent of an individual oscillator, which corresponds to the stability along the synchronization manifold.

\bibitem{Yook2005}
S.-H. Yook and H. Meyer-Ortmanns, e-print cond-mat/0507422.

\bibitem{Zheng2000}
M. Zhan, G. Hu, and J. Yang, Phys. Rev. E {\bf 62}, 2963 (2000);
Z. Zheng, G. Hu, and B. Hu, Phys. Rev. E {\bf 62}, 7501 (2000).

\bibitem{Hwang2005}
D.-U. Hwang {\it et al.}, 
%D.-U. Hwang, M. Chavez, A. Amann, and S. Boccaletti,
Phys. Rev. Lett. {\bf 94}, 138701 (2005).

\bibitem{braess}
D. Braess, Unternehmensforschung {\bf 12}, 258 (1968); A. D. Irvine, Int. Studies in Philosophy of Science 7, {\bf 145} (1993).

%\bibitem{ecoli}
%J.~S. Edwards and B.~O. Palsson, Proc. Natl. Acad. Sci. USA {\bf 97}, 5528 (2000);
%J.~S. Edwards, R.~U. Ibarra, and B.~O. Palsson, Nature Biotechnol. {\bf 19}, 125 (2001).
%%E. Almaas, B. Kovacs, T. Vicsek, Z.~N. Oltvai, and A.-L. Barabasi, Nature (London) {\bf 427}, 839 (2004).

\bibitem{Barthelemy2006}
M.~Barth\'elemy and A.~Flammini, e-print physics/0601203.

\bibitem{neuron}
J.~Jing and K.~R. Weiss.
Current Biology {\bf 15}, 1712 (2005);
R.~Q. Quiroga {\it et al.},
Nature (London) {\bf 435}, 1102 (2005).

%\bibitem{Jing2005}
%J.~Jing and K.~R. Weiss.
%Current Biology {\bf 15}, 1712 (2005).

%\bibitem{Quiroga2005}
%R.~Q. Quiroga {\it et al.},
%%R.~Q. Quiroga, L.~Reddy, G.~Kreiman, C.~Koch, and I.~Fried,
%Nature (London) {\bf 435}, 1102 (2005).

%\end{references}
\end{thebibliography}
\end{document}



cond-mat0605619.extracted_bib
\begin{thebibliography}{99}

\bibitem{recent}
Interesting recent works include:
%Self-Sustained Activity in a Small-World Network of Excitable Neurons
A. Roxin, H. Riecke, and S.~A. Solla, Phys. Rev. Lett. {\bf 92}, 198101 (2004);
%Communication and Synchronization in Disconnected Networks with Dynamic Topology: Moving Neighborhood Networks
J.~D. Skufca and E.~M. Bollt,  Math. Biosci. Eng. {\bf 1}, 347 (2004);
%Coherence in scale-free networks of chaotic maps
P.~G. Lind, J.~A.~C. Gallas, and H.~J. Herrmann, Phys. Rev. E {\bf 70}, 056207 (2004);
%Synchronization in power-law networks
L.~Kocarev and P.~Amato, Chaos {\bf 15}, 024101 (2005); 
%Clustering and the synchronization of oscillator networks
P.~N. McGraw and M. Menzinger, Phys. Rev. E {\bf 72}, 015101(R) (2005);
%Onset of synchronization in large networks of coupled oscillators 
J.~G. Restrepo, E. Ott, and B. R. Hunt,  Phys. Rev. E {\bf 71}, 036151 (2005);
%Frequency clustering of coupled phase oscillators on small-world networks,
L.~G. Morelli, H. Cerdeira, and D.~H. Zanette, Eur. Phys. J. B {\bf 43}, 243 (2005);
%
L.~Donetti, P.~Hurtado, and M.~A.~Mu\~noz, Phys. Rev. Lett. {\bf 95}, 188701 (2005).

\bibitem{wattsstrogatz1998}
D.~J. Watts and S.~H. Strogatz, Nature (London) {\bf 393}, 440 (1998).

\bibitem{watts1999}
D.~J. Watts, {\it Small Worlds} (Princeton University  Press, Princeton, 1999).

\bibitem{barabasi1999}
A.-L. Barab\'asi and  R. Albert, Science {\bf 286}, 509 (1999).

\bibitem{Newman2002}
M.~E.~J. Newman, Phys. Rev. Lett. {\bf 89}, 208701 (2002).

\bibitem{Motter2005b}
A.~E. Motter, C.~S. Zhou, and J. Kurths, Europhys. Lett. {\bf 69}, 334 (2005); Phys. Rev. E {\bf 71}, 016116 (2005). 

\bibitem{Motter2005}
A.~E. Motter, C.~S. Zhou, and J. Kurths, AIP Conference Proceedings {\bf 776}, 201 (2005).

\bibitem{Zhou2006}
C.~S. Zhou, A.~E. Motter, and J.~Kurths, Phys. Rev. Lett. {\bf 96}, 034101 (2006).

\bibitem{pecora1998}
L.~M. Pecora and T.~L. Carroll, Phys. Rev. Lett. {\bf 80}, 2109 (1998).

\bibitem{Nishikawa2003}
%T. Nishikawa, A.~E. Motter, Y.-C. Lai, and F.~C. Hoppensteadt, 
T. Nishikawa {\it et al.},
Phys. Rev. Lett. {\bf 91}, 014101 (2003).

\bibitem {Chavez2005}
%M. Chavez, D.-U. Hwang, A. Amann, H.~G.~E. Hentschel, and S. Boccaletti,
M. Chavez {\it et al.}, 
Phys. Rev. Lett. {\bf 94}, 218701 (2005).

\bibitem{fishcer:2005}
E. Fischer and U. Sauer, Nature Genetics {\bf 37}, 636 (2005).

\bibitem{korniss:2003}
G. Korniss {\it et al.},  Science {\bf 299},  677 (2003).

%\bibitem{connected}
%This follows directly from a result in C.~W.  Wu, Linear Algebra and Its Applications {\bf 402},  207 (2005).  
%In the case of undirected networks, the condition is equivalent to the usual notion of connectedness.

\bibitem{Wu2005}
C.~W.  Wu, Linear Algebra and Its Applications {\bf 402},  207 (2005).

\bibitem{conjugate}
In addition, it follows from the complex conjugate of Eq.~\eqref{eqn:eta1}
that $\Lambda(\alpha)$ is always symmetric about the real axis.

\bibitem{barahona2002}
M. Barahona and L.~M. Pecora, Phys. Rev. Lett. {\bf 89}, 054101 (2002).

\bibitem{non_convex}
If the number $n$ of oscillators is even, the same would hold true even
for non-convex regions.

%\bibitem{exceptions}
%Exceptional cases include networks with the global coupling (complete graph) and the star configuration.
 
%\bibitem{detail}
%The detailed proofs will be published elsewhere.

%\bibitem{zero}
%$\Lambda(\sigma\lambda_1)=\Lambda(0)$ is the largest Lyapunov exponent of an individual oscillator, which corresponds to the stability along the synchronization manifold.

\bibitem{Yook2005}
S.-H. Yook and H. Meyer-Ortmanns, e-print cond-mat/0507422.

\bibitem{Zheng2000}
M. Zhan, G. Hu, and J. Yang, Phys. Rev. E {\bf 62}, 2963 (2000);
Z. Zheng, G. Hu, and B. Hu, Phys. Rev. E {\bf 62}, 7501 (2000).

\bibitem{Hwang2005}
D.-U. Hwang {\it et al.}, 
%D.-U. Hwang, M. Chavez, A. Amann, and S. Boccaletti,
Phys. Rev. Lett. {\bf 94}, 138701 (2005).

\bibitem{braess}
D. Braess, Unternehmensforschung {\bf 12}, 258 (1968); A. D. Irvine, Int. Studies in Philosophy of Science 7, {\bf 145} (1993).

%\bibitem{ecoli}
%J.~S. Edwards and B.~O. Palsson, Proc. Natl. Acad. Sci. USA {\bf 97}, 5528 (2000);
%J.~S. Edwards, R.~U. Ibarra, and B.~O. Palsson, Nature Biotechnol. {\bf 19}, 125 (2001).
%%E. Almaas, B. Kovacs, T. Vicsek, Z.~N. Oltvai, and A.-L. Barabasi, Nature (London) {\bf 427}, 839 (2004).

\bibitem{Barthelemy2006}
M.~Barth\'elemy and A.~Flammini, e-print physics/0601203.

\bibitem{neuron}
J.~Jing and K.~R. Weiss.
Current Biology {\bf 15}, 1712 (2005);
R.~Q. Quiroga {\it et al.},
Nature (London) {\bf 435}, 1102 (2005).

%\bibitem{Jing2005}
%J.~Jing and K.~R. Weiss.
%Current Biology {\bf 15}, 1712 (2005).

%\bibitem{Quiroga2005}
%R.~Q. Quiroga {\it et al.},
%%R.~Q. Quiroga, L.~Reddy, G.~Kreiman, C.~Koch, and I.~Fried,
%Nature (London) {\bf 435}, 1102 (2005).

%\end{references}
\end{thebibliography}