\section{\bf Introduction}
{\noindent \bf Introduction---} The non-linear dynamics of ecological systems where multiple species interact is commonly described by the Lotka-Volterra equations~\cite{Yorke1973} (also known as the predator-prey equations), which contains a set of coupled first-order, non-linear, differential equations. It turns out that the Lotka-Volterra equations finds a broad range of applications beyond ecological systems, and has created a profound impact on physical sciences~\cite{Nutku1990,Kerner1990,Matsuda1992,Malcai2002,Mobilia2007}. In this work, we are interested in a specific class of the Lotka-Volterra equations, where the Malthusian growth/decay term is negligible and the transition coefficients are all anti-symmetric. The resulting equations (see Eq.~(\ref{ALV_eq})) are referred to as anti-symmetric Lotka-Volterra (ALV) equations~\cite{Knebel2015,Reichenbach2006,Vorberg2013,Zakharov1974,DiCera1988}. The ALV equations describe a variety of interesting processes, including predatory-prey oscillations in population biology~\cite{Reichenbach2006}, boson condensation far from equilibrium~\cite{Vorberg2013}, plasma oscillation~\cite{Zakharov1974}, kinetics of chemical reactions~\cite{DiCera1988}, etc. 

Here, we prove that there exists a novel bifurcation mechanism predicted by the ALV equations, formed by flipping the stability of a pair of fixed points (see Fig.~\ref{fig:phase}). Specifically, we study its implications on evolutionary networks as a concrete application; we found that the bifurcation mechanism can be triggered by varying the amount of quantum entanglement~\cite{Horodecki2009a} in the microscopic interaction in the evolutionary networks, which represents a quantum solution to the problem of promoting global cooperation in homogeneous networks.


Over the past few decades, evolutionary network theory~\cite{NM1992} has become the main paradigm of connecting non-equilibrium statistical physics to many other scientific areas including biology, psychology, economics, and behavioral sciences (see e.g. Ref.~\cite{Szabo2007} for a review). 
Real-world networks are realized in systems characterized by graphs in which the vertices (or nodes) represent players, and the edges represent the games played. One of the major goals of evolutionary network theory is to understand the main factors that affect the emergence of {\it global cooperation}. It is known that network topology is crucial for the long-time behaviors of evolutionary networks~\cite{NM1992,HS2005, Szabo2007,Gomez-Gardenes2007, RCS2009}. 
Particularly, it was recognized that {\it heterogeneity} (e.g. existence of large hubs) of networks plays a role in the emergence of cooperations of evolutionary games in networks~\cite{F1974}. 
However, for many {\it homogeneous} networks, which is of interest in this work, it remains a major challenge to achieve global cooperation~\cite{Szabo1998,Wu2006, Nowak2006, Vukov2008}. A central open question in this area is how one can increase the global payoffs of evolutionary games in homogeneous networks. 

The existence of the bifurcation mechanism implies that for a quantum model of evolutionary network, by varying the amount of quantum entanglement~\cite{Horodecki2009a} in the underlying interaction, one can unambiguously achieve the goal of enhancing global cooperation for homogenous networks, avoiding the need of including large hubs.

Moreover, our result implies that the global behaviors of the network can be insensitive to a certain amount of entanglement involved in the microscopic interaction, but they can be abruptly changed, when the entanglement in the microscopic interaction exceeds a critical value. Away from the critical point, we found that the equilibrium state forms a demographic condensations with very weak fluctuations, which resembles Bose-Einstein condensation in statistical physics.

In addition, we show that the non-linear dynamics of the anti-symmetric Lotka-Volterra equation can be constructed by a mean-field theory on the evolutionary networks, which describes the process as a dynamical phase transition in the network; depending on the amount of quantum entanglement involved, the equilibrium state can be driven from an unordered phase to one out of the two condensed phases. 

The origin of the bifurcation is found to be related to a critical scenario where a pair of distinct fixed points exchange their roles from being stable to unstable and vice versa. In our model of evolutionary network, this scenario is possible only if quantum entanglement are provided. Finally, our numerical simulations indicate that the emergence of bifurcation is robust against the change of the network topologies, and the critical values are in good agreement with the theoretical prediction from our mean-field analysis. These results suggest that the bifurcation mechanism could exist in systems other than evolutionary networks; our analysis serves as a guide for finding bifurcation points in other systems governed by the ALV equation.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.7\columnwidth]{fig_phase.pdf}
    \caption{Fixed-point stability and phase diagram obtained by numerical simulation of network of square lattice (see Eq.~(\ref{dynamical_eom_mf})). (a) The stability of two strategies, Defection and Quantum depends on the value of the parameter $\gamma$. (b) The color bar represents the value of the density of Quantum strategy. When $\gamma<\gamma_*$, the network is dominated by the Defection strategy. When $\gamma > \gamma_*$, the network is dominated by the Quantum strategy. The white dotted line comes from the prediction on the phase boundary by our mean-field theory (see Eq.~(\ref{eq:phase_boundary})).}
    \label{fig:phase}
\end{figure}
 

\medskip
{\noindent \bf Anti-symmetric Lotka-Volterra (ALV) equation---} The anti-symmetric Lotka-Volterra (ALV) equation~\cite{Zakharov1974,DiCera1988,Reichenbach2006,Vorberg2013,Knebel2015} is given by the following expression, 
\begin{equation}\label{ALV_eq}
  \frac{d}{{dt}} \ {\rho_i} = {\rho_i} \sum\limits_{j \ne i} {{A_{ij}} \ {\rho_j}} \ ,
\end{equation}
where $\rho_i \geq 0$ is the value of some quantity of interest (e.g. population of a species $i$). The elements of the anti-symmetric matrix, ${A_{ij}} = {w_{ij}} - {w_{ji}} = - A_{ji}$,
contain the transition rates~$w_{ij}$. We first locate the fixed points of the ALV equation: by imposing the normalization condition, $\sum\nolimits_i \ {{\rho _i}}  = 1$, and  setting all $d \rho_i / dt=0$ in the long-time limit, a set of fixed points $\bm {\rho _{{\text{fix}}}} = \{ \ {\rho _i^\infty} \ \} $, where $\rho _i^\infty  \equiv {\rho _i}\left( {t \to \infty } \right)$, can be located by choosing a final state dominated by one of the species, i.e, $\rho _i^\infty  = 1$ and $\rho _{j \ne i}^\infty  = 0$ for all $j \ne i$. It is sufficient to demonstrate our results in the situations involving only three species, i.e., $i=1,2,3$. We shall show that, except for some singular points, in general, there are exactly (i) one stable, (ii) one unstable, and (iii) one saddle fixed points. 

\medskip
{\noindent \bf Fixed-point stability of ALV equation---} Around one of the fixed points, the linearized differential equations is given by 
$  d {\rho _i}/dt = {\bm \nabla} ( {{\rho _i}\sum\nolimits_{j \ne i} {{A_{ij}} \ {{\rho} _j}} } ) \cdot \left( {{\bm \rho}  - {{\bm \rho}_{{\text{fix}}}}} \right)$,
where the gradient, ${\bm \nabla}  \equiv (\partial /\partial {\rho _1},\partial /\partial {\rho _2},...)$, is evaluated at the corresponding fixed point. Furthermore, we consider the corresponding Jacobian matrices, labeled by $ \{ J_{1}, J_{2}, J_{3} \}$ for the three fixed points, namely ${\left[ {1,0,0} \right]^T}$, ${\left[ {0,1,0} \right]^T}$, and ${\left[ {0,0,1} \right]^T}$, respectively, listed as follows:
\begin{equation}
\left \{ \left[ {\begin{array}{*{20}{c}}
  0&{{A_{12}}}&{{A_{13}}} \\ 
  0&{{A_{21}}}&0 \\ 
  0&0&{{A_{31}}} 
\end{array}} \right], \ \left[ {\begin{array}{*{20}{c}}
  {{A_{12}}}&0&0 \\ 
  {{A_{21}}}&0&{{A_{23}}} \\ 
  0&0&{{A_{32}}} 
\end{array}} \right],  \ \left[ {\begin{array}{*{20}{c}}
  0&{{A_{12}}}&{{A_{13}}} \\ 
  0&{{A_{21}}}&0 \\ 
  0&0&{{A_{31}}} 
\end{array}} \right] \right \} .
\end{equation}
The eigenvalue spectra $\lambda \left( {{J_i}} \right)$ of the matrices are found to be $\lambda \left( {{J_i}} \right) \ \in \ \{ \ {0, \ {A_{ji}}, \ {A_{ki}}} \ \} $
for distinct values of $j \ne k \ne i$. Recall that a fixed point is stable (unstable) when the remaining two eigenvalues are both negative (positive); otherwise it is an unstable fixed point. The zero-eigenvalue exist for all Jacobian matrices, which results from the constraint of the normalization condition.

Now, suppose each node (or player) $i$ is associated with a value, $P_i$ (e.g. energy or payoff), where the transition rates $w_{ij}$ are a monotonic function of the difference between the values of $P$'s, i.e., ${w_{ij}} = w( {{P_i} - {P_j}} )$. A common example for the transition rate is given by the Fermi function~\cite{Szabo2007},
\begin{equation}\label{Fermi_function}
{w_{ij}} = 1/{(1 + {e^{ - ({P_j} - {P_i})/T}})} \ ,
\end{equation}
where $T$ is a parameter for controlling the transition rates. We note that if $P_j - P_i>0$, then ${w_{ij}} > {w_{ji}}$ which implies that $A_{ij}$ ($A_{ji}$) is positive (negative), i.e., $A_{ij} > 0 $ and $A_{ji} < 0 $. 

In general, the values of $P$'s are non-degenerate in equilibrium, except for some singular points, which means that there exists an order, e.g., ${P_3} > {P_2} > {P_1}$. In this case, (i) the fixed point ${\left[ {0,0,1} \right]^T}$ is stable, as $A_{13}<0$ and $A_{23}<0$; (ii) ${\left[ {0,1,0} \right]^T}$ is a saddle point, as $A_{12} < 0$ but $A_{32} > 0$; (iii) ${\left[ {1,0,0} \right]^T}$ is unstable, as $A_{21} > 0$ and $A_{31} >0$. Similar results can be obtained by any permutation of the values of $P$'s. 

\medskip
{\noindent \bf Guiding principles for the bifurcation mechanism---} Let us focus on the following scenario: suppose there exist a parameter $\gamma$ (to be identified as the amount of entanglement for evolutionary networks) such that when it is smaller than a critical value~$\gamma_{*}$, i.e., $\gamma < \gamma_{*}$, we have one particular order of the $P$ values, e.g., $P_2 > P_3 > P_1$. Based on our analysis above, it means that the long-time population of the anti-symmetric Lotka-Volterra equation should be dominated by the state $[0,1,0]^T$. Furthermore, suppose, whenever we increase the value of $\gamma_{*}$ to cross the critical point, i.e.,  $\gamma > \gamma_{*}$, we have a different order, e.g., $P_3 > P_2 > P_1$; consequently, we shall be able to observe that the equilibrium state {\it abruptly} changes to another state $[0,0,1]^T$, even though the change the parameter $\gamma$ is smooth.

Instead of an abstract formalism, in the following, we provide a concrete instance on how such a bifurcation mechanism can indeed occur in the context of evolutionary networks, where we prove that the parameter~$\gamma$ can be characterized by the amount of entanglement in the microscopic interactions between the nodes in the network. In principle, for other dynamical systems governed by the anti-symmetric Lotka-Volterra equation, if the corresponding transition rates contain an adjustable parameter that can be varied as described in our model, a bifurcation point can also be located in the same manner, with or without entanglement.


\medskip
{\noindent \bf Evolutionary network as application---} The physical model of the evolutionary network is defined as follows~\cite{Szabo2007,HS2005}. (i) There are $N \gg 1$ rational and identical agents located on the sites in a network; (ii) They interact (formalized as a two-player game) repeatedly with their neighbors to gain/lose an income. (iii) After each game, the agents are allowed to change their strategies in order to increase their utility. In particular, they tend to learn their neighbor's strategies that has generated a higher income. Here a game (see below) is an abstract formulation of an interaction among $N$ players that have potential conflicting interests. The updating probability of a node $i$ to adopt the strategy ${\cal S}_j$ of a reference node $j$, is determined by a transition probability, $w_{ij} \equiv w( {{{\cal S}_i} \to {{\cal S}_j}} )$, which is a function of the difference of the payoffs $P_i - P_j$, and is defined to be identical to the Fermi function (see Eq.~(\ref{Fermi_function})).

Since the local games are played probabilistically, the macroscopic configurations form a probabilistic distribution; we define $Q({\bf n},t)$ to be the probability for the configuration $\bf n$ to appear at time $t$. The rate of change of $Q({\bf n},t)$ is determined by the inflow into and outflow from the configuration $\bf n$, which means that  the dynamics of the evolutionary network can be described by a master equation,
\begin{equation}\label{master_eq}
  \frac{d}{{dt}}Q\left( {{\bf n},t} \right) = \sum\limits_{{\bf n}'} {\left[ {Q\left( {{\bf n}',t} \right){W_{{\bf n}' \to {\bf n}}} - Q\left( {{\bf n},t} \right){W_{{\bf n} \to {\bf n}'}}} \right]} \ ,
\end{equation}
where ${W_{{\bf{n}} \to {\bf{n}}'}}$ is the configurational transition rate from the configuration $\bf n$ to $\bf n '$ and plays the main role for the dynamics.

The configurational transition rate ${W_{{\bf{n}} \to {\bf{n}}'}}$ comes from the microscopic transition rate $W^\mu_{x \to y}{\left( {\bf{n}} \right)}$, which describes an agent changes his strategy from $\mathcal S_x$ to $\mathcal S_y$ for a given macroscopic configuration $\bf n$. The configurational transition rate is proportional to the number $n_x$ of agents with $\mathcal S_x$. Explicitly, we can write ${W_{{\bf{n}} \to {\bf{n}}'}} = \sum\nolimits_{x,y} {{n_x}\;W^{\mu}_{x \to y} {\left( {\bf{n}} \right)} \ {\delta _{{\bf{n}}',{{\bf{n}}_{x \to y}}}}}$,
where $x$ and $y \in \left\{ {a,b,c, \cdots } \right\}$ and the vector ${{\bf{n}}_{x \to y}} = [\cdots, n_x {-}1, \cdots, n_y {+}1, \cdots ]^T$ labels the configuration that is obtaiend from $\bf n$ by changing $n_x \to n_x -1$ and $n_y \to n_y +1$.

\begin{figure}[t]
    \centering
    \includegraphics[width=1\columnwidth]{fig_prisoner.pdf}
    \caption{Microscopic interaction in the evolutionary network. (a) The players are given an entangled pair and three choices \{Cooperate, Defect, Quantum\} for each game; their choice will be sent to a quantum machine to probabilistically make the final decisions for them. (b) The games are played by all members in a network. (c) The payoff table depends on the amount of entanglement (see Eq.~(\ref{eq:PI00011011})).}
    \label{fig:prisoner}
\end{figure}
 


\medskip
{\noindent \bf ALV equation as a mean-field theory---} Furthermore, we can obtain a dynamical equation, $d\left\langle {{n_x}} \right\rangle /dt = \sum\nolimits_{{\bf{n}},{\bf{n}}{'}} {({n'_x} - {n_x}){W_{{\bf{n}} \to {\bf{n}}'}} Q\left( {{\bf{n}},t} \right)} $, for the average occupation of each strategy, where $\left\langle {{n_x}} \right\rangle  \equiv \sum\nolimits_{\bf{n}} {{n_x} \ Q\left( {{\bf{n}},t} \right)} $ (as $\left\langle f \right\rangle  = \sum\nolimits_{\bf{n}} {f\left( {\bf{n}} \right)} Q\left( {{\bf{n}},t} \right)$ in general). In terms of the microscopic rates, we have $d \langle {{n_x}} \rangle /dt  = \sum\nolimits_y {[\langle {{n_y}W^\mu_{y \to x} {\left( {\bf{n}} \right)}} \rangle  - \langle {{n_x} W^\mu_{x \to y} {\left( {\bf{n}} \right)}} \rangle ]}$,
which is so far an exact equation. Assuming the fluctuations around the mean value is small, we can the following approximation: $\langle {n_x} \ W^\mu_{x \to y} {\left( {\bf{n}} \right)}\rangle  \approx \langle {n_x}\rangle \ W^\mu_{x \to y} {\left( \langle {\bf{n}} \rangle \right)}$.
Furthermore, since the players are homogeneous and the transition is taken randomly among the neighboring players, the microscopic transition rate $W_{x \to y}^\mu  ({\bf n}) = {w_{x \to y}} \ {n_y}$ is proportional to the number of players $n_y$ with strategy ${\cal S}_y$. Here ${w_{x \to y}}$ is the transition rate between any pair of players, taken for the mean-field configuration. Therefore, we have $\langle {n_x} W_{x \to y}^\mu \left( {\bf{n}} \right)\rangle  \approx \left\langle {{n_x}} \right\rangle {w_{x \to y}} \langle {{n_y}} \rangle $ under the mean field approximation.

Note that the averaged value $\langle {\bf{n}} \rangle $ is now taken as the input for the microscopic transition rate. Now, we define ${\rho _x}\left( t \right) \equiv \left\langle {{n_x}} \right\rangle /N$ and write, 
\begin{equation}\label{dynamical_eom_mf}
\frac{d}{{dt}}{\rho _x}\left( t \right) = \sum\limits_y {\left[ {{\rho _y}\left( t \right){w_{y \to x}} \ {\rho _x}\left( t \right) - {\rho _x}\left( t \right){w_{x \to y}} \ {\rho _y}\left( t \right)} \right]} \ ,
\end{equation}
which is in exactly the form of the anti-symmetric Lotka-Volterra equation (see Eq.~(\ref{ALV_eq})). The remaining task is to connect this equation with the game and the update rules. 


\begin{figure*}[t]
    \centering
    \includegraphics[width=1.8\columnwidth]{fig_topo.pdf}
    \caption{Numerical simulation on the anti-symmetric Lotka-Volterra (ALV) equation in Eq.~(\ref{dynamical_eom_mf}) for frequencies versus entanglement in different network topologies, namely (a) square lattice (grid of $100 \times 100$), (b) small-world network ($100 \times 100$ network with $p=0.01$), and (c) Erd\H{o}s-R\'{e}nyi model (with average degree $4$). Each data point is taken by playing the games after $10, 000$ steps.}
    \label{fig:topo}
\end{figure*}


\medskip
{\noindent \bf Microscopic interaction---} 
In order to make our discuss concrete, we assume that the rewards of the players are given by the payoff table in the game of prisoner's dilemma, which is one of the most well-known two-player games and is adopted to explain how the emergence of cooperative or altruistic behavior can be resulted by selfish competitors~\cite{Szabo2007}. Specifically, each play has two options, either defection ({\bf D}) or cooperation
({\bf C}), and they have to make the decision simultaneously without communication. If a mutual cooperation is made, both players will be rewarded by an amount $R$. For mutual defection, both of them will be punished with $P$. If one cooperates and the other defects, then the cooperator gains the lowest payoff $S$ and the traitor gains temptation $T$. 

The values of the payoffs are required to have the following order: $T  > R  > P  > S $.
If one of the players (called A) chooses to cooperate, the payoff of the other player (B) is higher if he/she chooses to defect (as $T>R$). The same is true even if A chooses to defect; B still gain more by choosing defection (as $P>S$). If both players are rational, they should both choose defection, i.e., $({\bf D},{\bf D})$, which is the { \it Nash equilibrium}. However, this outcome does not give them the highest average payoffs, which creates the dilemma.



\medskip
{\noindent \bf Making decisions with entanglement---} Entanglement~\cite{Horodecki2009a} is essential for many tasks in quantum information processing e.g.~\cite{Xu2012,Zhang2012c}. In the microscopic interactions of the network, the players are initially shared with an entangled state, $\left| {{\psi _{\rm ini}}} \right\rangle  = J\left| 0 \right\rangle  \otimes \left| 0 \right\rangle$,
where $J = \exp ( {i{\textstyle{\gamma  \over 2}}{\sigma _y} \otimes {\sigma _y}} )$ is an entangling operator (see Fig.~\ref{fig:prisoner}). The degree of entanglement of the initial state is strengthened by increasing the value of the real parameter $\gamma \in [0,2\pi]$; when $\gamma =0$, no entanglement exists in the initial state, but the entanglement reaches its maximum when $\gamma = 2 \pi$. However, the shared entanglement along does not allow the players to communicate. For example, in the protocol of quantum teleportation, a classical communication is needed at the end.

The players in the quantum game~\cite{Eisert1998,Li2014a} are allowed to apply, {\it independently}, a (single-qubit) unitary operation $U$ of the form, ${U}(\theta,\phi) = \left( \begin{smallmatrix} c & s \\  - s & c^* \end{smallmatrix} \right)$, where $c =  e^{i\phi}\cos\theta/2$ and $s= \sin\theta/2$,
on their own qubit. However, we restrict the available choices for the angles are contained in a sets, namely $(\theta ,\phi ) \in \left\{ {(0,0),(\pi ,0),(0,\pi /2)} \right\}$. These three sets of choices correspond to three available strategies, ${\cal S} = \{ {\bf C} , {\bf D} , {\bf Q} \}$; they stand for $\bf C$ooperation, $\bf D$efection, and $\bf Q$uantum, given by $ {\bf C}  \Leftrightarrow  {U}(0,0)$, $ {\bf D}  \Leftrightarrow  {U}(\pi,0)$, and $ {\bf Q}  \Leftrightarrow  {U}(0,\pi/2)$.
After that, the qubits from the players $A$ and $B$ are sent to a quantum machine that applies a unitary operator $J^\dagger$ to the qubits, which yields a final state, $\left| {{\psi _{{\text{fin}}}}} \right\rangle  = {J^\dag } ({U_A} \otimes {U_B}) \left| {{\psi _{{\text{ini}}}}} \right\rangle$.
Finally, a joint quantum measurement, $\left\{ {{\Pi _{00}},{\Pi _{01}},{\Pi _{10}},{\Pi _{11}}} \right\}$, where ${\Pi _{ij}} \equiv \left| {ij} \right\rangle \left\langle {ij} \right|$ for $i,j \in \{ 0, 1\}$, is then applied to the final state to determine the expectation value of the individual utilities or payoffs $p_A$ and $p_B$ of the players, based on the payoff table of the classical game, i.e., 
\begin{equation}\label{eq:PI00011011}
p_{A} = p_B = R\left\langle {{\Pi _{00}}} \right\rangle  + S\left\langle {{\Pi _{01}}} \right\rangle  + T\left\langle {{\Pi _{10}}} \right\rangle  + P\left\langle {{\Pi _{11}}} \right\rangle \ ,
\end{equation}
where $\langle {{\Pi _{ij}}} \rangle  = {\left| {\left\langle {ij} \right.\left| {{\psi _{{\text{fin}}}}} \right\rangle } \right|^2}$, for $i,j \in \{ 0, 1\}$, represents the probability for the outcome $(i,j)$. 




\medskip
{\noindent \bf Fixed-points of evolutionary network---} We denote the population of the three strategies by a vector ${\bf n} = {\left[ {{n_C},{n_D},{n_Q}} \right]^T}$, where ${n_C} + {n_D} + {n_Q} = N$, and ${\rho _X}\left( t \right) \equiv \left\langle {{n_X}} \right\rangle /N$. In the mean-field approach, the transition rate is determined by Eq.~(\ref{Fermi_function}), where the average payoffs $P_X, P_Y \in \{ P_C, P_D, P_Q\}$ are determined by the quantum game of prisoner's dilemma through the following relation:
\begin{equation}\label{Ps_relatedby_rhos}
\left[ {\begin{array}{*{20}{c}}
  {{P_C}} \\
  {{P_D}} \\
  {{P_Q}}
\end{array}} \right] =  \left[ {\begin{array}{*{20}{c}}
  R&S&{{R_\lambda}} \\
  T&P&{{T_\lambda }} \\
  {{R_\lambda}}&{{S_\lambda}}&R
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  {{\rho _C}} \\
  {{\rho _D}} \\
  {{\rho _Q}}
\end{array}} \right] \ ,
\end{equation}
where $R_\lambda \equiv R \lambda + P(1-\lambda)$, $T_\lambda \equiv T\lambda + S (1- \lambda)$, and $S_\lambda \equiv S \lambda + T (1-\lambda) $. The antisymmetric coefficients $A_{XY}$ in the ALV equation is given by ${A_{XY}} = w({S_X} \to {S_Y}) - w({S_Y} \to {S_X}) \equiv \tanh {\Delta _{XY}}$,

Based on the analysis presented previously, the fixed points are given by ${[\rho _C^\infty ,\rho _D^\infty ,\rho _Q^\infty ]^T} = \{ {[0,1,0]^T},{[0,1,0]^T},{[0,0,1]^T}\}$. The stability of these fixed points are determined by the payoff table of the games. Recall that in general we have exactly one stable, one unstable, and one saddle fixed point. In particular, to make ${[0,0,1]^T}$ the stable fixed point, the entanglement parameter has to exceed a critical value, $\gamma  > \gamma_{*1} \equiv {\cos ^{ - 1}}\sqrt {(R - S)/(T - S)}$.

{\it Proof:} To make ${[0,0,1]^T}$ a stable point, we require ${P_Q} - {P_C} > 0$ and ${P_Q} - {P_D} > 0$. The first condition is always satisfied, as ${P_Q} - {P_C} = R - {R_\lambda } = R - R \ {\cos ^2}\gamma  - P \ {\sin ^2}\gamma$ is always positive for the order  $ T> R > P> S$ in the game of prisoner's dilemma. The second condition yields the following inequality, $  R \ {\cos ^2}\gamma  + P \ {\sin ^2}\gamma  < R $, which implies the result. $\blacksquare$ 

Similarly, to stabilize ${\left[ {0,1,0} \right]^T}$, it needs to satisfy $  S{\cos ^2}\gamma  + T{\sin ^2}\gamma  < P $, which means that ${\cos ^2}\gamma  > (T - P)/(T - S)$ or $\gamma  < \gamma_{*2} \equiv {\cos ^{ - 1}}\sqrt {(T - P)/(T - S)}$.
As a result, we now have two critical values, namely $\gamma_{*1}$ and $\gamma_{*2}$. The former stabilizes the fixed point at ${\left[ {0,0,1} \right]^T}$ whenever $\gamma > \gamma_{*1}$, and the latter stabilizes ${\left[ {0,1,0} \right]^T}$ whenever $\gamma < \gamma_{*2}$. If we choose to make the two critical points to coincide with each other, i.e., ${\gamma _{*1}} = {\gamma _{*2}}$. Equivalently, we obtain the critical condition in terms of the payoff values:
\begin{equation}
  T - P = R - S \ .
\end{equation}
As a result, by adjusting the entanglement parameter $\gamma$ across the critical point, we can observe a bifurcation occurs that changes abruptly the equilibrium state from the one dominated by the  defection strategy ({\bf D})  to the one dominated by the quantum strategy ({\bf Q}). To realize such a bifurcation, the payoffs of the game can be assigned as follows: $T = 1 + r$, $R=1$, $P=0$, and $S=-r$, for any $r>0$. In this case, the phase boundary is given by the following relation:
\begin{equation}\label{eq:phase_boundary}
r_* = (1 - {\cos ^2}\gamma_* )/(2{\cos ^2}\gamma_*  - 1) \ ,
\end{equation}
which is represented by the white dotted line in Fig.~\ref{fig:phase}. Our numerical simulation of the ALV equation in Eq.~(\ref{dynamical_eom_mf}) reveals that similar critical transition exists for different network topologies (see Fig.~\ref{fig:topo}), where  the bifurcation points are in good agreement with the mean-field prediction in Eq.~(\ref{eq:phase_boundary}). This completes our analysis on the bifurcation mechanism.

For completeness, let us look at the remaining fixed point at $[1,0,0]^T$. The relevant quantities are ${P_D} - {P_C} = T - R$ and ${P_Q} - {P_C} = {R_\lambda} - R = R{\cos ^2}\gamma  + P{\sin ^2}\gamma  - R$. Since $T>R$, it is always true that $P_D \ge P_C$. Consequently, the fixed point $[1,0,0]^T$ can never become a stable one. To investigate further, since $R>P$, it is also true that $R \ {\cos ^2} \gamma  + P \ {\sin ^2} \gamma  - R < 0$, which implies that $P_C > P_Q$. Therefore, the fixed point at $[1,0,0]^T$ is always a saddle fixed point. 


\medskip
{\noindent \bf Conclusion---}
In summary, we presented a novel bifurcation mechanism for non-linear dynamical systems governed by the anti-symmetric Lotka-Volterra equation. The bifurcation point occurs whenever the stable and unstable fixed points are exchanged by a control parameter. As an application, we focused on evolutionary networks, where the control parameter is characterized by the shared entanglement in the microscopic interaction. Furthermore, we performed numerical simulations to verify that such a bifurcation phenomenon is robust against the change in the network topology. In principle, the bifurcation mechanism can occur in other dynamical systems governed by the ALV equations as well, including predatory-prey oscillations in population biology, condensation of bosons far from equilibrium, plasma oscillation, kinetics of chemical reactions, where the dynamics is governed by the anti-symmetric Lotka-Volterra equation. 

\medskip
{\noindent \bf Acknowledgements---} We thank the support by the Grand Project
``Network Algorithms and Digital Information'' of the Institute of Software, Chinese Academy of Sciences, by an NSFC grant No. 61161130530 and No. 11405093, by a High-Tech Program (863) Grant No. 2012AA8113011, and by a China Basic Science Program (973) Grant No. 2014CB340302.

