text:
\section*{References}}

\begin{document}

\title{Autocorrelation Function Characterization of Continuous Time Markov Chains}

\author[mec]{G.~Rama Murthy }
\ead{rama.murthy@mechyd.ac.in}

\author[mcm]{D.G.~Down \corref{cor1}}
\ead{downd@mcmaster.ca}

\cortext[cor1]{Corresponding author}

\address[mec]{Department of Computer Science, Mahindra Ecole Centrale (MEC), Hyderabad, India}

\address[mcm]{Department of Computing and Software, McMaster University, Hamilton, Canada}

\begin{abstract} We study certain properties of the function space of
autocorrelation functions of Unit Continuous Time Markov Chains (CTMCs). It is shown
that under particular conditions, the $L^p$ norm of the autocorrelation function of
arbitrary finite state space CTMCs is infinite. Several interesting inferences are made for
point processes associated with CTMCs/ Discrete Time Markov Chains (DTMCs).
\end{abstract}

\maketitle

\section{Introduction}

Many natural/artificial phenomena are endowed with non-deterministic dynamic behavior. Stochastic/random processes with finite/countable/uncountable state space are utilized to model such dynamic phenomena. Random processes which are stationary and whose state space consists of two ($\{ +1, -1 \}$) values arise naturally in many applications e.g.\ bit transmission in communications systems. The characterization of  the autocorrelation function of such unit stochastic processes is considered an important problem \cite{Mcm}. Several interesting properties of such autocorrelation functions are studied in \cite{She,Mas}.

Wide sense stationary (or even strictly stationary) random processes naturally arise as stochastic models in a variety of applications. They also arise in time series models (AR, ARMA processes) of natural/artificial phenomena. In most such models, the autocorrelation function, $R(\tau)$ is integrable and hence the power spectral density (the Fourier transform of $R(\tau)$) exists. There are a number of applications (such as random telegraph signals) where the underlying stochastic process can be in only one of two states, so-called unit random processes. An example can be taken from detection theory. Suppose that a random process, $x(t)$ is provided as input to a threshold detector, where the output $y(t)$ is the sign of $x(t)$, i.e., $y(t)=+1$ if $x(t) \geq 0$ and $y(t)=-1$ if $x(t)<0$. Thus the output is a unit process and can be shown to be Markovian under some conditions on $x(t)$. In this case, the autocorrelation of $y(t)$ is in general not integrable.  More generally, quantization of a general random process leads to a finite-state random process, which is often Markovian, so the study of more general finite-state processes is also of interest. With this in mind, we are motivated to study the function space of finite state Markovian processes. Masry \cite{Mas} has studied the functional space properties of stationary unit random processes. However, the study of the integrability of $R(\tau)$ was not undertaken.  To the best of our knowledge, integrability, or more generally the $L^p$-norm of $R(\tau)$ for finite state Continuous Time Markov Chains (CTMCs) has not ben investigated. In this paper, we determine conditions under which the autocorrelation function is not integrable, and by extension conditions under which the $L^p$-norm of $R(\tau)$ approaches zero as $p\rightarrow \infty$.

To put our work into context, there is related work in three directions: the characterization of autocorrelation functions of random processes, the characterization of point processes, and the use of autocorrelation properties in the analysis of stochastic models, in particular the analysis of queues. For the characterization of autocorrelation functions, we point the reader to work in time series analysis \cite{deb11,deglac03,ino08} and in telecommunications \cite{eunjinhonson04}. Point process characterization has been studied in \cite{hagvanmol99,hua90,ivamerpla07}. Properties of autocorrelation functions have been employed to determine appropriate simulation strategies for queues \cite{whi91} and is a feature of modelling arrival traffic to queues, using Markovian Arrival Processes, see \cite{klelinloh03}, for example.

This paper is organized as follows. In Section~\ref{s2}, the autocorrelation function of a Unit ($\{ +1, -1 \}$ state) CTMC is computed and the structure of the function space is studied. In Section~\ref{s3}, the autocorrelation function of a symmetric state space CTMC (i.e., the state space is $\{ -N, -N+1,\ldots,-1, +1, \ldots+N-1, +N \}$) is computed and the finiteness of its $L^{p}$ norm is discussed. It is shown that under some conditions, the autocorrelation function is not integrable. In Section~\ref{s4}, various interesting inferences are made for point processes. Finally, the paper concludes in Section~\ref{s5}.

  \section{Auto-Correlation Function of Homogeneous Unit CTMC: Integrability}
  \label{s2}

In the following, we consider a homogeneous Continuous Time Markov Chain (CTMC) (with the state space being $\{ +1, -1 \}$) whose generator matrix is denoted by $Q$. We assume that the resulting stochastic process is wide sense stationary (and not necessarily strictly stationary).

The autocorrelation function is given by
\[ R(t,t+\tau) = E [X(t)X(t +\tau)]. \]
Since the CTMC is homogeneous, we have that
\[ R( t,t+\tau) = R(\tau) = E[X(0)X(\tau)]. \]
Since the CTMC is a unit random process (i.e., the state space is $\{ +1,-1 \}$), we have that
\[ R(\tau) = P \{ X(0) = X(\tau) \} - P \{ X(0)  \neq X(\tau) \}. \]
Combining this with
\[ P \{ X(0)  \neq X(\tau) \} = 1 - P\{ X(0) = X(\tau) \}  \]
yields
\begin{equation}
\label{e1}
R(\tau) = 2P\{ X(0) = X(\tau) \}  -1.
\end{equation}
It remains to compute $P \{ X(0)=X(\tau) \}$.
\begin{eqnarray}
\lefteqn{P\{ X(\tau) = X(0) \}} \nonumber\\
& = & P\{ X(\tau) = +1,X(0) = +1\} + P\{ X(\tau) = -1,X(0) = -1 \} \nonumber \\
& = & P\{ X(\tau) = +1 | X(0) = +1 \} P\{ X(0) = +1 \} + P\{ X(\tau) = -1 | X(0) = -1\} P\{ X(0) = -1 \}. \label{e2}
\end{eqnarray}
The conditional probabilities in the above expression are computed using
the transient probability distribution for a homogeneous CTMC. This computation requires the determination of the matrix exponential of $Q$, i.e., $e^{Qt}$. Using the
spectral representation, $e^{Qt}$ is computed below.

\begin{itemize}
\item Computation of $e^{Qt}$. Let the generator matrix of $\{ +1, -1 \}$ state CTMC (i.e., Unit CTMC) be given by
\[ Q = \begin{bmatrix}
   - \alpha & \alpha \\
  \beta & - \beta \\
  \end{bmatrix}.
\]
Since $Q$ is a rank one matrix, the eigenvalues are $\{ -( \alpha + \beta),0 \}$. Also, let the right eigenvectors be given by $\{ \bar{g}_1, \bar{g}_2 \}$ (i.e., column vectors) and let the left eigenvectors be given by $\{ \bar{f}_1,\bar{f}_2 \}$ (i.e., row vectors). Thus, we have that
\[ e^{Qt} = e^{-(\alpha + \beta)t} \bar{g}_1 \bar{f}_1 + \bar{g}_2 \bar{f}_2. \]
\end{itemize}

Now, we discuss the computation of eigenvectors.

\begin{itemize}
\item Computation of $\bar{g}_1$. By definition of right eigenvector, we have that
\[ Q\bar{g}_1 = \  - (\alpha + \beta) \bar{g}_1. \]
By letting $\bar{g}_1 = [ x \ y ]^T$ and using the expression for the generator matrix $Q$,
\[ \begin{bmatrix}
 - \alpha & \alpha \\
\beta & - \beta \\
\end{bmatrix} \begin{bmatrix}
x \\
y \\
\end{bmatrix} = - (\alpha + \beta)\ \begin{bmatrix}
x \\
y \\
\end{bmatrix}. \]
Solving the above linear equations, we have $y = - \frac{\beta}{\alpha} x$. Thus letting $x = 1$, we have $y = - \beta / \alpha$. Hence
\[ \bar{g}_1 = \begin{bmatrix}
1 \\
\frac{-\beta}{\alpha} \\
\end{bmatrix}. \]

\item Computation of $\bar{g}_2$. Since $Q$ is the generator of the CTMC, we have
\[ \begin{bmatrix}
   - \alpha & \alpha \\
  \beta & - \beta \\
  \end{bmatrix} \begin{bmatrix}
  1 \\
  1 \\
  \end{bmatrix}  = \begin{bmatrix}
  0 \\
  0 \\
  \end{bmatrix}. \]
Thus, the right eigenvector corresponding to the zero eigenvalue is given by
\[ \bar{g}_2 = \begin{bmatrix}
1 \\
1 \\
\end{bmatrix}. \]

\item Computation of $\bar{f}_1$. Since the generator matrix $Q$ is diagonalizable, we have
\[ \bar{f}_1 \bar{g}_1 = 1 \ \mbox{and} \ \bar{f}_1 \bar{g}_2 = 0. \]
Letting $\bar{f}_1 = [a \ b]^T$, $\bar{f}_1 \bar{g}_2  = 0$ implies $a = -b$.  Furthermore, $\bar{f}_1 \bar{g}_1 = 1$ implies that $b = -\frac{\alpha}{\alpha + \beta}$. Hence
\[
\bar{f}_1 = \left[ \begin{array}{cc}
\frac{\alpha}{\alpha + \beta} & - \frac{\alpha}{\alpha + \beta}
\end{array} \right].
\]

\item Computation of $\bar{f}_2$. Letting $\bar{f}_2 = [ c \ d]^T$, $\bar{f}_2 \bar{g}_1 = 0$ implies that $c = d \frac{\beta}{\alpha}$. Furthermore, $\bar{f}_2 \bar{g}_2 = 1$ implies $d = \frac{\alpha}{\alpha + \beta}$. Hence, we have that
\[ \bar{f}_2 = \left[ \begin{array}{cc}
\frac{\beta}{\alpha + \beta} & \frac{\alpha}{\alpha + \beta}
\end{array}  \right]. \]
\end{itemize}
Thus, we have effectively computed the matrix exponential $e^{Qt}$.

The transient probability distribution of our homogeneous, Unit CTMC is given by
$\bar{\pi}(\tau) = [P\{ X(\tau) = +1\} \ P\{ X(\tau) = -1 \}]$.
Using the matrix exponential, the transient probability distribution is
given by
\begin{eqnarray*}
\bar{\pi}(\tau) & = & \bar{\pi}(0)\ e^{Q\tau}\\
& = & \bar{\pi}(0)\ \left[ e^{-( \alpha + \beta)\tau} \bar{g}_1 \bar{f}_1 + \bar{g}_2 \bar{f}_2 \right]\\
& = & \bar{\pi}(0)\ \left[ \begin{matrix}
P\{ X(\tau) = +1 | X(0) = +1\} & P\{ X(\tau) = -1 | X(0) = +1\}  \\
P\{ X(\tau) = +1 | X(0) = -1\} & P\{ X(\tau) = -1 |X(0) = -1\} \\
\end{matrix} \right].
\end{eqnarray*}
We now provide an explicit expression for the matrix exponential $e^{Q \tau}$.
\begin{eqnarray*}
e^{Q\tau} & = & e^{-( \alpha + \beta)\tau} \begin{bmatrix}
1 \\
\frac{- \beta}{\alpha} \\
\end{bmatrix} \left[ \begin{matrix}
\frac{\alpha}{\alpha + \beta} & -\frac{\alpha}{\alpha + \beta} \\
\end{matrix}  \right] + \begin{bmatrix}
1 \\
1 \\
\end{bmatrix}\left[ \begin{matrix}
\frac{\beta}{\alpha + \beta} & \frac{\alpha}{\alpha + \beta} \\
\end{matrix}  \right]\\
& = & \begin{bmatrix}
e^{-( \alpha + \beta)\tau}\frac{\alpha}{\alpha + \beta} + \frac{\beta}{\alpha + \beta} & -e^{-( \alpha + \beta)\tau}\frac{\alpha}{\alpha + \beta}  + \frac{\alpha}{\alpha + \beta}  \\
 -e^{-( \alpha + \beta)\tau}\frac{\beta}{\alpha + \beta} + \frac{\beta}{\alpha + \beta} & e^{-(\alpha + \beta)\tau}\frac{\beta}{\alpha + \beta} + \frac{\alpha}{\alpha + \beta} \\
\end{bmatrix}.
\end{eqnarray*}
We utilize the expression for $e^{Q\tau}$ in (\ref{e2}) to compute $P\{ X(\tau) = X(0) \}$.
\begin{eqnarray}
\lefteqn{P\{ X(\tau) = X(0) \}} \nonumber\\
& = & \ P\{ X(\tau) = +1 | X(0) = +1\}  P\{ X(0) = +1 \} + P\{ X(\tau) = -1 | X(0) = - 1\}  P\{ X(0) = -1 \} \nonumber\\
& = &
\left[ e^{- (\alpha + \beta)\tau}\frac{\alpha}{\alpha + \beta} \right] P\{ X(0) = +1 \}  + \left[ e^{-(\alpha + \beta)\tau}\frac{\beta}{\alpha + \beta} \right] P\{ X(0) = -1  \} \nonumber \\
& & + \frac{\beta}{\alpha + \beta}P\{ X(0) = +1\  \}  + \frac{\alpha}{\alpha + \beta} P\{ X(0) = -1 \} \label{e3}
\end{eqnarray}

Now, we compute the equilibrium probability distribution of the unit CTMC. Let the equilibrium probability distribution be denoted by the row vector
$\bar{\pi}$. Thus
$\bar{\pi} = [\pi_{1} \ \pi_{- 1}]$ and $\overline{\pi}Q = 0$. Solving
the linear equations, we have $\pi_{1} = \frac{\beta}{\alpha + \beta}$ and $\pi_{- 1} =
\frac{\alpha}{\alpha + \beta}$.

\textbf{Note:} We now assume that the initial probability distribution
equals the equilibrium probability distribution, i.e., $\overline{\pi}(0) = \overline{\pi} = [\pi_{1} \ \pi_{- 1}]$.
In this case, the transient probability distribution is also equal to the equilibrium probability distribution, which allows us to write
\[ P \{ X (\tau) = X(0) \} = \frac{\alpha^{2} + \beta^{2}}{(\alpha + \beta)^{2}}. \]
Now substituting for $P\{ X(\tau) = X(0) \}$
in (\ref{e1}), we have that
\[ R(\tau) = 2P\{ X(\tau) = X(0)  \}  - 1  = 2\left[ \frac{\alpha^{2} + \beta^{2}}{(\alpha + \beta)^{2}} \right] - 1. \]
Simplifying, we have that
\begin{equation} \label{e:R}
R(\tau) = \left[ \frac{(\alpha - \beta)^{2}}{(\alpha + \beta)^{2}} \right].
\end{equation}
Thus, the autocorrelation function is a constant when the initial
probability distribution equals the equilibrium probability
distribution. Further, if $\alpha = \beta$, we have that $R(\tau)=0$, i.e., it is identically zero. Since the initial probability distribution is assumed to equal the
equilibrium probability distribution, with $\alpha = \beta$, we have
\[ E[X(0)] = (1)\frac{\beta}{\alpha + \beta} +(-1) \frac{\alpha}{\alpha + \beta} = 0, \]
i.e., with mean zero initial random variable, (\ref{e:R}) holds for
the autocorrelation function.

Suppose, we consider the case where the initial probability distribution
is not equal to the equilibrium probability distribution and is
arbitrary, i.e.,
\[ P\{ X(0) = + 1\} = q \quad \mbox{and} \quad P\{ X(0) = -1\} = 1 - q. \]
Substituting in (\ref{e3}), we have that
\begin{eqnarray*}
\lefteqn{P\{ X(\tau) = X(0)\} }\\
& = & \left[ e^{-(\alpha + \beta)\tau}\frac{\alpha}{\alpha + \beta} \right] q  + \left[ e^{-(\alpha + \beta)\tau}\frac{\beta}{\alpha + \beta} \right] (1-q) + \frac{\beta}{\alpha + \beta} q  + \frac{\alpha}{\alpha + \beta} (1 - q)\\
& = &
e^{-(\alpha + \beta)\tau}\ \left( \frac{\alpha - \beta}{\alpha + \beta} \right)q + e^{-(\alpha + \beta)\tau}\left( \frac{\beta}{\alpha + \beta} \right) + \left( \frac{\beta - \alpha}{\alpha + \beta} \right)q + \frac{\alpha}{\alpha + \beta}.
\end{eqnarray*}
Substituting in $R(\tau) = 2P\{ X(\tau) = X(0) \}  - 1$,
\[ R(\tau) = 2\left[ e^{-(\alpha + \beta)\tau} \left( \frac{\alpha - \beta}{\alpha + \beta} \right)q + e^{-(\alpha + \beta)\tau}\left( \frac{\beta}{\alpha + \beta} \right) \right] + 2 \left[ \left( \frac{\beta - \alpha}{\alpha + \beta} \right)q + \frac{\alpha}{\alpha + \beta} \right] - 1. \]
Thus
$R(\tau)$ is a sum of a function of $\tau$ and a constant term, which we write as $R(\tau)=f(\tau)+c$, where
\[ c = 2\left[ \left( \frac{\beta - \alpha}{\alpha + \beta} \right)q \right] + \frac{\alpha - \beta}{\alpha + \beta}. \]
If $q = \frac{1}{2}$, we have that $c$ is zero and thus
\[
R(\tau) = e^{-(\alpha + \beta)\tau}\ \left[ \left( \frac{\alpha - \beta}{\alpha + \beta} \right) + 2\left( \frac{\beta}{\alpha + \beta} \right) \right]
= e^{-(\alpha + \beta)\tau}.
\]
Now, if $\alpha = \beta$, and since the
autocorrelation function is symmetric, we have
\[ R(\tau) = e^{-2\alpha |\tau|}. \]

As discussed earlier, an interesting problem is to
characterize the autocorrelation function of Unit CTMCs, i.e., $\{ +1, -1 \}$-state CTMCs.
The following lemma presents one possible characterization of the function
space of autocorrelation functions of Unit CTMCs

\begin{lemma}
Consider a Unit CTMC with $\alpha \neq \beta$. The
autocorrelation function of such a Unit CTMC is not in
$L^{p}[R(\tau)]$ for any $p \geq 1$ (the $L^p$ norm of the autocorrelation function is infinite) . However, as $p$ tends to $\infty$, the $L^{p}$ - norm of the autocorrelation function, $R(\tau)$, approaches a finite constant. Further the $L^{\infty}$-norm is equal to one.
\end{lemma}

{\bf Proof.} In the case of unit CTMCs with $\alpha \neq \beta$ and the initial probability distribution is equal to the equilibrium probability distribution,
the expression for autocorrelation function is identically (for all
time) a non-zero constant which is less than one in  magnitude. Thus, such a function is not in $L^{p}(R)$ for any
$p \geq 1$. Thus, such an autocorrelation function is necessarily not integrable.
Since the constant term is less than one in magnitude, it readily
follows that the $p^{th}$ power of it approaches zero. Hence the
$L^{p}$-norm of the autocorrelation function, $R(\tau)$, approaches a finite constant.

In the case of arbitrary initial probability distribution, if
$q \neq \frac{1}{2}$, then $c$ is non-zero. Further if
$q \leq \frac{1}{2}$, $0 \leq c \leq 1$ and
$q \geq \frac{1}{2}$, $-1 \leq c \leq 0$. Thus, even in
this case, the autocorrelation function of such a Unit CTMC is not in
$L^{p}(R)$ for any $p \geq 1$. Hence,
it is necessarily not integrable. Also, as the constant term is strictly
less than one in magnitude, its $p^{th}$ power approaches zero, the $L^{p}$-norm of the autocorrelation function, $R(\tau)$, approaches a finite constant. \hfill $\diamondsuit$

In the following discussion, we generalize the above results. We
consider symmetric state space CTMCs (with state space $\{ -N,
-N+1,\ldots,-1, +1, \ldots, +N-1, +N \}$) as well as arbitrary state space CTMCs. It is shown that the existence
of an equilibrium probability distribution ensures that the expression for
autocorrelation function has a constant part that is not
zero under some conditions.

\section{Auto-Correlation Function of Homogeneous Finite
  State Space CTMC}
\label{s3}

\subsection{Asymmetric State Space}

We now prove that for any finite state space CTMC, the autocorrelation function is not integrable and in
fact the $L^{p}$-norm of the autocorrelation function, $R(\tau)$,
is infinite for any $p \geq 1$. Without loss of generality, the state space of the CTMC is assumed to be $\{ +1, +2, \ldots, +N \}$.

For a homogeneous CTMC, we have that $R(t,t + \tau) = R(\tau) = E[X(0)X(\tau)]$. Since the CTMC under consideration has finite state space, we have that
\begin{eqnarray*}
R(\tau) & = & \sum_{i = 1}^{N}\sum_{j = 1}^{N} ij P\{ X(0) = i,X(\tau)=j  \}\\
& = & \sum_{i = 1}^{N}\sum_{j = 1}^{N} ij P\{ X(\tau) = j | X(0) = i \} P\{ X(0) = i \}\\
& = & \sum_{i = 1}^{N}\sum_{j = 1}^{N} ij \left( e^{Q\tau} \right)_{ij} P\{ X(0) = i  \}.
\end{eqnarray*}
But, we have that
\[ e^{Q\tau} = \sum_{k = 1}^{N}e^{\gamma_{k}\tau} E_{k}, \]
where $E_{k}$ is the residue matrix such that $E_{k} = \bar{f}_k\ \bar{g}_k$ with $\bar{f}_k$ being the right eigenvector of $Q$ corresponding to the eigenvalue $\gamma_k$ and $\bar{g}_k$ being the left eigenvector of $Q$ corresponding to the eigenvalue $\gamma_k$. Let
$|\gamma_1| \geq |\gamma_{2}| \geq \cdots \geq | \gamma_{N} | = 0$ (since $Q$ is the generator matrix, we have that $\gamma_{N} = 0$).

\begin{eqnarray}
R(\tau) & = & \sum_{i = 1}^{N} \sum_{j = 1}^{N} ij \left( e^{Q\tau} \right)_{ij} P\{ X(0) = i \} \nonumber \\
& = & \sum_{i = 1}^{N} \sum_{j = 1}^{N} ij \left[ \sum_{k = 1}^{N}e^{\gamma_{k}\tau} E_{k}\ \right]_{ij} P\{ X(0) = i  \} \nonumber \\
& = & \sum_{i = 1}^{N} \sum_{j = 1}^{N} ij P\{ X(0) = i \} \left[ \sum_{k = 1}^{N - 1}e^{\gamma_{k}\tau} E_{k} \right]_{ij}
+ \sum_{i = 1}^{N}\sum_{j = 1}^{N}
ij P\{ X(0) = i  \}[ \bar{f}_N \ \bar{g}_N]_{ij} \label{e:star}
\end{eqnarray}
We would like to write (\ref{e:star}) as $f(\tau)+c$, where
\[
c = \sum_{i = 1}^{N}\sum_{j = 1}^{N} ij P\{ X(0) = i \} [\bar{f}_N \ \bar{g}_N]_{ij}.
\]
We must justify that $c$ is not a function of $\tau$. To do this, it
is first noted that $\bar{f}_N$ is the equilibrium
probability vector of the CTMC, i.e., $[ \pi_{1} \ \pi_{2} \cdots \pi_{N}]$.
It is also the left eigenvector of the generator matrix $Q$ corresponding to the
eigenvalue zero. Also
$\bar{g}_N = [ 1 \ 1\cdots 1]^{T}$, i.e.,
a column vector of all ones. It is the right
eigenvector of the generator matrix $Q$ corresponding to the eigenvalue zero. By the Perron-Frobenius Theorem, all of the components of the vectors $\{
\bar{f}_N, \bar{g}_N \}$ are non-negative. Thus,
$c$ involves a sum of non-negative quantities and does not depend on
$\tau$. It readily follows that $[ \bar{f}_N \ \bar{g}_N]_{ij} = \pi_{j}$. Let $P\{ X(0) = i \} = q_{i}$. Hence, we have that
\begin{eqnarray*}
c & = &  \sum_{i = 1}^{N}\sum_{j = 1}^{N} ij q_{i} [\bar{f}_N \ \bar{g}_N]_{ij} \\
& = & \left( \sum_{i = 1}^{N} iq_{i} \right) \left( \sum_{j = 1}^{N} j\pi_{j} \right) \\
& = & E[X(0)]E[Z],
\end{eqnarray*}
where $Z$ is the random variable associated with the equilibrium
probability distribution. Thus, $c$ is zero if and only if either $E[X(0)]$ or $E[Z]$ is zero (or both). For our choice of state space, this is not possible.

\textbf{Note:} If the initial probability distribution equals the
equilibrium probability distribution, then $c = E[X(0)]^2 = E[Z]^2$.

\textbf{Note:} In general, the state space of a finite state space CTMC
could assume negative values (as in the case of a Unit CTMC). It can be
readily verified that if the state space of CTMC is the set $\{ -N, -N+1,
\ldots, -1, +1, +2,\ldots, +N \}$ and the equilibrium probability
distribution is the uniform distribution, then $c$ equals zero.
Thus, such a condition is sufficient for $c$ to be zero. We
conjecture that the condition is also necessary for $c$ to be zero.

Now, $f(\tau)$ is a sum of decaying exponentials. It can be easily
verified that $f(\tau)$ is integrable. More
generally, $f(\tau)$ corresponds to a function which is in
$L^{p}(R)$ for $p \geq 1$. But, when $c$ is
non-zero, then
$R(\tau)$ is not in $L^{p}(R)$
for any $p \geq 1$.
If $c \neq 0$, then $\int (R(\tau))^{p} d\tau$ is infinite for every $p \geq 1$. Further if $|c| < 1$, then $\int (R(\tau))^{p}d\tau$ approaches zero as $p \rightarrow \infty$. The following lemma readily follows from the above discussion.

\begin{lemma}
Consider an arbitrary finite state CTMC with $c$
being non-zero. The autocorrelation function of such a unit CTMC is
not in $L^{p}(R)$ (where $R$ is the
real valued lag) for any $p \geq 1$,  i.e., such an autocorrelation
function is necessarily not integrable. However, as
$p$ tends to $\infty$, the $L^{p}$-norm of the autocorrelation function,
$R(\tau)$ approaches a finite constant.
\end{lemma}

{\bf Note:} We can provide further details about $c$, in the case that the
random process is a Unit CTMC.

Case (A): $X(0) = Z$, i.e., the initial probability distribution equals the
equilibrium probability distribution.
\[
c = E[X(0)]^2 = E[Z]^2 = (1)\frac{\beta}{\alpha + \beta}\  + (-1)\frac{\alpha}{\alpha + \beta} = \left( \frac{\alpha - \beta}{\alpha + \beta} \right)^{2}.
\]

CASE (B) : $X(0) \neq Z$, i.e., the equilibrium probability distribution is arbitrary, i.e., $P(X(0) = +1) = q$ and $P(X(0) = -1) = 1 - q$. Hence, we have that
\begin{eqnarray*}
E[X(0)] & = & q(1) + (1 - q)(-1) = 2q-1\\
E[Z] & = & \frac{\beta - \alpha}{\alpha + \beta}.
\end{eqnarray*}
Hence $E[X(0)]E[Z] = 2\left[ \left( \frac{\beta - \alpha}{\alpha + \beta} \right)q \right] + \frac{\alpha - \beta}{\alpha + \beta}$
which is consistent with our earlier result.

\textbf{Claim:} We now interpret the above results for arbitrary
Continuous Time Markov Chains (CTMCs) for which an equilibrium
distribution exists. Since $R(\tau) = E[X(0)X(\tau)]$,
we have that the following asymptotic result holds true.
\[ \lim_{\tau \rightarrow \infty} R(\tau) = E[X(0)]E[Z], \]
where $Z$ is the equilibrium random variable.

This result agrees with the fact that asymptotically the initial random
variable $X(0)$ and the equilibrium random variable are independent. Also, these results can be easily generalized to countable state space CTMCs.

\section{Finite State Space Continuous Time Markov Chains: Point
  Processes}
  \label{s4}

It is well known that the interarrival times of a Poisson process are exponentially distributed random variables. Also, the sojourn times in every state of a finite state CTMC are exponentially distributed random variables. This observation has been explored in \cite{cin75}, for example, to establish that when successive visits to a state of a CTMC are stitched together, a Poisson process naturally results. Hence, an arbitrary finite state CTMC can be viewed as a superposition of point processes. From a practical viewpoint, the superposition of point processes naturally arises in applications, such as packet streams in packet multiplexers. Such packet streams have been modelled in \cite{sriwhi86}, for example. Several versatile point processes have also been studied in \cite{hefluc86,neu79}, amongst others. Such Markovian point processes are actively utilized in queueing theoretic applications. In discrete time, sojourn times in any state of a DTMC are geometrically distributed random variables. In this sense, a DTMC can be viewed as a superposition of discrete time point processes (which could be dependent). Thus, homogeneous CTMCs/DTMCs provide interesting models  of superpositions of point processes. In other words, many well known results related to equilibrium/transient probability distributions of such Markov chains, effectively provide novel inferences related to the superposition of packet streams in communication networks.
Our main goal in this section is to  to derive several such inferences, considering CTMCs to illustrate the idea.

\subsection{Equilibrium/Transient Inferences related to Point Processes}

  Let $\overline{\pi} =
  [\pi_{1} \ \pi_{2} \cdots \pi_{N} ]$ be the
  vector of equilibrium probabilities.

\begin{enumerate}
\item
  Probability that in equilibrium the arrival is from
  the $j^{th}$ Point Process = $\pi_{j}$, for $1 \leq \ j \leq N$.
\item Probability that at time $\tau$ (under transient conditions), the
  arrival is from the $j^{th}$ point process, $P\{ X(\tau) = j  \} = [\overline{\pi}(0)e^{Q\tau}]_{j}$,
 i.e., $j^{th}$ component of such a vector for $1 \leq \ j \leq N$.
\item
  Conditional Probability that at time $\tau$ (under transient
  conditions), the arrival is from the
  $j^{th}$ point process given that at time zero  the  process was  in state $i$ is equal to
$[e^{Q\tau} ]_{ij}$ for $1 \leq  i \leq N$ and $1 \leq  j \leq N$.

  Many such inferences can be made related to the point processes in
  superposition given the interesting quantities associated with the
  transient/equilibrium probability distribution of corresponding DTMC/CTMC models.
\end{enumerate}

\section{Conclusion}
\label{s5}

We have computed the autocorrelation function of a unit CTMC
and the conditions for integrability (more generally finiteness of the
$L^{p}$-norm) were established. More generally, the function space structure of
  arbitrary finite state space CTMC was explored. Interesting inferences
  related to point processes (in a superposition point process) were
  made based on their relationship to finite state space Markov chains.

\vspace{.1in}

\noindent {\bf Acknowledgment.} The first author is supported by internal funding for research and development from Mahindra \'{E}cole Centrale, Hyderabad (for the year 2018-19). The second author is supported by NSERC under the Discovery Grant program.


text:
\section{INTRODUCTION}

The majority of stars form in clusters of 100 or more members \citep{lada03}. The dynamics and distribution of stars within clusters have important ramifications for
planet formation, stellar multiplicity, mass segregation, cluster disruption, and any other property of stars and planets that can be influenced by interactions between stars in the dense environment of a young cluster.
For example, \citet{adams01} show that close encounters with neighboring stars can disrupt the outer parts of planetary systems (however, see \citealt{megeath08}). At an earlier stage, protoplanetary disks in dense clusters are subject to dynamical disruption (e.g.\ \citealt{spurzem09}), photoevaporation (e.g.\ \citealt{adams04}), and contamination with short-lived radio-nuclides \citep[e.g.][]{gounelle09, takigawa08, ouellette07}. On larger scales, radiation from neighboring stars in a cluster is an important feedback effect that can modify the stellar IMF \citep{adams08, offner09}, while gravitational interactions among stars lead to mass segregation \citep{bonnell98, mcmillan07, allison09}. All of these phenomena depend critically on the kinematics of newly formed stars, since this determines the frequency of star-star encounters. Kinematics also influence the fraction of star clusters that remain bound rather than dissolving into the field \citep[e.g.][]{baumgardt07}.

Observations of the kinematics of young stars are challenging because the stars spend much of their youth enshrouded by gas, making spectroscopy difficult. However, measurements have been possible in a few regions, most notably the Orion Nebula Cluster (ONC), a particularly interesting and well-studied young star cluster containing thousands of optically revealed stars as well as traces of remnant molecular gas. The ONC is observed to be centrally concentrated and mass segregated, with the most massive stars residing in the center (e.g.~\citealt{hillenbrand98}). Radial velocity measurements of stars within the ONC indicate a substantial stellar velocity gradient and significant correlation between the local gas and star velocities \citep{furesz08, tobin09}. As a result, \citet{tobin09} suggest that the ONC is kinematically young and not dynamically relaxed, which they interpret as requiring that the parent cloud be in a state of dynamical collapse \citep{hartmann01, bonnell03}, rather than evolving in quasi-equilibrum \citep{tan06a, kandt07, huff07}. \citet{proszkow09} use N-body simulations to demonstrate that collapsing clusters with subvirial initial conditions can indeed generate such observational signatures.

However, it is far from clear how to extrapolate from observations of stellar kinematics to the conditions within the host giant molecular cloud (GMC), particularly since not all components of a GMC share the same velocity dispersion. Instead, the dispersion of the line-of-sight (centroid) gas velocity decreases with density (e.g.~\citealt{blitz86, heyer01}). Gas velocity dispersions in prestellar cores are comparable to the sound speed and are much less than the global cloud virial velocity \citep{goodman93, andre07, kirk07, roso08}. Furthermore, the dispersion of the bulk velocities of the cores themselves is found to be subvirial \citep{andre07}, indicating the the relative motion between cores is
small. Molecular gas in star forming regions such as Perseus and Ophiuchus also exhibits large-scale linear velocity gradients similar to those observed in the Orion stellar population \citep{andre07, kirk07, roso08}.  Perseus, in particular, resembles an elongated filament with a gradient of $\sim0.3$ km s$^{-1}$ pc$^{-1}$ along its length, similar to the ONC.

These observations of dense gas suggest that one cannot simply assume that, because stars are subvirial, the parent cloud must be too. Indeed, for gas the opposite appears to be the case. All of these subvirial signatures for dense components coexist with large-scale velocity dispersions that are consistent with virial balance \citep[e.g.][]{blitz07a, heyer08a, bolatto08a}. Simulations of turbulent, virialized clouds are able to reproduce the observed kinematics of dense gas \citep{padoan01, walsh04, Offner08b}, while simulations of clouds undergoing subvirial collapse have trouble producing quiescent cores like those we observe \citep{Offner08b}. Nor do clouds appear to convert large fractions of their mass into stars in a single dynamical time, as expected for a subvirial collapse  \citep[e.g.][]{klessen00, bonnell03}. Instead, observations both of individual star-forming regions \citep{evans09} and entire galaxies \citep{kandt07} show that, even for gas as dense as the cloud that formed the ONC, $\la 5\%$ of the mass per dynamical time can form stars.

The complicated relationship between the kinematics of dense and diffuse gas in molecular clouds suggests that it would be profitable to perform a systematic investigation of how the kinematics of newborn stars relate to the properties of the clouds from which they formed. Such an investigation is also of broader interest, since, as noted above, the kinematics of young clusters have a strong influence on planet formation, mass segregation, and other aspects of star and star cluster evolution. Thus it is important to develop first-principles predictions for the expected motion of newborn stars as a function of parent cloud properties, a subject that has not been investigated prior to this work. In this Letter, we use turbulent, self-gravitating, adaptive mesh refinement (AMR) hydrodynamics simulations of star cluster formation to study the properties of the star clusters that form in clouds in a variety of dynamical states.  In particular, we examine the stellar velocity dispersion, projected velocities, and velocity gradients. We compare the simulation results with recent observations of the ONC and discuss the relevance to current modeling of star clusters. In \S \ref{parameters} we describe our simulation methods and parameters. In  \S \ref{results} we present our results, and we summarize our findings in \S \ref{conclusions}.

\section{SIMULATION PARAMETERS} \label{parameters}

For our study, we perform three simulations of an isothermal, non-magnetized molecular cloud using the ORION AMR code (e.g.~\citealt{truelove98, klein99}). We reexamine the driven and decaying turbulent simulations described in \citet{Offner08a}, and we perform an additional simulation with higher resolution and different turbulent conditions.

We generate the initial state of the runs by driving turbulent motions in the gas for two crossing times without gravity. At $t=0$, a turbulent steady state has been achieved and the gravitational energy is comparable to the turbulent energy in the box. The cloud virial parameter is then
\begin{equation}
\alpha = \frac{5 \sigma^2R}{GM } = 1.67
\end{equation}
\citep{bertoldi92}.
In the simulation we refer to as U1 energy injection is halted after the initial driving phase and the turbulence decays
on a dynamical timescale, initiating a global collapse of the type posited by \citet{tobin09}.
In the simulation D1 we continue turbulent driving so that the Mach number and virial parameter are approximately constant;
this is a very rough representation of the behavior expected in a cloud where feedback drives turbulence and prevents global collapse \citep[e.g.][]{li06, nakamura07}.
The third simulation, D2, begins with the same physical initial conditions, but we adopt a driving pattern with $k=1..2$ and drive the turbulence at a constant energy injection rate. As gravitational motions become significant, the net kinetic energy in the box increases and the virial parameter rises.
The parameters for the three runs are summarized in Table \ref{table1}.

Isothermal self-gravitating gas is scale free, so it is possible to scale the simulation values to observed cloud properties using the characteristic Jeans length and Mach number (see scaling relations in \citealt{Offner08b}). However, the main quantity of interest, the ratio of the stellar dispersion to the gas dispersion, is independent of the physical values we adopt. To estimate physical stellar velocities, we adopt a gas temperature of $T=20$ K and  a gas density of $\rho = 9.74\times 10^{-21}$ g cm$^3$, which are similar to the conditions observed in Orion \citep{johnstone01, johnstone06}. This corresponds to a box length of 2 pc and a cloud mass of 1150 $\ifmmode {\rm M}_{\mathord\odot}\else $M_{\mathord\odot}$\fi$.

We run the simulations with gravity for
one mean-density
freefall time, $t_{\rm ff}$.  We introduce a sink particle whenever the Jeans conditions is exceeded on the finest AMR level \citep{krumholz04}. These particles represent stars, serving as an indicator of the stellar cluster velocity dispersion  imparted by the initial turbulent conditions.

Our simulations neglect radiative feedback, which likely has negligible influence on the initial stellar velocity. The simulations also neglect magnetic fields, which may affect the details of the simulated turbulence and exact distribution of gas densities. However, strong magnetic fields are more likely to further constrain the motion of the dense gas and forming stars rather than enhance it (relative to non-magnetized gas).

\begin{deluxetable}{ccccccc}
\tablewidth{0pt}
\tablecolumns{7}
\tablecaption{Simulation Parameters\tablenotemark{a} \label{table1}}
\tablehead{
\colhead{} & \colhead{$N_0^3$\tablenotemark{b}} & \colhead{$\Delta x_{\rm min}$\tablenotemark{c}} & \colhead{Driving Method}  & \colhead{$k_{\rm min}..k_{\rm max}$}& ${\mathcal M}_f$\tablenotemark{d} & $\alpha_f$\tablenotemark{e}}
\startdata
U1 & 128$^3$	& 200 & No Driving  & 3..4 & 3.9 &  1.1 \\
D1 & 128$^3$	& 200  & Constant $\mathcal{M}$ & 3..4  & 8.5  & 1.7 \\
D2 & 512$^3$	& 100  &  Constant $dE/dt$ & 1..2 & 12.0& 3.4  \\
\enddata
\tablenotetext{a}{All simulations have the same initial $L_0$, $\rho_0$, $T_0$, ${\mathcal M}$. }
\tablenotetext{b} {Number of cells on the base grid.}
\tablenotetext{c} {Minimum cell size in AU.}
\tablenotetext{d} {Final 3D Mach number at $1t_{\rm ff}$.}
\tablenotetext{e} {Final virial parameter at $1t_{\rm ff}$.}
\end{deluxetable}

\section{RESULTS} \label{results}

In this section, we analyze the distribution and velocity information of the stars that form in each simulation. The simulation with decaying turbulence, U1, produces only 20 stars. D2 with the highest amount of turbulence and highest resolution forms 74 stars. In our simulations we do not resolve close binaries, and so in some cases the stellar properties represent the net mass and momentum of a stellar system. At the final time, the runs have efficiencies of $\sim$ 10-15 \%.
Because we do not include protostellar winds, we neglect ejection of mass from prestellar cores, and thus likely overestimate both the stellar masses and the overall efficiency by a factor of
$3-4$
\citep{enoch08, rathborne09}. However, this shift in stars' masses is unlikely to affect their kinematics.

\subsection{Stellar Velocity Dispersion}

Figure \ref{sigmavst} shows the 3D stellar velocity dispersion as a function of time.
The larger driving scale in D2 works to delay the onset of collapse so that stars form much later  than in the D1 and U1 calculations.
For U1 and D1, the stellar velocity dispersion is sonic to transonic. Once a sufficient number of stars are formed, the velocity dispersion trend appears to be mostly  uncorrelated with the simulation Mach number. In the undriven case, the stellar velocity dispersion declines slightly, while remaining nearly fixed in the simulation with constant Mach number driving. In D2, the dispersion decreases slightly with increasing Mach number.
However, in all cases the stellar dispersion is significantly less than the global gas velocity dispersion:
$\sigma_*/\sigma_{\rm gas} \simeq 0.2$. This corresponds to final virial parameters of
$\sim$ 0.10, 0.09 and 0.10 for U1, D1, and D2, respectively, using the total mass in stars and gas. If all the remaining gas were to disperse, leaving only stellar mass, then the virial parameter of the stars would be measured to be $\sim$ 0.8, 0.6, and 1.1 for the runs U1, D1, and D2, respectively.

This indicates that only if the initial gas cloud is significantly supervirial does the virial parameter of the stellar component become greater than unity, and that star clusters forming in turbulent virialized clouds naturally begin with subvirial velocities. Consequently, the subvirial nature of the stellar velocity dispersion is not necessarily an indication of collapse.
Indeed, \citet{proszkow09} find that they are able to reproduce the kinematic signature of Orion assuming an initial virial ratio of $\sim 0.1$ (for stars plus gas), the value produced by simulation D1, in which the gas cloud is turbulent and virialized. Similarly, \citet{allison09}  assume an initial cluster virial parameter of 0.8, much less subvirial than we find for our simulated clouds with virial parameters of $\alpha \simeq 1$.
If the ratio of the stellar dispersion to the gas dispersion influences how quickly the cluster undergoes mass segregation, then most clusters with substantially subvirial initial stellar velocities dispersions may segregate on short timescales.

\begin{figure}
\plotone{f1.eps}
\caption{The top plot shows the 3D star velocity dispersion as a function of time, while the bottom plot shows the dispersion normalized to the gas velocity dispersion.
\label{sigmavst}}
\end{figure}

\subsection{Projected Stellar Velocities }

In this section we investigate the stellar velocities projected along the line-of-sight to the cluster. For the most detailed examination, we use the final time of the D2 simulation, which has the best statistics. Figure \ref{plotvelbin} shows the velocities projected along the three cardinal sight-lines.  The velocities are somewhat clustered with only a couple instances of large deviations.  We have further binned the projected velocities into groups containing 9 stars as a function of position and we have overlaid error bars indicating the velocity dispersion in each bin. These plots look similar to Figure 9 in \citet{proszkow09} of the ONC, which shows evidence of a velocity gradient across the cluster.

We derive the stellar velocity gradient by fitting the line-of-sight velocities with the function $V = V_0 + \nabla V_x X+  \nabla V_y Y$.
Although the data may not be exactly linear in nature (e.g. the ONC exhibits a linear gradient mainly in the cluster center), we can characterize the change in the velocities across the cloud.
Since we use periodic boundary conditions the velocity bins tend to wrap around the domain. Thus, we can better compare to the data by dividing each projection into smaller boxes and investigating the velocity trend in each box separately.
In Figure \ref{grad}, we show the velocity gradients, $|\nabla V| = (\nabla V_x ^2 + \nabla V_y ^2)^{1/2}$, along each projection for the entire domain and for boxes half and a quarter of the domain. We include gradients at three different times and plot ONC data from \citet{tobin09} fitted in the same manner. To limit statistical fluctuations, we only include gradients fitted to 4 or more stars.
The full and half box fits have median velocity gradient magnitudes of 1 km s$^{-1}$ pc$^{-1}$ and 2.1 km s$^{-1}$ pc$^{-1}$, respectively. The magnitude of the gradients increases in both the simulations and the ONC for smaller boxes due to smaller statistics, where stars separated by several shock fronts may easily acquire significantly different velocities.
For comparison,  \citet{proszkow09} find that the ONC has a radial velocity gradient of 1.1 km s$^{-1}$ pc over half of its length, which corresponds to an overall gradient of $\sim$ 0.5 km s$^{-1}$ pc.
This value is similar to what we find for our simulated virialized clusters. Thus velocity gradients are not necessarily a hallmark of a subvirial turbulence or a strongly collapsing cloud.
Instead, large-scale velocity gradients appear to be a signature of the fact that both our simulations and real molecular clouds have most of their turbulent power in large-wavelength modes.

We find that the linear fits of the projected velocities fluctuate somewhat as a function of time, but they tend to fall within 0.2-2 km s$^{-1}$ pc$^{-1}$ on the largest scales.
These values are similar to gradients in the gas centroid velocity observed in regions like Perseus and Ophiuchus \citep{andre07, roso08}.

\subsection{Gas Velocities }

Observations of clusters frequently find that the star and gas velocities are correlated (e.g., \citealt{furesz08}). We can reproduce such observations by making simulated observations of the gas in particular molecular tracers. Following the procedure of \citet{Offner08b}, we generate a position-position-velocity cube of the gas intensity in the $^{13}$CO (J = 1 $\rightarrow$ 0) molecular line along one line-of-sight in the D2 run. We solve for the line emission, assuming that the gas is in statistical equilibrium and that radiative pumping is negligible.\footnotemark
\footnotetext{We obtain the molecular data from the Leiden Atomic and Molecular database \citep{schoier05}.}
We adopt a cloud distance of 400 pc and 26'' beam size. To model the telescope resolution, we smear each velocity channel with the Gaussian beam.

Figure \ref{gas_kin} shows the line-of-sight $^{13}$CO gas and star velocities;
following Furesz et al.~(2008, their figure 10), to produce this
figure from the 3-dimensional position-position-velocity cube, we have
integrated along one direction on the plane of the sky, so the
intensity shown in a given $(x, v)$ pixel is the sum over all $y$
values for that $x$ and $v$. Our simulated image shows a strong
correlation between gas and star velocities, similar to that seen in
the star and gas velocity maps presented by \citet{furesz08}. For
the projection shown in Figure \ref{gas_kin}, the mean stellar velocity is $\sim
-0.8$ km s$^{-1}$, somewhat to the left of the mean $^{13}$CO gas
velocity at 0.003 km s$^{-1}$. Most of the disagreement in the
histogram occurs because not all regions traced by $^{13}$CO have
formed stars. Projections along other directions show similar
correlations between gas and star velocities.

The strong correlation between gas and star velocities might seem
surprising, given that the gas and stars have such different velocity
dispersions. The explanation is that the anti-correlation between
density and velocity found in turbulent gas serves to artificially
enhance the gas-star correlation. Figure \ref{gas_kin} is constructed by summing
over the $y$ direction, so each pixel represents a density-weighted
average velocity. The dispersion, however, is a density-weighted
average of velocity squared, which, since velocity and density are
anti-correlated, is much more heavily weighted toward low density,
high velocity gas than what is shown in Figure \ref{gas_kin}. This low density gas
is also that which is least correlated with the stars. The bias is
exacerbated because $^{13}$CO is thermally excited only in relatively
dense gas that is closely associated with forming stars. It is
subthermally excited in the lower density gas that dominates the
dispersion. Thus, both averaging and excitation serve to improve the
correlation between $^{13}$CO emission and stellar velocities, while
at the same time obscuring the low density gas that carries the bulk
of the kinetic energy.

\begin{figure}
\plotone{f2}
\caption{The figure shows the projected line-of-sight velocities in the three cardinal directions as a function of linear position in the sky for the D2 run at $1t_{\rm ff}$. The error bars show the standard deviation of the velocities binned with a constant number of 10 per bin.
\label{plotvelbin}}
\end{figure}

\begin{figure}
\plotone{f3C}
\caption{The figure shows the distributions of gradients, $|\nabla V|$,  for sight lines along the x, y, and z directions as a function of box size. The ONC stellar data has been provided by \citet{tobin09}. The different datasets are offset slightly for clarity.
\label{grad}}
\end{figure}

\begin{figure}
\plotone{f4}
\caption{The figure shows the position of the stars as a function of the line-of-sight velocity (crosses) in run D2. The gas velocity is overlaid, where the colorbar indicates the
relative  $^{13}$CO intensity averaged over the y coordinate.
 We have assumed that cloud is at a distance of 400 pc and is observed by a 26'' beam. The histograms on the right shows the distribution of gas (solid) and star (dashed) velocities normalized to the total intensity and star number, respectively.
\label{gas_kin}}
\end{figure}

\section{CONCLUSIONS} \label{conclusions}

In this letter, we have investigated the kinematic properties of stars forming in simulated  turbulent self-gravitating clouds. These calculations represent the initial conditions of stars in clusters prior to dynamic evolution and cloud dispersal. We demonstrate that subvirial stellar velocity dispersions arise naturally from clouds in virial equilibrium. Thus, long lived clouds satisfying a quasi-steady state are not ruled out by observations of clusters with small stellar velocity dispersions. We find that gas virial parameters of order unity produce a star to gas velocity dispersion ratio of approximately 0.2.

We find that turbulent initial conditions easily generate stellar velocity gradients of the magnitude observed in star-forming regions.
Using a planar fit,
we demonstrate typical gradients of 1 km s$^{-1}$ pc$^{-1}$, similar to the observed gradient in the ONC.
We also show that for young clusters there is a strong similarity in the simulations between the dense gas traced by $^{13}$CO and the star velocities.
Thus, trends in the stellar velocities as a function of position may be indicative of
the dominance of large-scale turbulent modes,
rather than evidence of global collapse or cloud rotation as suggested by some authors.

 \acknowledgments{
 Support for this work was provided by the US Department of Energy at the Lawrence Livermore National Laboratory under contract  B-542762 (S.~S.~R.~O., C.~E.~H.),
the Alfred P.\ Sloan Foundation (M. R. K.), NASA/JPL through the Spitzer Theoretical Research Program (M. R. K.), the National Science Foundation through grants AST-0807739 (M.~R.~K.) and AST-0901055(S.~S.~R.~O). Computational resources were provided by the NSF San Diego Supercomputing Center through NPACI program grant UCB267.}

redpajama_set_name: RedPajamaArXiv


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

    %

    %
    %

    %

    %

Recurrent neural networks (RNN) and architectures based on RNNs like LSTM ~\cite{Hochreiter:97} has been used to process sequential data more than a decade. Recently, alternative architectures such as convolutional networks~\cite{Dauphin:17,Gehring:17} and transformer model ~\cite{Vaswani:17} have been used extensively and achieved the state of the art result in diverse natural language processing (NLP) tasks. Specifically, pre-trained models such as the OpenAI transformer~\cite{Radford:18} and BERT~\cite{Devlin:18} which are based on transformer architecture, have significantly improved accuracy on different benchmarks.

In this paper, we are introducing a new dataset which we call \textit{ParagraphOrdering}, and test the ability of the mentioned models on this newly introduced dataset. We have got inspiration from "Learning and Using the Arrow of Time" paper ~\cite{Wei:18} for defining our task. They sought to understand the arrow of time in the videos; Given ordered frames from the video, whether the video is playing backward or forward. They hypothesized that the deep learning algorithm should have the good grasp of the physics principle (e.g. water flows downward) to be able to predict the frame orders in time.

Getting inspiration from this work, we have defined a similar task in the domain of NLP. Given two paragraphs, whether the second paragraph comes really after the first one or the order has been reversed. It is the way of learning the arrow of times in the stories and can be very beneficial in neural story generation tasks. Moreover, this is a self-supervised task, which means the labels come from the text itself.

\begin{table}[bp]
	\begin{center}
		\begin{tabularx}{\linewidth}{lX}
			\hline Label & 1\\
			\hline First Paragraph & Now they were walking through the trees, one of them carrying him in its huge arms, quite gently. He was scarcely conscious of his surroundings. It was becoming more and more difficult to breathe.\\
			\hline Second Paragraph & Then he felt himself laid down on something soft and dry. The water was not falling on him now. He opened his eyes.\\
			\hline
		\end{tabularx}
	\end{center}
	\caption{\label{Example} A single example of \textit{ParagraphOrdering} dataset.  }
\end{table}

\section{Paragraph Ordering Dataset}

We have prepared a dataset, \textit{ParagraphOrdreing}, which consists of around 300,000 paragraph pairs. We collected our data from Project Gutenberg. We have written an API for gathering and pre-processing in order to have the appropriate format for the defined task.\footnote{API for downloading the dataset:
\url{https://github.com/ShenakhtPajouh/transposition-data}.
The implementation of different algorithms: \url{https://github.com/ShenakhtPajouh/transposition-simple}}
Each example contains two paragraphs and a label which determines whether the second paragraph comes really after the first paragraph (true order with label 1) or the order has been reversed (Table~\ref{Example}). The detailed statistics of the data can be found in Table~\ref{Statistics}.

\begin{table}[h]
	\begin{center}
		\begin{tabular}{lr}
			\hline \#Train Samples & 294265 \\
			\#Test Samples & 32697  \\
			\hline Unique Paragraphs & 239803 \\

			\hline Average Number of Tokens & 160.39 \\
			Average Number of Sentences & 9.31\\
			\hline
		\end{tabular}
	\end{center}
	\caption{\label{Statistics} Statistics of \textit{ParagraphOrdering} dataset.  }
\end{table}

\section{Approach}
\label{sec:approach}

Different approaches have been used to solve this task. The best result belongs to classifying order of paragraphs using pre-trained BERT model. It achieves around $84\%$ accuracy on test set which outperforms other models significantly.

\begin{table}[h]
\begin{center}
\begin{tabular}{|l|r|}
\hline \bf Model & \bf Accuracy ($\pm0.01$) \\ \hline
LSTM+Feed-Forward & 0.518  \\
LSTM+Gated CNN+Feed-Forward & 0.524  \\
BERT Features(512 tokens)+Feed-Forward & 0.639  \\
BERT Classifier(30 tokens / 15 tokens from each paragraph) & 0.681 \\
BERT Classifier(128 tokens / 64 tokens from each paragraph) & 0.717 \\
BERT Classifier(256 tokens / 128 tokens from each paragraph) & 0.843 \\
\hline
\end{tabular}
\end{center}
\caption{\label{Accuracy} Accuracy on Test set.  }
\end{table}

\subsection{Encoding with LSTM and Gated CNN}
\label{sect:lstmgated}

In this method, paragraphs are encoded separately, and the concatenation of the resulted encoding is going through the classifier. First, each paragraph is encoded with LSTM. The hidden state at the end of each sentence is extracted, and the resulting matrix is going through gated CNN~\cite{Dauphin:17} for extraction of single encoding for each paragraph. The accuracy is barely above $50\%$, which depicts that this method is not very promising.

\subsection{Fine-tuning BERT}
\label{ssec:BERT}
We have used a pre-trained BERT in two different ways. First, as a feature extractor without fine-tuning, and second, by fine-tuning the weights during training. The classification is completely based on the BERT paper, i.e., we represent the first and second paragraph as a single packed sequence, with the first paragraph using the A embedding and the second paragraph using the B embedding. In the case of feature extraction, the network weights freeze and CLS token are fed to the classifier. In the case of fine-tuning, we have used different numbers for maximum sequence length to test the capability of BERT in this task. First, just the last sentence of the first paragraph and the beginning sentence of the second paragraph has been used for classification. We wanted to know whether two sentences are enough for ordering classification or not.
After that, we increased the number of tokens and accuracy respectively increases. We found this method very promising and the accuracy significantly increases with respect to previous methods (Table~\ref{Accuracy}). This result reveals fine-tuning pre-trained BERT can approximately learn the order of the paragraphs and arrow of the time in the stories.


text:
\section{Following Gibbs states}
{How to follow adiabatically a given Gibbs state? Consider
  first the two ``up'' and ``down'' equilibrium states in an Ising
  ferromagnet at low temperature. We can force the system to be in the
  Gibbs state of choice by fixing the all negative or all positive
  boundary conditions. Even far away from the boundaries, the system
  will stay in the selected state for all $T<T_c$ (above the Curie
  point $T_c$ any boundary condition will result in a trivial
  paramagnetic state). By solving the thermodynamics conditioned to
  the boundaries, we can thus obtain the adiabatic evolution of each of
  the two states.

  What boundary conditions should be applied in glassy systems where
  the structure of Gibbs states is very complicated? The answer is
  provided by the following gedanken experiment \cite{BB}: Consider an
  equilibrium configuration of the system at temperature~$T_p$. Now
  freeze the whole system {\it except} a large hole in it. This hole
  is now a sub-system with a boundary condition {\it typical} for
  temperature~$T_p$. If the system is in a well-defined state, then no
  matter the size of the hole, it will always remain correlated to the
  boundaries and stay in the same state. One may now change the
  temperature and study the adiabatic evolution of this state. This
  can be achieved through Monte-Carlo simulations in any model. We
  shall here instead concentrate on mean-field systems, where this
  construction allows for an analytic treatment in the spirit of the
  cavity method \cite{Reconstruction,PLANTING}.  }

\section{Mean field glassy models}
We focus on two of the most studied mean field glassy systems: the
{Ising} $p$-spin \cite{p-spin} and the Potts glass \cite{POTTSGLASS}
models. These are also of fundamental importance in computer science
where they are known as the XOR-SAT \cite{XORSAT} and coloring
\cite{COLORING} problems. Consider a graph defined by its
vertices $i=\{1,\dots,N\}$ and edges $(i,j)\in {\cal E}$, the coloring
Hamiltonian reads
\be
{\cal H}(\{s\}) = \sum_{(i,j) \in {\cal E}} \delta(s_i,s_j)\, ,
\label{H_col}
\end{equation}} \newcommand{\bea}{\begin{eqnarray}
where $s=1,\ldots,q$ are the values of the Potts spins.  The
{Ising} $p$-spin is defined on a hyper-graph with $N$
vertices and $M$ $p$-body interactions (or constraints, if only zero
energy configurations are of interest) with the Hamiltonian
\be {\cal H}(\{s\}) = - \sum_{a=1}^M J_a \prod_{i\in \partial a} s_i\, ,
\label{H_pspin}
\end{equation}} \newcommand{\bea}{\begin{eqnarray}
where $s=\pm1$, $\partial a$ is the set spins involved in interaction
$a$, and $J_a=\pm 1$ are chosen uniformly at random. For XOR-SAT, one
defines instead the number of unsatisfied constraints $E_{\rm
  xor}=(M+E_{\rm p-spin})/2$, hence for XOR-SAT also the temperature
is divided by a factor $2$ with respect to the $p$-spin model. With
these definitions, the XOR-SAT and coloring problems are said to be
satisfied if the ground state energy is zero. In both cases, we will
consider the lattice to be a random (hyper-)graph with fixed degree~{$c$},
i.e. every variable being involved in $c$ interactions, in the
thermodynamic limit, $N \to \infty$. We also consider the large
connectivity $c\to N^{p-1}/(p-1)!$ limit of (\ref{H_pspin}) with $J_a
= \pm \sqrt{p!}/(\sqrt{2N^{p-1}})$: the ``fully-connected'' $p$-spin
model \cite{p-spin}.

\section{Cavity equations for following states}
{We now derive the state-following equations for $T_p\ge T_K$
  in the XOR-SAT problem. This derivation can be generalized for
  $T_p<T_K$, and for any model where the cavity approach \cite{cavity}
  can be applied, and this shall be detailed elsewhere
  \cite{longversion}.}

As a large random (hyper-)graph is locally tree-like, let us thus
first consider the problem on a large (hyper-)tree (see
Fig.~\ref{fig0}). Once the proper boundary conditions are chosen,
computations on the tree give correct results for large random graphs:
this is the basis of the cavity approach.

\begin{figure}[t]
\begin{center}
\includegraphics[width=8.5cm]{tree_new.eps}
\end{center}
\caption{(color online) a) Recursive construction of an equilibrium
  configuration at temperature~$T_p$ in XOR-SAT. Given the tree, and a
  random choice of interactions (full square $J=1$, empty square
  $J=-1$), one starts from the root, and chooses iteratively the
  configuration of ancestors (full circles $s=1$, empty circle $s=-1$)
  randomly such that it satisfies the constraints with probability
  $1-\epsilon(T_p)$ (here $\epsilon=3/7$). Violated constraints have
  dashed/red borders.  b) The problem can also be Gauge transformed
  into a fully polarized configuration with
  all $s=1$ but where the $J$'s are chosen from $P_p(J)=\epsilon(T_p)
  \delta (J+1) + [1-\epsilon(T_p)] \delta (J-1)$.
\label{fig0}}
\end{figure}

We derive the method of following states for the p-spin
Hamiltonian (\ref{H_pspin}).  We concentrate on an equilibrium
state at $T_p$, with $T_p\!  \ge\!  T_K$, where the equilibrium
solution is given by the replica symmetric (RS) cavity method, or
equivalently by the fixed point of the Bethe-Peierls (or Belief
Propagation, BP) recursion \cite{XORSAT}:
\begin{align}
\beta h^{i\to a} &= \sum_{b\in \partial i \setminus a} {\rm
  atanh}{\big[ \tanh{(\beta J_b)} \prod_{j\in \partial b\setminus i}
  \tanh{(\beta h^{j\to b})} \big]}{ \equiv} \nonumber \\
&{\equiv \beta {\cal F}(\{h^{j \to b}\},\beta)} \, ,
\label{eq_BP}
\end{align}
where $h^{i \to a}$ is an effective cavity field seen by the spin $i$ due to all
its neighboring spins $j$ except those connected to the interaction $a$; by $\partial i \setminus a$ we denote all interaction in which spin $i$ is involved except $a$.
For $T_p\!\ge\!T_K$ the effective field $h^{i\to a}\!=\!0$ for all edges $ia$, and the fraction of violated constraints is hence
$\epsilon(T_p)=(1+e^{2/T_p})^{-1}$ \cite{XORSAT}.

One can generate an equilibrium configuration on a tree once the
effective fields are known{\cite{Reconstruction} using the
  iterative procedure} described in Fig.~\ref{fig0}. The values
assigned to the variables on the leaves are then fixed and the measure
they induce in the bulk of the tree defines {an} equilibrium
Gibbs state at temperature~$T_p$. As long as $T_p\ge T_K$, the
computation in the bulk of the tree describes correctly the properties
of the problem on the random graph.

All the properties of the Gibbs state are thus obtained by solving the
BP equations initialized  {in this} boundary
configuration. When $T=T_p\ge T_K$, the results of the usual 1RSB
calculation are recovered, as first discussed in the context of
reconstruction on trees \cite{Reconstruction}. For $T=T_p>T_d$ BP will
converge back to $h^{i\to a}\!=\!0$ for all edges $ia$, but when $T_d
\ge T=T_p\ge T_K$, the configuration we picked lies in one of the
exponentially many equilibrium Gibbs states and the BP fixed point
thus describes one of them.

Now is the new crucial turn: Since the boundary conditions define the
equilibrium Gibbs state at $T_p$, we can use the BP equations
(\ref{eq_BP}) initialized in the boundary condition but with a {\it
  different} temperature $\beta=1/T \neq 1/ T_p$. The resulting fixed
point now describes the properties of {\it the very same state} but at
a {\it different temperature} $T\neq T_p$.

This line of reasoning translates readily in a set of coupled
recursive cavity equations, which are a two-temperatures extension
of the reconstruction formalism \cite{Reconstruction,PLANTING}. Two
distributions of fields $P_s(h)$, with $s=\pm 1$ depending on whether
the site was set $\pm 1$ in the broadcasting, are given by
\bea P_s(h) =\sum_{J_a} P(J_a) \sum_{\{s_i\}} \frac{ e^{\beta_p J_a s
    \prod_{i} s_i} }{2^{p-1}\cosh{\beta_p}} \times  \nonumber \\  \int
\prod_{i=1}^{p-1} \prod_{j_i=1}^{c-1} {\rm d}P_{s_i}(h^{j_i}) \,
\delta[h-{\cal F}(\{h^{j_i}\},\beta)]\, , \label{pop_pl} \end{eqnarray}} \newcommand{\mc}{\mathcal
where the delta function ensures that cavity field $h$ is related to
the fields $h^{j_i}$ via eq.~(\ref{eq_BP}).
The term in front of the integral describes the
properties of the equilibrium configuration at inverse temperature
$\beta_p$. Eq.~(\ref{pop_pl}) describes adiabatic evolution of a Gibbs
state that is one of the equilibrium ones at $\beta_p<\beta_K$ in the $p$-spin model with
distribution of interaction strength
$P(J_a)=\rho\delta(J-1)+(1-\rho)\delta (J+1)$, it can be solved
numerically using the population dynamics technique \cite{cavity}. The
(Bethe) free energy of the same state at the new temperature~$T$ reads
\bea
\beta f_{\beta_p}(\beta) \!\!\!\!&&\!\!\!\! = \frac{c-1}{2} \sum_{s_i} \int
\prod_{i=1}^c {\rm d}P_{s_i}(h^i)
\log{Z^{i}} -  \frac{c}{p}   \sum_{J_a} P(J_a) \nonumber \\  \!\!\!\!&&\!\!\!\!
 \sum_{\{s_i\}} \frac{ e^{\beta_p J_a \prod_{i}
    s_i}}{2^p \cosh{\beta_p}} \int\prod_{i=1}^p \prod_{j_i=1}^{c}
{\rm d}P_{s_i}(h^{j_i})
\log{Z^{a+\partial a}} \, , \label{free_pl} \end{eqnarray}} \newcommand{\mc}{\mathcal
where the $Z$'s are
\bea
Z^{a+\partial a}              \!\!\!\!&=&\!\!\!\!  \frac{\cosh(\beta J_a) \prod_{i\in\partial a} 2 \cosh(\beta h^{i\to a})}
{\prod_{i\in\partial a} \prod_{b\in\partial i \setminus a}
2 \cosh(\beta u^{b\to i})}[1+\tanh{(\beta u^a)}]\, , \nonumber \\
Z^{i} \!\!\!\!&=&\!\!\!\!  2 \cosh{(\beta \sum_{b\in \partial i}
  u^{b\to i})} \prod_{{b\in \partial i}} \frac{1}{2\cosh{\beta u^{b\to
      i}}}\, , \end{eqnarray}} \newcommand{\mc}{\mathcal
 where $\tanh{(\beta u^{b\to i})} = \tanh{(\beta
  J_b)} \prod_{j\in\partial b \setminus i} \tanh{(\beta h^{j\to b})}$.

\begin{figure*}[!ht]
\begin{center}
\hspace{-0.4cm}
\includegraphics[width=6.44cm]{PSPIN_FOLLOW_1RSB_letter.eps}
\hspace{-0.75cm}
\includegraphics[width=6.44cm]{FOLLOW_XOR_3_3.eps}
\hspace{-0.75cm}
\includegraphics[width=6.44cm]{FOLLOW_COL_9_4.eps}
\put(-510,80){\begin{sideways} Energy~$e$ \end{sideways}}
\end{center}
\caption{(color online) Behavior of Gibbs states in mean field glassy
  systems: the energy~$e$ is plotted as a function of temperature~$T$.
  The equilibrium energy (blue) is shown, together with the energy of
  several states (red) bellow the dynamical transition where we follow
  them out-of-equilibrium. States can be followed a) upon heating
  until a spinodal point (black-dotted curve) and b) upon cooling
  until they reach their bottom at zero temperature. However, the
  states closer to $T_d$ undergo a FRSB transition upon cooling and
  divide into many marginally stable sub-states at a temperature
  $T_m<T_p$ (red-dotted curve).  Left: the fully-connected 3-spin
  model with its dynamical $T_d$, Kauzmann $T_K$ and Gardner $T_G$
  (equilibrium states becoming FRSB \cite{Gardner}) transitions. We
  have used the 1RSB formalism to follow all the states down to zero
  temperature. The uppermost point at $T=0$ thus gives a 1RSB lower
  bound on the limiting energy of an adiabatically slow simulated
  annealing. Inset: the behavior of states for $T_p<T_K$, explicitly
  demonstrating the presence of temperature chaos. Middle: the XOR-SAT
  problem for $p=3,c=3$ where $T_K=0$. Again, the higher energy states
  become unstable, here we used only the RS formalism and thus were
  not able to continue all of them to zero temperature. In this model
  all $T_p>0$ states finish at a finite bottom energy at $T=0$.
  Right: The $4$-coloring of graphs with degree $c=9$. The situation
  is similar to the XORSAT case except that here all the states have
  their energy decaying very fast to zero when $T\to 0$. {Asterixes
  (green) represent the results of adiabatic simulations starting from
  an equilibrated configuration on a $N=10^5$ graph using
  Monte-Carlo (in XORSAT) and BP (in coloring)
  evolution.} \label{fig1} \vspace{-2mm} }
\end{figure*}

\section{Gauge transformation} In the $p$-spin model, the above
equations can be further simplified by exploiting a Gauge
invariance. For any spin $i$, the Gauge transformation $ s_i \to -s_i$ and
$J_{ij} \to - J_{ij}$ for all $j \in \partial i$
keeps the Hamiltonian eq.~(\ref{H_pspin}) invariant. As shown in
Fig.~{\ref{fig0}}, this allows to transform the equilibrium spin configuration into a uniform one (all $s=1$), the disorder distribution then changes to $P_p(J)$ (see Fig.~{\ref{fig0}}).
Since all $s=1$, there is no need to distinguish between the $+1$ and
the $-1$ sites, and eq.~(\ref{pop_pl}) reduces to the usual replica symmetric cavity equation for a problem with mixed ferromagnetic/anti-ferromagnetic interactions at temperature~$T$ initialized in the uniformly positive state, where the distribution of interactions is given by the
Nishimori-like \cite{Nishimori}
condition $P(J)=\epsilon(T_p) \delta (J+1) + [1-\epsilon(T_p)] \delta (J-1)$.

The Gauge invariance have thus transformed the task of following an
equilibrium state in a glassy model into simply solving a known
ferromagnetically biased model with the standard cavity approach.  One
can show in particular that following states in the fully connected
$p$-spin model for $T_p\ge T_K$ is equivalent to solving the $p$-spin
model with an additional effective ferromagnetic coupling $\langle J_a
\rangle N^{p-1}/p! = J_0^{\rm eff} =\beta_p^2/2$, and one can thus
readily take the solution of the p-spin in the literature,
e.g. \cite{p-spin}, to obtain properties of equilibrium states.

\section{Energy landscape} We now present the results of the above
formalism for the fully connected $3$-spin problem, the XOR-SAT with
$c=3, p=3$ and the $4$-coloring of graphs with degree
$c=9$. Fig.~\ref{fig1} shows the energy density $e(T)$ for several
Gibbs states, that are the equilibrium ones at $T=T_p$, as they become
out-of-equilibrium at $T\neq T_p$. Although such plots were often
presented as sketches in previous works, analytical results were so
far available only for few very simple models
\cite{kurchan,spherical}.

We confirmed that glassy equilibrium Gibbs states exist only for $T_p
\le T_d$ while for $T_p>T_d$ we saw only the liquid solution. As these
states are heated, they can be followed until a well-defined spinodal
temperature~$T_s(T_p)$ that grows as $T_p$ decreases (this
is reminiscent of the Kovacs effect in glassy materials
\cite{kovacs}). Interestingly we find that $T_s(T_p=T_d)=T_d$, i.e. an
equilibrium state at $T_d$ disappears (melts) for any increase of
temperature, this is at variance with the (unphysical) behavior in
spherical models where such a state exists until much larger
$T$~\cite{spherical}.

We also consider the state evolution upon cooling. As anticipated
based on the statistical features of the energy landscape
\cite{ISOCOMPLEXITY} and the study of spherical models
\cite{kurchan,spherical}, we find that for $T_p$ near enough the
dynamical temperature~$T_d$, the states undergo a FRSB transition: at
some $T_m<T_p$ the states decompose into many marginally stable
sub-states \cite{Gardner}. In such a case the exact solution for
adiabatic evolution requires the FRSB approach \cite{MPV} and our RS
approach yields only a lower bound on the true energy. Moreover, for
temperatures slightly below the FRSB instability, the RS solution
undergoes an unphysical spinodal transition and we are thus unable to
obtain even the RS lower bound. For the fully connected $p$-spin model
we have therefore used the 1RSB formalism (using the mapping onto the
ferromagnetic-biased model) which allowed us to follow states until
$T=0$ (see Fig.~\ref{fig1}, right panel).  The FRSB solution is
numerically much more involved, making the exact analysis obviously more
difficult.

Note that following the evolution of a state that is the equilibrium
one at $T_p=T_d$ is particularly interesting. An adiabatically slow
annealing is able to equilibrate down to $T_d$ and the evolution of
the state at $T_d$ is thus giving the asymptotic behavior of the
simulated annealing algorithm.

{To assess the validity of our approach, we have performed
  numerical simulations. For models such as XOR-SAT and coloring for
  $T\ge T_K$, where the annealed average is equal to the quenched
  average, it is possible to use the quiet planting trick
  \cite{PLANTING} to generate an equilibrated configuration together
  with a typical random graph: Starting with a random configuration
  one simply creates the (hyper-)graph randomly such that the
  configuration has $\epsilon(T_p) N$ violated constraints. As seen in
  Fig.~\ref{fig1}, when initialized in $T_p \in [T_d,T_K]$, the
  evolution upon heating and cooling, simulated both by BP and slow
  Monte-Carlo simulations, follows precisely our predictions.}

Finally, we have also consider the adiabatic following of states
for $T_p<T_K$. Their behavior is depicted in the inset of the left
panel in Fig.~\ref{fig1}.  At variance with the situation in spherical
models \cite{spherical}, these states become out-of-equilibrium as
soon as the temperature is changed. The equilibrium configurations for
$T<T_K$ thus {\it do not} belong to a single state, but instead to a
succession of many different states whose free energies cross as $T$
changes: This demonstrates explicitly the presence of temperature
chaos in the static glass phase $T<T_K$ \cite{chaos}.

\section{Comparison with previous approaches}
{How do our results compare with previous heuristic
  approaches to adiabatic annealings? Two arguments have been mainly}
proposed. The first one is the marginality criterion \cite{DYNAMIC}
according to which the energy {reached by a slow annealing at
  zero temperature} can be computed by sampling a typical energy
minima at a given energy~$e$, and then choosing $e$ such that this
minima is marginally stable with respect to the replica symmetry
breaking.  This uniform minima-sampling argument is, however,
unjustified, since the dynamics always goes in out-of-equilibrium
states bellow $T_d$. Indeed, as already shown by \cite{ISOCOMPLEXITY},
the marginality criterion is not correct and its result not
consistent.

A refined approach called iso-complexity was thus introduced in
\cite{ISOCOMPLEXITY}.  {It is} proposed to count the number
of equilibrium states at a given $T_p$, and then to consider the
energies at $T<T_p$ for which the number of states is equal to the one
at $T_p$. Iso-complexity leads indeed to a lower bound on adiabatic
annealings, because in order to end up at lower energies one would
have to be exponentially lucky. {With} our formalism (where
we explicitly follow states) we checked that adiabatic annealings
always ends up at higher energies than the iso-complexity ones. This
is illustrated in Fig.~\ref{fig_iso}, where we show the energy of the
bottoms of states that are the equilibrium ones at temperature~$T_p$
and compare it to the iso-complexity result (that gives strictly lower
values).
Unfortunately the RSB instability within states
mentioned previously prevents us from estimating the asymptotic energy
when $T_p$ is close to $T_d$.

\begin{figure}[!ht]
  \hspace{0.2cm}
  \includegraphics[width=8.2cm]{ISO.eps}
\put(-240,80){\begin{sideways} Energy~$e$ \end{sideways}}
\caption{\label{fig_iso} (color online) Comparison between the exact
  adiabatic evolution of states and the iso-complexity lower
  bound. The black line is the energy of the bottoms of states that
  were the equilibrium ones at temperature $T_p \in [T_K , T_d]$ in
  XOR-SAT for $c=3$, $p=3$ (left, the inset is a zoom) and for $c=4$,
  $p=3$ (right). The red (upper-most) line is the equilibrium energy
  at $T_p$. The blue line is the iso-complexity lower bound
  \cite{ISOCOMPLEXITY}.}
\end{figure}

\section{Canyons versus valleys}
An important class of mean field glassy models are the constraint
satisfaction problems where one searches for a configuration
satisfying all the constraints. As opposed to early predictions
\cite{SP}, it has been observed that glassiness does not prevent
simple algorithms from finding a ground state
\cite{COLORING,Heuristic}. The following state method allows to
understand this fact and to shed light on the energy landscape of
these problems. Indeed, we see in Fig.~\ref{fig1} that although in the
XOR-SAT problem all the typical finite $T_p$ equilibrium states have
their bottoms at positive energies, the situation is different in
$4$-coloring with $c=9$ where all depicted states descent to zero
energy when the temperature is lowered.

\begin{figure}[t]
\includegraphics[width=\linewidth]{PLOT_STATES.eps}
\caption{(color online) Data from Fig.~\ref{fig1} plotted in order to
  visualize the energy landscape. The energy~$e$ is plotted against
  the entropy $s=\beta (e-f)$ for different equilibrium states. We
  plotted $s(e)/2$ and $-s(e)/2$ such that the width corresponds to
  the logarithm of the number of configurations at energy~$e$ for the
  Gibbs state. Left: XOR-SAT with $c=3$, $p=3$.  Right: 4-coloring of
  random graphs with $c=9$. The black curve corresponds to the
  equilibrium total entropy. The red curves are different equilibrium
  states, corresponding to $T_p=0.15,0.2,0.24$ (left), and $T_p=0.12$
  on (right), the energies corresponding to $T_p$ are depicted by
  horizontal black dashed lines. The left side states, with their
  finite energy bottoms, remind us of the valley in Fig. \ref{fig00}, while the right side reminds of deep canyons that all
  reach the ground state energy. \label{fig4} }
\end{figure}

Looking back to Fig.~\ref{fig00}, we see that the landscape has many
valleys with bottoms at finite energies, but also canyon-shaped states
that reach the ground state energy. We now define two types of glassy
landscape: (1) In the canyons-dominated landscape, {a typical
  (equilibrium) state at $T_p=T_d$ has its bottom at zero energy
  while} (2) in the valleys-dominated landscape {a typical
  equilibrium state at $T_d$ has its bottom} at strictly positive
energy.

In order to {\it quantitatively} observe canyons and valleys, we have
computed the "shape" of the states (see Fig.~\ref{fig4}). For XOR-SAT at $c=3$, $p=3$ we observe the
standard picture of {\it valleys} with bottoms at positive energy. In
4-coloring of graph with degree $c=9$, however, the states indeed have
a {\it canyon-like} shape and go down to zero energy. While an
adiabatic annealing would be stuck at finite energy in the first case,
it would instead reach a solution (although not an equilibrium one) in
the second one. This is not to say that $4$-coloring of graphs of
degree $c=9$ is really easy (since we are speaking of an infinitely
slow annealing procedure) but rather to explain based on analytical
calculations why it is sometimes possible to find solutions even in
the clustered glassy phase using simple local search algorithms, as
observed in \cite{COLORING,Heuristic}.

In problems such as graph coloring or satisfiability of Boolean
formulas, there will thus be a sharp transition {(in general
  different from the clustering and the satisfiability transitions)}
as the constraint density is increased, where the energy landscape
changes from canyon-dominated to valley-dominated one: this
{transition} marks the onset of difficulty of the problem for
an ideally slow annealing, and most likely also for other stochastic
local search algorithms. This point can be in principle computed by
the following state formalism, however, the replica-symmetry-breaking
instability discussed above complicates the numerical resolution of
the corresponding equations and we will thus discuss it elsewhere
\cite{longversion}. {We also show explicitely in \cite{longversion} that this transition is upper bounded by the so-called rigidity transition point where frozen variables appear in the equilibrium ground state configurations \cite{COLORING,frozen}. This further supports the conjecture of \cite{COLORING} that solutions with frozen variables are really hard to find.}

{Another observation can be made from Fig.~\ref{fig1}: Even
  if one equilibrates the system at $T=T_d$, the state soon becomes
  unstable towards FRSB upon cooling. Therefore {\it any} cooling
  procedure will end up {\it at best} in far from equilibrium} FRSB
states. This shows how futile are the attempts to study equilibrium
predictions, such as the appearance of clustering or BP fixed points,
starting from solutions obtained by heuristics solvers that performing
a kind of annealing in the landscape.  {Instead typical
  configuration {\it must} be obtained. This can be achieve by
  Monte-Carlo, or using exhaustive search \cite{EXAUST} for small
  instances, or by planting\cite{PLANTING} for larger ones (as we did
  in Fig.\ref{fig1}.)}

\section{Conclusions} We have described how to follow adiabatically
Gibbs states in glassy mean field models, and answered some
long-standing questions on their energy landscape: We have discussed
the residual energy after an adiabatically slow annealing, the
behavior of out-of-equilibrium states, and demonstrated the presence
of temperature chaos. We have also found new features of the energy
landscape, and identified a transition from a canyons-dominated
landscape to a valleys-dominated one.

{The following state method presented here has a wide range
  of applications and we believe that many mean fields model will
  profit from being revisited in these directions.}