\subsection{Gaussian process regression}

Consider a set of $N$ input-output pairs ${\cal D} = \{\vect{x}_i, y_i\}_{i=1}^N$, where $\vect{x}_i \in \Real^D$ are $D$-dimensional input vectors and $y_i\in\Real$ are scalar outputs. We assume that the observed data can be described by the following model,
\begin{equation}
y_i = f(\mathbf{x}_i) + \varepsilon_i,
\label{eq:regmodel}
\end{equation}
in which $f$ represents an unobservable \emph{latent function} and $\varepsilon_i \sim \mathcal{N}(0,\sigma_n^2)$ is zero-mean Gaussian noise.

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution \cite{rasmussen2006gaussian}.
To indicate that a random function $f(\mathbf{x})$ follows a Gaussian process we write it as
\begin{equation*}
f(\mathbf{x}) \sim \GP(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}')).
\end{equation*}
All values of $f$ at any locations $\mathbf{x}$ are jointly normally distributed, with $m(\mathbf{x})$ and $k(\mathbf{x},\mathbf{x}')$ representing the mean function and covariance function, respectively.

In a Bayesian regression setting, we are interested in inferring the predictive distribution of a new, unseen output $y_{\ast}$ given the corresponding input $\mathbf{x}_{\ast}$ and the data ${\cal D}$.
In particular, we take a Gaussian process as the prior over the latent function, and the vector of observations $[y_1,\dots,y_n]^\top$ is related to the latent function through the likelihood function $p(\vect{y}|f)$. When the observations are contaminated with zero-mean Gaussian noise, as in Eq.~\eqref{eq:regmodel},
there exists a closed-form solution for the \emph{posterior distribution} over functions, i.e. the distribution over the unknown function $f(x)$ after incorporating all the observed data. Specifically, the posterior  of the function at any new location $\mathbf{x}_{\ast}$ is described by
\begin{equation*}
p(f_{\ast}|\mathbf{x}_{\ast},{\cal D}) = \Normal(\hat f_\ast,\hat \sigma_\ast^2).
\end{equation*}
When $m(\mathbf{x})=0$, which is a very common assumption, we obtain the following expressions \cite{rasmussen2006gaussian}
\begin{subequations}
\begin{align}
\hat f_\ast &= \vect{k}_\ast^\top (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \vect{y} \label{eq:pred_mean_gpml}\\
\hat \sigma_\ast^2 & = k_{\ast\ast} - \vect{k}_\ast^\top (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \vect{k}_\ast^\top \label{eq:pred_var_gpml}
\end{align}
\label{eq:pred_distr_gpml}%
\end{subequations}
where the covariances (or \emph{kernel}) matrices $\mathbf{K}$ contain the elements $[\mathbf{K}]_{ij} = k(\mathbf{x}_i,\mathbf{x}_j)$ and we have introduced the shorthand notations $\vect{k}_\ast = [k(\mathbf{x}_1,\mathbf{x}_\ast),\dots,k(\mathbf{x}_N,\mathbf{x}_\ast)]^\top$ and $k_{\ast\ast} = k(\mathbf{x}_\ast,\mathbf{x}_\ast)$.
The matrix inversion involved in Eqs.~\eqref{eq:pred_mean_gpml} and \eqref{eq:pred_var_gpml} leads to $\mathcal{O}(N^3)$ complexity.

\subsection{Incremental GP updates}

In an \emph{online} scenario, the data pairs are made available on a one-at-a-time basis, i.e. $(\mathbf{x}_t,y_t)$ arrives at time $t$. Instead of recalculating the predictive distribution entirely once a new data pair $(\mathbf{x}_{t+1},y_{t+1})$ arrives, i.e. by solving \eqref{eq:pred_distr_gpml}, it is more interesting to perform an incremental update.
In this section we briefly review the sequential updates for online GP regression as presented in \cite{vanvaerenbergh2012kernel}. In order to avoid the unbounded growth of the involved matrices, the online learning process is typically coupled with a \emph{sparsification} procedure.

At the $t$-th iteration of the online GP, the model contains the variables
\begin{equation*}
\mathcal{M}_t = \{{\cal D}_t,\vect{\mu}_t,\mat{\Sigma}_t,\mathbf{Q}_t\},
\end{equation*}
where ${\cal D}_t$ is the observed data set that contains the data pairs $\{\mathbf{x}_i,y_i\}_{i=1}^t$; $\vect{\mu}_t$ and $\mat{\Sigma}_t$ are the mean and covariance matrices of the posterior $p(\mathbf{f}_t|{\cal D}_t) = \Normal(\vect{\mu}_t,\mat{\Sigma}_t)$; and $\mathbf{Q}_t = \mathbf{K}_t^{-1}$ is the inverse covariance matrix corresponding to ${\cal D}_t$.

When a new data pair $(\mathbf{x}_{t+1},y_{t+1})$ is obtained, the posterior distribution
\begin{equation*}
p(\mathbf{f}_{t+1}|{\cal D}_{t+1}) = \Normal(\mathbf{f}_{t+1}| \vect{\mu}_{t+1}, \mat{\Sigma}_{t+1})
\end{equation*}
is updated as
\begin{subequations}
\begin{align}
\vect{\mu}_{t+1} & =
\begin{bmatrix} \vect{\mu}_t \\ \hat{y}_{t+1} \end{bmatrix}
+ \frac{y_{t+1} - \hat{y}_{t+1}}{\hat{\sigma}_{yt+1}^2}
\begin{bmatrix} \mathbf{h}_{t+1} \\ \hat{\sigma}_{ft+1}^2, \end{bmatrix}, \label{eq:pmean} \\
\mat{\Sigma}_{t+1} & =
\begin{bmatrix} \mat{\Sigma}_t & \mathbf{h}_{t+1}\\ \mathbf{h}_{t+1}^\top & \hat{\sigma}_{ft+1}^2 \end{bmatrix}
- \frac{1}{\hat{\sigma}_{yt+1}^2}
\begin{bmatrix} \mathbf{h}_{t+1} \\ \hat{\sigma}_{ft+1}^2 \end{bmatrix}
\begin{bmatrix} \mathbf{h}_{t+1} \\ \hat{\sigma}_{ft+1}^2 \end{bmatrix}^\top, \label{eq:pcov}
\end{align}
\end{subequations}
where
\begin{align*}
\mathbf{q}_{t+1} & = \mathbf{Q}_t \vect{k}_{t+1}, \\
\mathbf{h}_{t+1} & = \mat{\Sigma}_t \mathbf{q}_{t+1},
\end{align*}
and the vector $\vect{k}_{t+1}$ has elements $[\vect{k}_{t+1}]_i = k(\mathbf{x}_i,\mathbf{x}_{t+1})$. Furthermore, the output variance is obtained as
\begin{equation*}
\hat{\sigma}_{yt+1}^2 = \sigma_n^2 + \hat{\sigma}_{ft+1}^2,
\label{eq:sigma_y}
\end{equation*}
and the variance of the latent function evaluations is
\begin{align*}
\hat{\sigma}_{ft+1}^2 & = k_{t+1} + \vect{k}_{t+1}^\top (\mathbf{Q}_t\mat{\Sigma}_t\mathbf{Q}_t-\mathbf{Q}_t)\vect{k}_{t+1}\\
& = \gamma_{t+1}^2 + \mathbf{q}_{t+1}^T \mathbf{h}_{t+1}.
\end{align*}
In order to further simplify equations, the variable
\begin{equation*}
\gamma_{t+1}^2 = k_{t+1} - \vect{k}_{t+1}^T \mathbf{Q}_t \vect{k}_{t+1}
\end{equation*}
is introduced.
The inverse kernel matrix $\mathbf{Q}_t$ can be updated efficiently through
\begin{equation}
\mathbf{Q}_{t+1} =
\begin{bmatrix}
\mathbf{Q}_t & \mathbf{0} \\ \mathbf{0}^\top & 0
\end{bmatrix}
+ \frac{1}{\gamma_{t+1}^2}
\begin{bmatrix}
\mathbf{q}_{t+1} \\ -1
\end{bmatrix}
\begin{bmatrix}
\mathbf{q}_{t+1} \\ -1
\end{bmatrix}^\top.
\label{eq:q_update}
\end{equation}
After applying Eqs.~\eqref{eq:pmean}, \eqref{eq:pcov}, and \eqref{eq:q_update}, we obtain the updated model $\mathcal{M}_{t+1} = \{{\cal D}_{t+1},\vect{\mu}_{t+1},\mat{\Sigma}_{t+1},\mathbf{Q}_{t+1}\}$.

At each step $t$ of the learning process, the predictive distribution of a new observation $y_{t+1}$ given all past data is a Gaussian
$
p(y_{t+1}|{\cal D}_t) = \Normal(y_{t+1}|\hat y_{t+1},\hat \sigma_{yt+1}^2)
$
with
\begin{subequations}
\begin{align}
\hat y_{t+1} & = \mathbf{q}_{t+1}^\top \vect{\mu}_t = \vect{k}_{t+1}^\top \mathbf{Q}_t \vect{\mu}_t \label{eq:pred_mean_ogp} \\
\hat \sigma_{yt+1}^2 & = \sigma_n^2 + k_{t+1} +
\vect{k}_{t+1}^\top (\mathbf{Q}_t\mat{\Sigma}_t\mathbf{Q}_t-\mathbf{Q}_t)\vect{k}_{t+1}. \label{eq:pred_cov_ogp}
\end{align}
\label{eq:pred_distr_ogp}
\end{subequations}
For more details we refer the reader to \cite{vanvaerenbergh2012kernel,lazaro2011bayesian}.

\subsection{Kernel recursive least-squares}
The kernel-methods counterpart of online GP regression is kernel recursive least-squares (KRLS, see for instance \cite{engel2004kernel}), which obtains the KRR solution recursively. After receiving $t$ data points, the kernel weights obtained by KRLS are those that solve the batch problem
\begin{equation}
\vect{\alpha}_t = (\mathbf{K}_t + \sigma_n^2\mat{I})^{-1}\vect{y}_t.
\label{eq:krls_alpha}
\end{equation}
The same weights can be obtained in online GP regression by computing
\begin{equation}
\vect{\alpha}_t = \mathbf{K}_t^{-1}\vect{\mu}_t = \mathbf{Q}_t\vect{\mu}_t,
\end{equation}
which follows from Eq.~\eqref{eq:pred_mean_ogp}.

\subsection{KLMS algorithms}

Similar to online GP regression, updating the KRLS estimate in Eq.~\eqref{eq:krls_alpha} with a new data point requires quadratic complexity. Kernel least-mean-squares (KLMS) algorithms alleviate this computational burden by performing stochastic gradient descent of the mean square error, resulting in linear complexity per time step \cite{vanvaerenbergh2014online}.

KLMS algorithms can be categorized into two classes, depending on which kernel weights they update in each iteration to account for the prediction error\footnote{Both approaches are different approximations of the same update formulation in the \emph{kernel feature space}; see \cite{vanvaerenbergh2014online} for details.}. We outline both approaches briefly in the remainder of this section.

\subsubsection{Type-I KLMS: concentrating the novelty}

The first type of KLMS algorithm updates only one coefficient in order to compensate for the prediction error, at each time step. When the weight vector is allowed to grow, this update takes the form
\begin{equation}
\vect{\alpha}_{t+1}^{(I)} = \begin{bmatrix}
\vect{\alpha}_t\\ \eta e_{t+1}
\end{bmatrix},
\label{eq:klms_update}
\end{equation}
in which $e_{t+1} = y_{t+1} - \hat y_{t+1}$ is the instantaneous error.

Eq.~\eqref{eq:klms_update} represents the basic update of what is known as the KLMS algorithm, proposed in \cite{liu2008kernel}. In order to avoid the infinite growth of $\vect{\alpha}_t$, a more sophisticated version of this algorithm was presented in \cite{chen2012quantized}, known as Quantized Kernel Least Mean Square (QKLMS). When QKLMS receives a datum similar to a previously seen datum, for instance the $i$-th base it has stored, it does not expand $\vect{\alpha}_t$ but instead updates the corresponding weight $\alpha_i$.

\subsubsection{Type-II KLMS: spreading the novelty}

A different strategy consists in updating all coefficients of $\vect{\alpha}_t$ in each iteration. This approach is followed for instance by the Kernel Normalized Least Mean Square (KNLMS) algorithm \cite{richard2009online}, whose update reads
\begin{equation}
\vect{\alpha}_{t+1}^{(II)} = \begin{bmatrix}
\vect{\alpha}_t\\ 0
\end{bmatrix} + \eta \frac{e_{t+1}}{\epsilon + k_{t+1}^2 + \|\vect{k}_{t+1}\|^2} \begin{bmatrix}
\vect{k}_{t+1}\\ k_{t+1}
\end{bmatrix}
\label{eq:knlms_update}
\end{equation}
when $\vect{\alpha}_t$ is allowed to grow. Note that the rule \eqref{eq:knlms_update} updates all coefficients in each iteration. In order to avoid unbounded growth, KNLMS follows a \emph{coherence} criterion that promotes sparsity.

\subsection{Specific cases: $\beta=0$ and $\beta=1$}

Setting $\beta=0$ in Eq.~\eqref{eq:update_alpha_subs} yields the update rule
\begin{equation}
\vect{\alpha}_{t+1}^{\beta=0} = \begin{bmatrix}
\vect{\alpha}_t\\ \frac{1}{\sigma_n^2 + k_{t+1}} e_{t+1}
\end{bmatrix}.
\label{eq:update_alpha_subs1}
\end{equation}
This rule is very similar to the KLMS update \eqref{eq:klms_update}, and even identical when the learning rate is set to $\eta = 1/(\sigma_n^2 + k_{t+1})$. Note, furthermore, that $\beta=0$ indicates that the posterior covariance from Eq.~\eqref{eq:pcov_simple} reduces to
\begin{equation}
\mat{\Sigma}_t^{\beta=0} = \mathbf{K}_t.
\end{equation}
In other words, for $\beta=0$ we obtain a KLMS model that implies a fixed posterior covariance, equal to the prior covariance.
The predictive variance, Eq.~\eqref{eq:beta_sigma}, simplifies to
\begin{equation}
\hat\sigma^2_{ft+1}|^{\beta=0} = k_{t+1}.
\label{eq:pred_cov_beta0}
\end{equation}
Under the adopted framework, this is the model that underlies algorithms such as KLMS from \cite{liu2008kernel} and QKLMS from \cite{chen2012quantized}.

By setting $\beta=1$ in Eq.~\eqref{eq:update_alpha_subs}, the update rule becomes
\begin{equation}
\vect{\alpha}_{t+1}^{\beta=1} = \begin{bmatrix}
\vect{\alpha}_t\\ 0
\end{bmatrix} + \frac{e_{t+1}}{\sigma_n^2 + k_{t+1} + \|\vect{k}_{t+1}\|^2} \begin{bmatrix}
\vect{k}_{t+1}\\ 1
\end{bmatrix}
\label{eq:update_alpha_subs2}
\end{equation}
which is very similar to the update of KNLMS, shown in Eq.~\eqref{eq:knlms_update}. The implied posterior covariance then reads
\begin{equation}
\mat{\Sigma}_t^{\beta=1} = \mathbf{K}_t \mathbf{K}_t + \mathbf{K}_t,
\end{equation}
and the predictive covariance from Eq.~\eqref{eq:beta_sigma} now becomes
\begin{equation}
\hat\sigma^2_{ft+1}|^{\beta=1} = k_{t+1} + \vect{k}_{t+1}^\top \vect{k}_{t+1}.
\label{eq:pred_cov_beta1}
\end{equation}
Interestingly, this KLMS model implies a predictive covariance, for each point in space, that is larger than the prior covariance $k_{t+1}$. Furthermore, a closer inspection of the second term in Eq.~\eqref{eq:pred_cov_beta1} shows that the predictive variance grows as more training data is processed.
According to the GP framework, this is the model that underlies type-II KLMS algorithms, in particular KNLMS from \cite{richard2009online}, for which $\beta = 1$, and in general any algorithm that corresponds to $\beta>0$.

As seen through the adopted probabilistic perspective, the capability of Type-II KLMS algorithms to update all coefficients with each new data point is related to an increase in prediction uncertainty.
These effects are illustrated in Fig. \ref{fig:uncertainty}, which compares the interpreted predictive uncertainty for the different algorithms. GP regression exhibits the expected behavior, i.e. uncertainty shrinks around the observed data points. Type-I KLMS ($\beta=0$) assumes a fixed predictive covariance, regardless of the number of observed data and the distances between them. The behavior of type-II KLMS (for instance with $\beta=1$) is rather counterintuitive: its update procedure increases the predictive variance as more data points are observed, and this happens in the neighborhoods of those points. While an in-depth study is needed to extract solid conclusions from this observation, it is interesting to note that a similar behavior has been observed in the linear recursive least squares (RLS) algorithm with forgetting factor, see \cite{vanvaerenbergh2012kernel}.

\begin{figure}[t!]
\centering
\includegraphics[width=\linewidth]{uncertainty_plot_3}
\includegraphics[width=\linewidth]{uncertainty_plot_8}
\includegraphics[width=\linewidth]{uncertainty_plot_25}
\caption{Comparison of the predictive variance of three algorithms, for $3$ data points (top plot, data is marked as blue crosses), $8$ data points (middle), and $25$ data points (bottom). The predictive mean of GP regression is indicated as the black curve that passes through the observations. The grey zone marks the GP mean plus/minus two standard deviations $\hat \sigma_{y}$, corresponding to the GP's $95\%$ confidence interval. The dashed line marks the confidence interval for type-I KLMS ($\beta=0$), and the dash-dot line marks the confidence interval for type-II KLMS with $\beta=1$.
}
\label{fig:uncertainty}
\end{figure}

\subsection{Online regression on stationary data}

In the first experiment, we wish to study the effect of different values of $\beta$ in the KLMS algorithm from Eq.~\eqref{eq:update_alpha_subs}, which we will denote by $\beta$-KLMS. We evaluate this algorithm and three established kernel adaptive filtering algorithms on the stationary KIN40K benchmark.\footnote{Available at \url{http://www.cs.toronto.edu/~delve/data/datasets.html}} This data set is obtained from the forward kinematics of an $8$-link all-revolute robot arm, and it represents a very difficult regression problem. We randomly select $5000$ data points for online training, and $5000$ points for testing the regression.

The algorithms are considered in their evergrowing version here, in order to highlight only the influence of $\beta$. The KRLS-T algorithm, however, which implements the full online GP regression, is given a limited memory of $500$ bases, for computational reasons. A Gaussian kernel was used for all algorithms, with parameters determined offline by standard GP regression, as detailed in \cite{vanvaerenbergh2012kernel}.

The results are shown in Fig.~\ref{fig:kin40k}. Each point of the learning curves corresponds to the test error on the entire test set. We observe that for $\beta=1$, the $\beta$-KLMS algorithm obtains similar performance to KNLMS from \cite{richard2009online}. For $\beta=0$, the performance is almost identical to that of KLMS from \cite{liu2008kernel}.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{kin40k_nmse}
\caption{Performance comparison for online prediction on the KIN40K benchmark regression problem.} \label{fig:kin40k}
\end{figure}

\subsection{Online prediction of time series with a model switch}

In Fig. \ref{fig:reconvergence} we repeat a standard experiment from the kernel adaptive filtering literature. In particular, we consider a nonlinear system comprised of a linear channel followed by a nonlinearity, and the online regression task consists in predicting the next sample of the time series produced at its output. An abrupt change in the linear channel is triggered after $500$ time steps, in order to test the algorithms' ability to reconverge. The experiment is repeated $5$ times with random channel coefficients, and the average results are shown in Fig. \ref{fig:reconvergence}. The same kernel was used for each algorithm.

In this experiment we observe that, as is usual in the LMS and KLMS literature, the $\beta$-KLMS model exhibits tracking behavior although it is based on a static data model. Furthermore, it does so without the need of including an additional parameter to control the update step size.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{reconvergence}
\caption{Performance comparison for online prediction on a nonlinear channel with an abrupt change.} \label{fig:reconvergence}
\end{figure}

\section{Introduction}
\label{sec:intro}
\input{1_introduction.tex}

\section{Online GP regression}
\label{sec:ogp}
\input{2_ogp.tex}

\section{Kernel adaptive filtering and KLMS}
\label{sec:klms}
\input{3_klms.tex}

\section{From online GP to KLMS}
\label{sec:approx}
\input{4_approximation.tex}

\section{Experiments}
\label{sec:exp}
\input{6_experiments.tex}

\section{Conclusions}
\label{sec:concl}
\input{7_conclusions.tex}

\bibliographystyle{IEEEbib}


\section{Introduction}
\label{intro}

Uniaxially anisotropic antiferromagnets in a
magnetic field have been studied quite extensively in the
past, both experimentally
and theoretically. Typically, they display, at low
temperatures, the antiferromagnetic (AF) phase and, when increasing
the field, the spin--flop (SF) phase \cite{Neel,Gorter}. In addition, more
complicated structures, like biconical \cite{KNF} ones, have been
observed. The various ordered phases may
lead to interesting multicritical
behavior, including bi-- and tetracritical points \cite{KNF,LF,FN,Aha2}.

Experimentally, several antiferromagnets with uniaxial anisotropy
have been investigated, three--dimensional magnets
\cite{SF,RG,FK,OU,FPB}, such as  MnF$_2$ and GdAlO$_3$, as
well as  quasi two--dimensional magnets \cite{SCD,T,BS,BSL}.

Much of the theoretical work is based on analyzing the prototypical
classical Heisenberg model with uniaxial anisotropy, the XXZ
model, in a field, plus, possibly, further anisotropy
terms, such as single--ion anisotropies. Especially, mean--field
approximation \cite{Gorter,LF}, Monte Carlo
simulations \cite{LanBin,Holt1,hws,Zhou2,hs,sel,htl,fp}, and renormalization
group calculations \cite{KNF,FN,CPV,Folk,Eich} have been applied.

In this article, we shall deal with the anisotropic
XY antiferromagnet in a field. The $x$--axis is taken to be the
easy axis and the field acts on the $x$--component
of the spins. The model is a variant of the
much studied uniaxially anisotropic three--component XXZ antiferromagnet
in a field, with the field along the easy axis, the $z$--direction. In a
previous paper, the two--dimensional version, on the square lattice, of
the XY model had been studied, applying ground state considerations
and Monte Carlo techniques \cite{hs}. At the field
separating the AF and SF phases, the
ground state has been found to be highly degenerate due to the
presence of biconical (or, more precisely, non--collinear,
biangular or bidirectional, BD) structures. At finite
temperatures, $T >0$, these structures seem to
lead to a disordered phase between the AF and SF phases, similar
to the situation in the two--dimensional XXZ Heisenberg
antiferromagnet. Here, we shall briefly
reconsider this case. Our main
emphasis will be on the anisotropic XY model on the simple cubic
lattice. To our knowledge, no prior analysis exists. In particular, the
existence and nature of the possible bicritical point, at which
the AF, SF, and paramagnetic phases are expected to merge, in analogy to
the XXZ antiferromagnet, will be studied.

The paper is organised as follows: In the next section, the
anisotropic XY model will be introduced and ground state properties will
be discussed. Then, results for the model on the square lattice will
be presented, followed by our large--scale simulation
findings for the three-dimensional case. A short summary concludes
the article.

\section{The model}
The anisotropic XY antiferromagnet in a field is described by the
Hamiltonian

\begin{equation}
{\cal H}_{\mathrm{XY}} = J \sum\limits_{i,j} \left[ \, S_i^x S_j^x + \Delta S_i^y S_j^y \, \right] \; - \; H \sum\limits_{i} S_i^x
\end{equation}

\noindent
where the first sum runs over all pairs of neighbouring
sites, $i$ and $j$, of the lattice, with the second sum running
over all lattice sites. $J > 0$ is the antiferromagnetic coupling
constant, $\Delta <1$ is the anisotropy parameter, and $H$ is
the external field along the easy
axis, the $x$-axis. $S_i^{\alpha}$, $\alpha$=$x$ and $y$,
are the two components of the classical spin vector of length one.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov.eps}
}
\caption{Ground state configurations: (a) antiferromagnetic (AF),
 (b,d) spin--flop (SF), and (c) biconical as well as (e)
 bidirectional (BD) configurations. (a,d,e) occur in the anisotropic
 XY, (a,b,c) in the uniaxially anisotropic XXZ antiferromagnet. The
 tilt angles are defined with respect to
 the easy axis, being the $x$--axis in the XY and the $z$--axis in
 the XXZ case.}
\label{fig:1}
\end{figure}

We shall consider mainly the model on a simple cubic lattice ($d=3$), but
some results for the square lattice ($d=2$) will also be discussed
comparing them to previous findings \cite{hs}. In both cases, the
ground states, at zero temperature, $T=0$, see Figure 1, may be
determined in a straightforward way \cite{LF,hs,TM}. For small
fields, $H < H_{c1}$, the AF
structure is stable, in which neighbouring sites belong to
different sublattices, $A$ and $B$, where the $x$--components
of $A$-- and $B$--spins point in opposite directions. The $y$--component
vanishes at all sites. At intermediate
fields, $H_{c1} <H< H_{c2}$, one encounters
the SF structure, where the $x$--components of the spins on
both sublattices are equal, being smaller than one, while the
antiferromagnetic sublattice
structure shows up in the $y$--components.  At large
fields, $H > H_{c2}$, all spins
are aligned parallel to the field, with vanishing $y$--component
of the spin vectors. At $H=H_{c1}$, the
AF, the SF, and BD configurations form ground states. In
the highly degenerate BD configurations, the spins on the two
sublattices, $A$ and $B$, are tilted with respect to the field
direction, the $x$--axis. The resulting two tilt angles, say, $\Theta_A$ and
$\Theta_B$, are interrelated interpolating
continuously between the AF and SF structures. The exact relation
depends on the anisotropy parameter $\Delta$. Of
course, in the SF limit, one
has $\Theta_A= \Theta_B= \Theta_{SF}(H) >0$.

At non-zero temperatures, one expects, among others, phase
transitions in the Ising universality class from the
disordered phase to the AF and SF phases. The
longitudinal staggered magnetization, describing the antiferromagnetic
ordering of the $x$--component of the spins on the
$A$-- and $B$--sublattices, is
the order parameter in the AF phase. The transversal
staggered magnetization describes the
antiferromagnetic ordering of the $y$--component of the
spins in the SF phase, in which the $x$--component has the same value on
both sublattices. Accordingly, a possible bicritical
point is expected to belong to the $n=2$ or XY universality
class \cite{KNF,FN}. This behaviour
is in marked contrast to
the situation in the anisotropic XXZ antiferromagnet in a
field, along the $z$--axis. There
the SF phase is described by an antiferromagnetic
ordering of the two components perpendicular to the
$z$--component of the spins, implying a transition to the
paramagnetic phase
in the XY universality class. Of course, the transition from
the paramagnetic to the AF phase belongs to the Ising
universality class, in the XXZ case as well. Thence, a
possible bicritical point would fall into the $n=3$ or
Heisenberg universality class \cite{KNF,FN}.

To study the phase diagram, in the $(k_BT/J,H/J)$ plane, of the
2d and 3d XY antiferromagnets, we did extensive Monte Carlo (MC)
simulations, using the standard Metropolis algorithm \cite{LB}. Lattices
with $L^d$ sites were considered, with $L$ ranging from 10 to 200 for
$d=2$, and from 8 to 40 for $d=3$. In all cases, full periodic
boundary conditions were employed. As usual, finite--size extrapolations
were done to obtain estimates for the thermodynamic limit. Typically, runs
of, at least, $10^7$ MC steps per site (MCS) were performed, averaging
over a few realizations, by using different random numbers, to estimate
error bars. Here, error bars are usually smaller than the size of
the symbols shown in the figures. Of course, close to the
phase transition, larger lattices may
require longer runs to take into account critical fluctuations and
critical slowing down \cite{LB}.

We recorded, among others, quantities related to
the longitudinal staggered magnetisation $M_{st}^x$
defined by $M_{st}^x= (S_A^x-S_B^x)/2$, and the corresponding
transversal staggered magnetization for the $y$ component of the
spins. In particular, thermal averages over the second moments
of the magnetizations are expected to signal the transitions to
the AF and SF phases. Actually, these transitions
may be detected quite easily and reliably by the
Binder cumulant \cite{Binder}

\begin{equation}
  U^{x,y} = 1 - <(M_{st}^{x,y})^4>/(3 <(M_{st}^{x,y})^2>)
\end{equation}

\noindent
where the brackets denote thermal averages. The transition
from the disordered phase to the AF phase can be determined by
$U^x$, and that to the SF phase can be determined
by $U^y$. To detect a possible bicritical point in
the $n=2$ or XY universality class, we also computed the
analoguous Binder cumulant $U^{xy}$, invoking the second and
fourth moments of the total staggered magnetization
$M_{st}= \sqrt{((M_{st}^x)^2 + (M_{st}^y)^2)}$.

We also recorded the longitudinal and transversal
staggered susceptibilties, obtained, as usual, from the fluctuations
of the corresponding magnetizations. Further interesting information on
the thermal behavior of the model follows from the specific
heat, $C$, determined from the energy fluctuations, as
well as histograms. In particular, we computed probability functions
of the tilt angles, such as the probability $p(\Theta)$ for encountering
the tilt angle $\Theta$ at an arbitrary site and the
probability $p_2(\Theta_A,\Theta_B)$ for finding the two
angles, $\Theta_A$ and $\Theta_B$, at
neighbouring sites of the lattice, as before \cite{hs}.

\section{Phase diagrams}

We shall first briefly discuss previous and new MC findings for the 2d
case, presenting evidence for a narrow disordered phase separating
the AF and SF phases. The main emphasis will be on the 3d model, for
which a qualitatively different topology of the phase diagram is
observed. In fact, our simulation data suggest the existence
of a bicritical point, at which the AF, SF, and paramagnetic phases
meet. In both cases, we set the anisotropy parameter $\Delta = 0.8$.

\subsection{Anisotropic XY antiferromagnet on a square lattice}

The phase diagram of the anisotropic XY antiferromagnet in a field on
a square lattice has been determined before \cite{hs}, comprising
the AF, SF, and paramagnetic phases. It is depicted in Figure
2, setting the anisotropy parameter $\Delta$ equal to 0.8. The
characteristic fields for the ground states are $H_{c1}/J = 2.4$
and $H_{c2}/J = 7.2$. The transition lines of the AF and SF phases
to the disordered phase, having confirmed to be in the
Ising universality class, approach each other rather closely at
$k_BT/J \approx 0.68$ and $H/J \approx 2.47$.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov1.eps}
}
\caption{Phase diagram of the anisotropic XY antiferromagnet
   on the square lattice with $\Delta=0.8$. From Reference 20.}
\label{fig:2}
\end{figure}

In principle, the two lines might meet at a bicritical point in
the XY or Kosterlitz--Thouless \cite{KT} universality class, with
a line of transitions of first order between the AF and SF phases
at lower temperatures. Note that such a topology is excluded for the
two--dimensional XXZ antiferromagnet, where the bicritical
point would belong to the Heisenberg universality class. In
fact, the existence of such a point at $T>0$ is
excluded by the well--known theorem of Mermin and Wagner \cite{MW}.

In our present study, we searched for possible evidence, whether
there is a bicritical point in the anisotropic XY model
on a square lattice. In particular, we did MC simulations at the
temperatures $k_BT/J$ = 0.7, 0.6, and 0.4.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov2.eps}
}
\caption{Finite--size dependence of the crossing points
 of the longitudinal Binder cumulants $U_{cr}^x(L_{eff})$
 at $k_BT/J= 0.6$ (squares) and $0.7$ (circles). The value of the
 critical cumulant for the
 Ising case $U^*=0.61069...$ is indicated.}
\label{fig:3}
\end{figure}

The Binder cumulants $U^{x,y}$ turned out to be very useful in
investigating the low--temperature
region. We monitored the crossing points
$U_{cr}(L_{eff})$ of the
cumulants for two successive square lattices with linear dimensions
$L_1$ and $L_2$, assigning an effective length
$L_{eff}= (L_1+L_2)/2$. In general \cite{Binder}, one expects
$U_{cr}$ to occur in the thermodynamic limit at the phase
transition, with the critical value $U^*$. In the
case of phase transitions in
the 2d Ising universality class, employing full periodic
boundary conditions and considering systems with spatially
isotropic spin--spin correlations \cite{cd,sesh}, one has
$U^*$= 0.61069 ..\cite{sesh,kabl}.

In Figure 3, results for the crossing points of the
longitudinal Binder cumulant $U^x$ at $k_BT/J$ = 0.6
and 0.7, with lattice sizes ranging from $L$= 20 to 200, are
depicted. The typical critical value, $U^*$= 0.61069..,
seems to be approached for both temperatures, with smaller
finite--size effects at the higher
temperature. Similar observations hold for the transversal
cumulant. Accordingly, the two transitions are seen to
be Ising--like. Note that at $k_BT/J$= 0.6, the two
distinct transitions occur at about $H/J$ = 2.4505
and 2.4525, as estimated from
the size--dependence of the crossing points, $U_{cr}^x$ and
$U_{cr}^y$. The disordered phase between the AF and SF
phases becomes more and more narrow, as the temperature
is lowered, see also Figure 2.

The existence of that narrow disordered phase had been
previously \cite{hs} inferred from histograms of the
probability $p_2$. $p_2$ describes the
relation between the tilt angles at neighbouring
sites, as stated above. Actually, at the temperature $k_BT/J = 0.4$, in an
extremely small range of fields, the
dominant configurations have neither an AF nor SF character, but
they are of BD type. The disordered phase had been argued to
be due to BD fluctuations caused by the highly degenerate
ground state at $H_{c1}$.

Indeed, the present simulations provide additional evidence
against a transition of first order between the AF and
SF phases at such low temperatures. For example, the
longitudinal and transversal Binder cumulants are expected to
become strongly negative near a first--order
transition \cite{BiLa}, even diverging to minus infinity in the
thermodynamic limit. We do not observe such a behavior.

Moreover, we studied the maximal specific heat $C_{max}$, close
to the AF--SF transition as a function of the lattice
size, $L$. At a first order transition in two
dimensions,  $C_{max}$ is predicted \cite{FiBe} to grow, for
sufficiently large lattices, like $C_{max} \propto L^2$. On
the other hand, at
an Ising--like transition in two dimensions, there is only
a logarithmic increase. Actually, we
fitted the simulation data for $C_{max}(L)$ to a
power law, $C_{max} \propto L^X$. We determined the
effective exponent $X_{eff}(L) = d \ln C_{max}/ d \ln L$. Simulations,
at $k_BT/J$ = 0.4 and 0.6, are performed at discrete
values of $L$, ranging from 20 to 200. Data for successive
sizes, $L_1$ and $L_2$, are fitted to
the discretized effective exponent $X_{eff}(L_e)$, with
$L_e = \sqrt{(L_1\times L_2)}$. We find rather small
effective exponents, being, at most, about 0.25, far from
the quadratic behavior characterising a transition of first--order.

In conclusion, for the anisotropic XY antiferromagnet on a square
lattice, we find no evidence for a direct transition of first
order between the AF and SF phases down to $k_B/J = 0.4$. Accordingly, there
is no indication for a bicritical point.

\subsection{Anisotropic XY antiferromagnet on a simple cubic lattice}

The crucial part of the phase diagram for the 3d
anisotropic XY model in a field is shown in Figure 4. It summarizes
present MC findings, setting the anisotropy parameter $\Delta$
equal to 0.8. At zero temperature, the AF configuration is stable up to
$H_{c1}/J = 3.6$ and the SF structure is stable up to
$H_{c2}/J = 10.8$.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov3.eps}
}
\caption{Phase diagram of the XY antiferromagnet
  with $\Delta=0.8$ on the simple cubic lattice.}
\label{fig:4}
\end{figure}

The transition lines between the disordered paramagnetic
phase and the two orderded phases AF and SF phases are
believed to belong
to the Ising universality class. In fact, our estimates
for the critical exponents  of the staggered magnetizations and
susceptibilities as well as of the specific heat agree with
the Ising values \cite{Vic}. Furthermore, for both
transitions, the critical Binder
cumulants, $U^*$, are found to be close to the
characteristic value for isotropic phase transitions
of 3d Ising type, $U^* \approx 0.465$ \cite{Hb}.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov4.eps}
}
\caption{Maximal specific heat, $C_{max}$, versus lattice size, $L$, near the
 AF--SF transition at $k_BT/J$ = 1.0 (circles), 1.3 (squares),
 1.45 (diamonds), 1.50 (triangles left), and 1.55 (triangles up). The
 straight line corresponds to $C_{max} \propto L^3$.}
\label{fig:5}
\end{figure}

The most interesting aspect of the phase diagram concerns
the existence, location, and
characteristics of the bicritical point, at which the AF, SF, and
paramagnetic phases may meet.

To identify a, possibly, first--order transition between the AF and
SF phases at low temperatures, we monitored the size--dependence
of the maximum in the specific heat, $C_{max}$. A first
transition of first order is expected \cite{FiBe} to be signalled by a
divergence of the form $C_{max} \propto L^3$, for sufficiently large
lattice sizes. In contrast, at a continuous transition in three
dimensions, one expects a power--law
behavior with an exponent $\alpha/\nu$, with
$\alpha$ and $\nu$ being the standard critical
exponents. For an Ising--type
transition in 3d \cite{Vic}, one obtains
$\alpha/\nu \approx 0.175$. As depicted in Figure 5, the
size--dependence typical for a first--order transition is observed to be
approximated closely at low temperatures. For instance, at
$k_BT/J$ = 1.0, the cubic size dependence, characteristic
for a first--order transition, is
approached already for $L$ = 16. Increasing the
temperature, the crossover towards the (almost) cubic behaviour is
shifted to larger lattices, being about $L$ = 30 at
$k_BT/J$ = 1.45. Obviously, at
higher temperatures, as shown in Figure 5 for $k_BT/J = 1.5$, it becomes
more and more demanding to reach the asymptotic regime. Presumably, the
difficulty is related to
the increase in the correlation length at the
first--order transition when getting closer to the
bicritical point. Actually, at the bicritical point in
the XY universality class, one expects no divergence in
$C$, but a cusp--like singularity, with the corresponding
standard critical exponent $\alpha$ being slightly
negative \cite{Vic}. Indeed, see Figure 5, at
$k_BT/J = 1.55$, $C_{max}$ increases only slightly
with the lattice size, indicating, probably, the closeness of the
bicritical point.

Another signature of the first--order transition at low temperatures
is provided by the longitudinal and transversal Binder
cumulants, $U^x$ and $U^y$. In fact, both quantities display, close
to the boundary between the AF and SF phases, negative
minima, choosing temperatures $k_BT/J$ ranging
from 1.0 to 1.50. Lattices with $L$ up to 36
were considered. The minima became more pronounced, when
lowering the temperature, reflecting the fact that then the transition
gets more strongly of first order, with a smaller
correlation length, in agreement with the behavior of the
specific heat, as discussed above. By increasing the lattice
size, the minima are getting deeper. In fact, in the
thermodynamic limit, the cumulant
is predicted to be minus infinity at the phase
transition of first order \cite{BiLa}.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov5.eps}
}
\caption{Binder cumulant $U^{xy}$ at constant temperature, $k_BT/J$ = 1.3,
  as a function of the field $H/J$ near the AF--SF transition for
  lattices with $L$ =16 (circles), 20 (squares) and 30 (diamonds).}
\label{fig:6}
\end{figure}

The data on the specific heat $C$ suggest that the transition between
the AF and SF phases is of first order up to temperatures of, at least,
about $k_BT/J \approx 1.5$.  To locate the bicrital point
and to characterize its nature, the Binder cumulant
$U^{xy}$ turns out to be quite useful as well. Its
critical value, at a
transition point in the XY--universality class in
three dimensions, has
been estimated \cite{HT} to be $U^{(xy)*} \approx 0.586$.

Of course, $U^{xy}$ is expected to approach, for large
lattices, in both ordered phases the limiting value
2/3. When the transition between the AF and SF
phases is of first order, as
suggested by the specific heat, the cumulant
displays, fixing the temperature
and varying the field (or vice versa), a minimum close to
the transition, $U^{xy}_{min}$, as exemplified in Figure 6. When
increasing the lattice size, the minimum is, eventually,
lowered, falling clearly below the value, 0.586, characterising
the XY universality class. Results of the present simulations
on the size--dependence of the height of the minimum for various
temperatures and fields are depicted in Figure 7.

\begin{figure}
\resizebox{0.95\columnwidth}{!}{%
  \includegraphics{figcov6.eps}
}
\caption{Minimum of the xy--cumulant $U^{xy}_{min}$ as a function
  of lattice size $L$ for various temperatures $k_BT/J$ = 1.3 (circles),
  1.45 (squares), 1.5 (diamonds), and 1.55 (triangles left) as well as for
  fixed field $H/J$ = 4.03 (triangles up). The critical cumulant in
  the XY universality class occurs at $U^{(xy)*} \approx 0.586$ \cite{HT}.}
\label{fig:7}
\end{figure}

Perhaps most interestingly, Figure 7 shows a rather drastic
change in the size dependence of the height of the minimum, when
getting closer to the bicritical point. Eventually, $U^{xy}_{min}$ starts
to increase for large lattices. It seems to approach, for
sufficiently large lattices, the
characteristic value of the 3d XY universality
class, $U^{(xy)*} \approx 0.586$, in the
temperature range between about $k_BT/J$ = 1.51 and 1.55. The
lower bound is inferred from the simulation data when fixing
the field at $H/J$ = 4.03, with
the corresponding phase transition occurring at
$k_BT/J \approx 1.518$. Accordingly, the present
simulations are consistent with the existence of a bicritical
point of XY type in the three--dimensional model. The
result confirms the previous theoretical prediction \cite{KNF,FN}
on the nature of the bicritical point, based
on renormalization group calculations. Note that at even
higher temperatures, either the two transition lines of the AF and SF
phases to the paramagnetic phase are well separated or
there is only the transition to the AF phase, as depicted
in Figure 3. Then, the transition is no longer signalled by a
minimum in $U^{xy}$.

Of course, to explore, in detail, finite--size
effects in the crossover region in the vicinity of
the bicritical point, data of very high accuracy for larger
lattices would be desirable. This feature, however, is beyond
the scope of our study.

\section{Summary}

We conclude that the topology of the phase diagram, in
the ($k_BT/J, H/J$) plane, of the anisotropic XY antiferromagnet
depends significantly on the lattice dimension. For the square
lattice, in agreement with a previous study, we find no evidence
for a bicritical point. Instead, in between the AF
and SF phases, there is a narrow intervening paramagnetic phase down
to the lowest temperatures we considered. In contrast, for the simple
cubic lattice, we find a line of direct transitions of
first order between the two ordered phases, leading, eventually, to
a bicritical point.

We locate the bicritical point
at $k_BT/J \approx 1.53 \pm 0.02$, based on finite--size analyzes
for the specific heat and Binder cumulants. The point belongs to
the XY universality class.

The qualitatively different phase diagrams in
two and three dimensions may be explained
by the bidirectional configurations, leading
to the highly degenerate ground state at the field separating the
AF and SF structures. In two dimensions, these configurations
seem to surpress the direct transition between the AF and SF phases, while
they are thermally less relevant in three dimensions. Such a
feature has been seen to hold in the related XXZ Heisenberg antiferromagnet
in a field as well. In that case, the existence of a bicritical
point (belonging to the Heisenberg universality class) in two dimensions
is already excluded by the Mermin--Wagner theorem.

We should like to thank especially Kurt Binder, David Landau,
Reinhard Folk, and Martin Holtschneider for useful information, remarks,
and discussions on the topic of this article. The article is dedicated
to Wolfhard Janke on the occasion of his sixtieth
birthday. We thank him for very helpful and pleasant
interactions.


\section{Introduction}

The respiration signal is of scientific relevance because it is controlled by the autonomic nervous system and the central nervous system. This means that respiration can be consciously controlled by the subject (unlike, for instance, heart rate). Different parts of the brain are responsible for autonomic and behavioural breathing: voluntary breathing is controlled by the cerebral cortex, while autonomous breathing is controlled by the brain stem. There are, however, projections from the cortex to the brain stem that allow the higher centres to influence metabolic breathing patterns. The autonomic respiratory pattern is therefore a complex interaction among the brain stem, the limbic system, and the cortex \citep{Homma2008}.

While respiration is an interesting signal from a theoretical perspective, there are some practical issues that limit its use in an HCI setting. Traditional respiratory measurement techniques, such as spirometers, which force the user to breathe through a tube, are often overly intrusive and cumbersome, while non-invasive methods can be too imprecise. The fact that respiration can be consciously controlled can act as a confounding factor. Also, respiration is highly susceptible to artefacts produced by speech and movement \citep{wientjes1992}.

\section{Background}
Although respiration (respiration rate) is an important topic in health science, in physiology studies in a games or learning context it is less well researched than are many other psychophysiological signals, suffering from neglect even in studies of the cardiovascular system, with which it interacts intimately \citep{Nesterov2005}. This is partly an interpretative issue, as respiratory control is both voluntary and involuntary \citep{Harver2000}, but equally it is a function of the uses served by respiration -- the heart and brain together account for less than 3\% of body weight but more than 30\% of oxygen usage. Accordingly, states of low mental workload and high metabolic rate may be reflected in the respiration rate similarly to states of high mental workload and low metabolic rate. However, it has been shown that increased respiration rate results when subjects are sitting and attentively wakeful (e.g., listening to a story), as opposed to sitting without paying attention or with closed eyes. Thus, in similar metabolic situations, increased respiration indicates an increase in attention \citep{Harver2000}.
In consequence, given an attentionally demanding situation with static metabolic demand, respiration should increase in line with the subject's engagement. In other words, the more attention paid (barring metabolic variation), the higher the respiration rate.

In its effects, respiration is closely interwoven with most other physiological signals of interest, because it provides the oxygen that underlies localised mental activity, the energising of the musculature, and heart rate response \citep{Harver2000}. Respiration affects heart rate variability through a process called `respiratory sinus arrhythmia' (RSA) whereby the inter-beat interval of the heart is shortened during inspiration and prolonged during expiration \citep{Yasuma2004}. It may therefore be necessary to control for breathing in experiments that include HRV.

In addition to the respiratory rate, there are several (more specific) indices that can be derived from the respiratory pattern, such as inspiratory time, inspiratory volume, and inspiratory pause. A study by \cite{Boiten1998} explored the effect of film scenes, reaction time, and cold pressor tasks on affect-related respiratory patterns. The study showed that there were clear effects on respiratory patterns whenever the films elicited amusement or disgust, in that amusement induced a decrease in inspiratory time and tidal volume, while disgust would elicit prolonged inspiratory pauses (breath-holding). The reaction time task induced a relatively fast, shallow, and regular breathing pattern, while the cold pressor brought on a substantial increase in respiratory volume, an increase in post-inspiratory pause duration, and a large amount of breath-to-breath variability in the pattern of breathing. Accordingly, in the design of HCI applications it can be useful to pay attention to more specific details of respiration than respiration rate alone (see \citet{Kreibig2010} for a comprehensive list of possible respiratory pattern indices).

Frequent sighing has been shown to be associated with anxiety. \cite{blechert2006identifying}, for instance, studied the physiological correlates of anxiety and found that sighing frequency increased 150$\%$ when the participants thought that they might receive an electric shock.

\section{Methods}
Several methods exist for measuring respiration. The flow of air in and out of the lungs can be measured directly by means of a spirometer. However, while very accurate, this apparatus forces the user to wear a nose clip and to breathe through a tube, which makes it unsuitable for most HCI purposes. Therefore, instead of measuring respiration directly, it is often more convenient to measure the movement of the chest and the abdomen, by using respiratory inductance plethysmography (RIP). In this method, an elastic band is wrapped around the abdomen and around the chest. It is also possible to use only a single belt, around the abdomen. However, to estimate tidal volume (i.e., the total volume of air inspired and expired with each breath), it is necessary to use two belts, enabling both the abdomen and the chest to be measured (see Figure~\ref{fig.resp_belts}). Nonetheless, because of the physiological uniqueness of every individual, the measurement of tidal volume with two belts requires calibration of the system for each user, through use of a spirometer \citep{wientjes1992}.

\begin{figure}[!ht]
	\centering
	\includegraphics[scale=0.5]{respiration_belts}
	\caption{Schematic diagram of the canonical locations for respiration belt sensors, on the chest and at the abdomen. Respiration belts are typically a flexible cable transducer (pictured here as a red sinusoid) attached to an elasticated cloth belt (here, a grey rectangle), terminating in connectors to the amplifier (represented by red and black squares).}
	\label{fig.resp_belts}
\end{figure}

\section{Applications}
In an HCI context, respiration has been used mostly as an explicit control signal. ~\citet{Moraveji2011} studied whether they could influence the respiration of users during an information-processing task.
Reducing breathing rate has been shown to reduce stress and anxiety, so the aim was to build a system that helps users to regulate their breathing to this end. It displayed a semi-translucent bar in the bottom half of the screen that served as guidance for the user. Users were instructed to inhale when the bar goes up and exhale as it moves down. The system was designed to illustrate a breath rate 20$\%$ below the user's baseline. It indeed reduced users' breathing rate -- the authors found that users were able to decrease their breathing rate significantly while still able to continue performing information work such as research, writing, and programming.

The use of respiration in a more entertainment-oriented context was studied by ~\citet{Kuikkaniemi2010}, who explored the use of biofeedback in computer gaming. They developed a first-person shooter game wherein the fire rate and recoil were linked to the player's respiration. In their experiment, they studied how biofeedback was perceived when the players were aware of the adaptation as compared to when they were unaware of the biofeedback. The researchers established that conscious control through respiration was more immersive and rewarding, thereby highlighting the nature of respiration as a signal that can be dealt with as either an implicit or an explicit form of interaction.

\section{Conclusions}
Respiration is a signal that is rarely used on its own. More often, is serves as an auxiliary signal either to complement information from other signals, such as heart rate in the case of respiratory sinus arrhythmia, or to serve as part of a machine learning approach that draws together information from several physiological signals. When used alone, respiration is usually employed as an explicit control mechanism, offering an additional input channel that can complement manual input.

\bibliographystyle{plainnat}


\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.