\section*{Program Summary}
\noindent
{\em Title of program\/}: {\tt t3evol}\\
{\em Version\/}: 1.0 \\
{\em Catalogue identifier\/}: \\
{\em Program obtainable from\/}: {\tt http://arxiv.org/archive/hep-ph} or its mirrors by downloading the source of this document.\\
{\em Distribution format\/}: tar.gz \\
{\em E-mail\/}: {\tt bjoern.pirnay@physik.uni-r.de} \\
{\em Licensing provisions\/}: GNU General Public License. \\
{\em Computers\/}: all. \\
{\em Operating systems\/}: all. \\
{\em Program language\/}: {C++}. \\
{\em Other programs called\/}: none. \\
{\em External files needed\/}: none. \\
{\em No.\ of lines in distributed program, including test data, etc.\/}: 7018\\
{\em No.\ of bytes in distributed program, including test data, etc.\/}: 269445\\
{\em Keywords\/}: Twist-three correlations, evolution equations.\\
{\em Nature of the physical problem\/}: Solving the evolution equations for the twist-three antiquark-gluon-quark and triple-gluon correlation functions with leading order kernels.\\
{\em Method of solution\/}: Series expansion of the formal solution, each term in this series typically requires several integrations, which are performed numerically on a discretized grid.\\
{\em Restrictions on complexity of the problem\/}: Momentum grids of size $N_p>300$ (see text) are not feasible as they will result in a very long computation time.\\
{\em Memory required to execute\/}: Strongly dependent on the number of grid points; for $N_p=101$ (see text) up to $15\,\text{MB}$. \\
{\em Typical running time\/}: Strongly dependent on the number of grid points; for $N_p=101$ (see text) on a PC with a 3.40 GHz processor the time is $<\unit{15}{\second}$.
\end{abstract}%

\end{frontmatter}%

\section{Introduction}

By virtue of factorization theorems, observables in hadronic processes are typically represented as a convolution of a short distance partonic part with long distance distributions.
The former can be calculated reliably in perturbative Quantum Chromodynamics (QCD), while the latter do not admit such treatment.
Both parts depend on a separation (factorization) scale $\mu$, and although the long distance parton distributions are non-perturbative, their scale dependence can be calculated perturbatively.
The corresponding differential equations, dubbed evolution equations, are indispensable in high energy phenomenology.

A prime example for this concept are the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations for the usual collinear parton distribution functions (PDFs).
Typically one invents a model or a suitable parametrization for the PDFs at some low reference scale and evolves them to the relevant scales at which experimental data is available.
One may fit the parameters such that the data is described sufficiently well.
This approach has become customary for the DGLAP setup, where the evolution is governed by the anomalous dimensions of leading twist operators.
Many programs are available to perform the evolution numerically, see e.g.~\cite{Vogt:2004ns,Cafarella:2003jr,Salam:2008qg,Botje:2010ay}.

However, large single spin asymmetries in polarized reactions have been consistently observed in a large variety of experiments, cf.~\cite{Liang:2000gz} for an overview, which can not be explained by the leading twist PDFs, see~\cite{Liang:2000gz,Anselmino:1994gn,Barone:2001sp}.
One possible explanation of the asymmetries, that goes beyond the simple PDF description, relies on large parton correlations inside the nucleon.
Factorization theorems for this framework have been worked out, cf.~\cite{Efremov:1981sh,Qiu:1991wg}.

These correlations are defined via matrix elements of a set of twist-three operators and the corresponding evolution equations have been derived in leading order of the strong coupling, see~\cite{Bukhvostov:1984as,Bukhvostov:1985rn,Balitsky:1987bk,Braun:2009mi,Ma:2012xn}.
The situation is more complicated compared to the twist-two sector, since in general one has to deal with functions of two (or more) momentum fractions.
There seems to be a tendency to consider only certain values of the momentum fractions, like the soft-gluon-pole limit~\cite{Kang:2008ey,Vogelsang:2009pj,Schafer:2012ra,Kang:2012em}.
However the corresponding evolution equations are incomplete.
For a closed set of equations all ingredients (evolution kernels and correlation functions) have to be given as functions of arbitrary momenta.
Technically, the solution of the evolution equations gets more involved compared to e.g.\ the DGLAP-scenario.
On the conceptual level, the invention of a parametrization or a model for all possible momenta is rather in its infancy.
First steps have been made in~\cite{Braun:2011aw}.
Consequently, many existing data analyses either do not include evolution at all or revert to a simplified leading twist evolution.
Either of these treatments lack theoretical justification.
The general case will be addressed here.

In this work a computer program, named {\tt t3evol}, is presented that may either be used directly in the data analysis or serve as a guideline toward a correct treatment of the QCD evolution.
The presentation is organized as follows.
Sec.~\ref{sec:formulation} recapitulates the structure of the twist-three evolution and introduces necessary notations.
The method of solution is discussed in Sec.~\ref{sec:method}.
The actual program, its configuration and certain possible modifications to it are described in Sec.~\ref{sec:description}.
Sec.~\ref{sec:output} is devoted to the explanation of the output generated by {\tt t3evol}.
Conclusions are drawn in Sec.~\ref{sec:conclusion}.

\section{Formulation of the problem\label{sec:formulation}}

The evolution equations for the twist-three correlation functions have been derived in~\cite{Bukhvostov:1984as,Bukhvostov:1985rn,Balitsky:1987bk} and reexamined in~\cite{Kang:2008ey,Braun:2009mi,Schafer:2012ra,Ma:2012xn,Kang:2012em} in different operator bases.
In one of the choices one considers operators
\begin{align}
 S^\pm_{q,\sigma}(z) &=  \bar q(z_1) t^a\bigl(i\tilde F_{\sigma+}^a(z_2)\pm F_{\sigma+}^a(z_2)\gamma_5\bigr)\gamma_+ q(z_3)\,,\notag\\
 F^\pm(z) &= 2g C^{abc}_\pm s^\sigma_T F_+^{\,\,\nu,a}(z_1)F_{+\sigma}^{b}(z_2)F_{+\nu}^{c}(z_3)\,.
\end{align}
Here, the subscript ``$+$'' stands for the contraction of a Lorentz index with a light-like vector $n$.
All fields are assumed to ``live'' on the light-ray, $\phi(z_i) \equiv \phi(nz_i)$.
The color tensors are defined as $C^{abc}_+ = if^{abc}$ and $C^{abc}_- = d^{abc}$.
The vector $s_T^\sigma = \varepsilon^{\sigma\mu\nu\lambda}s_\mu n_\nu \bar n_\lambda$\footnote{The following convention is used: $\varepsilon^{0123}=-1$.} is a projection of the nucleon spin $s$ ($s^2=-1$) transverse to the light-cone spanned by $n$ and a second light-like vector $\bar n$.

For the purpose of this work it is convenient to use $C$-even and $C$-odd operators:
\begin{align}
 \mathbb S^\pm_q(z) &= i s^\sigma_T \bigl(S^+_{q,\sigma}(z) \pm  P_{13}S^-_{q,\sigma}(z) \bigr)\,,\notag\\
 \mathbb F^\pm(z) &= (1 \mp P_{23} \pm P_{12})F^\pm(z)\,,
\end{align}
where $P_{ij}$ denotes the permutation operator acting on the position $z_i$ and $z_j$.
The matrix elements between two nucleon states define the twist-three correlation functions $\mathfrak{S}^\pm_q $ and $\mathcal{F}^\pm$ in momentum fraction space:
\begin{align}
 \bra{p,s_T} \mathbb S^\pm_q(z)\ket{p,s_T} &= 2p_+^2\int\mathcal{D}x\,e^{-ip_+\sum_kx_kz_k} \mathfrak{S}^\pm_q(x)\,,\notag\\
 \bra{p,s_T} \mathbb F^\pm_q(z)\ket{p,s_T} &= 2p_+^3\int\mathcal{D}x\,e^{-ip_+\sum_kx_kz_k} \mathcal{F}^\pm(x)\,,
\end{align}
where the integration measure $\mathcal{D}x$ is given by
\begin{align}
 \mathcal{D}x &= dx_1dx_2dx_3\,\delta(x_1+x_2+x_3)\,.\label{eq:intmeasure}
\end{align}
The definitions of $\mathfrak{S}^\pm_q (x)\equiv \mathfrak{S}^\pm_q (x_1,x_2,x_3)$ and $\mathcal{F}^\pm(x)\equiv \mathcal{F}^\pm(x_1,x_2,x_3)$ coincide with those of Ref.~\cite{Braun:2009mi}.

In the following we will formulate the problem in terms of $\mathfrak{S}^\pm_q$ and $\mathcal{F}^\pm$, since the mixing  under renormalization between the $\mathfrak S$-type and $\mathcal F$-type functions takes a rather compact form.
Generally the notation of~\cite{Braun:2009mi} is used and the conversion to other notations can be achieved by taking an appropriate linear combination.
From now on it is tacitly assumed that all functions correspond to the matrix elements of operators between two \emph{proton} states.

$\mathfrak{S}^\pm_q$ and $\mathcal{F}^\pm$ are real and effectively functions of two momentum fractions, see Eq.~\eqref{eq:intmeasure}, and the scale $\mu^2$.
They obey coupled sets of evolution equations
\begin{align}
 \frac{\partial}{\partial\ln\mu^2}\begin{pmatrix}
                                   \mathfrak{S}^\pm\\
                                   \mathcal{F}^\pm
                                  \end{pmatrix}
&=
-\frac{\alpha_s(\mu^2)}{2\pi}\begin{pmatrix}
                              \mathbb{H}_{QQ}^\pm & \mathbb{H}_{QF}^\pm\\
                              \mathbb{H}_{FQ}^\pm & \mathbb{H}_{FF}^\pm
                             \end{pmatrix}
\begin{pmatrix}
                                   \mathfrak{S}^\pm\\
                                   \mathcal{F}^\pm
                                  \end{pmatrix}\,.\label{eq:evolutioneq}
\end{align}
Note that Eq.~\eqref{eq:evolutioneq} is written for the flavor singlet distribution $\mathfrak{S}^\pm = \mathfrak{S}^\pm_u+\mathfrak{S}^\pm_d$.
The corresponding flavor-nonsinglet distribution obeys a simpler equation, which decouples from the gluonic distributions.
Without loss of generality we assume the evolution to have the structure of Eq.~\eqref{eq:evolutioneq}.
The program will evolve both flavor sectors separately to disentangle up- and down-quarks.
The only quantity on the rhs.\ of Eq.~\eqref{eq:evolutioneq} that depends explicitly on $\mu^2$ is the running coupling $\alpha_s$.
At leading order it is given by
\begin{align}
 \alpha_s(\mu^2) &= \frac{4\pi}{b_0\ln(\mu^2/\Lambda_{\text{QCD}}^2)}\,,\label{eq:runningalphas}
\end{align}
where $b_0 = 11/3\,N_c-2/3\, n_f$ for QCD with $N_c$ colors and $n_f$ flavors. $\Lambda_{\text{QCD}}$ is the dimensional transmutation scale.
By changing variables from $\mu^2$ to the dimensionless \emph{evolution time} $t$,
\begin{align}
 t &= -\frac{2}{b_0}\ln\biggl(\frac{\alpha_s(\mu^2)}{\alpha_s(\mu_0^2)}\biggr)\label{eq:evoltime}\,,
\end{align}
where $\mu_0$ is some reference scale, one can recast Eq.~\eqref{eq:evolutioneq} into a more convenient form,
\begin{align}
 \biggl(\frac{\partial}{\partial t}+\mathbb{H}^\pm\biggr) \begin{pmatrix}
                                   \mathfrak{S}^\pm\\
                                   \mathcal{F}^\pm
                                  \end{pmatrix}
                                  &= 0\,, \label{eq:evolutioneq2}
\end{align}
where
\begin{align}
 \mathbb{H}^\pm &= \begin{pmatrix}
                              \mathbb{H}_{QQ}^\pm & \mathbb{H}_{QF}^\pm\\
                              \mathbb{H}_{FQ}^\pm & \mathbb{H}_{FF}^\pm
                             \end{pmatrix}\,.
\end{align}
For illustrative purposes we restrict ourselves to $\alpha_s$ at leading order.

Suppose at some initial time $t_0=0$ (corresponding to the initial scale $\mu_0$) one specifies the functions $\mathfrak{S}^\pm$ and $\mathcal{F}^\pm$, denoted by $\mathfrak{S}^\pm_0 = \mathfrak{S}^\pm(t=0),\mathcal{F}^\pm_0 = \mathcal{F}^\pm(t=0)$.
Then the formal solution of Eq.~\eqref{eq:evolutioneq2} will be given by exponentiation of the evolution kernels,
\begin{align}
 \begin{pmatrix}
  \mathfrak{S}^\pm(t)\\
  \mathcal{F}^\pm(t)
 \end{pmatrix}
 &=
 \exp\bigl(-t\,\mathbb{H}^\pm\bigr)\begin{pmatrix}
  \mathfrak{S}^\pm_0\\
  \mathcal{F}^\pm_0
 \end{pmatrix}\,.\label{eq:solution}
\end{align}
This expression will be the starting point for the numerical treatment.
Exact approaches like diagonalizing $\mathbb{H}^\pm$ will not be pursued here.
At best, methods analogous to~\cite{Derkachov:1999ze} may only be applicable in certain limits, e.g. at large-$N_c$.

\section{Method of solution\label{sec:method}}

For small deviations of $t$ from the initial $t_0$ the series expansion of the solution in Eq.~\eqref{eq:solution} can be reliably truncated at some finite order.
The remaining ingredient and most expensive operation is the action of a generic ``Hamiltonian'' $\mathbb{H}$.
We give a brief and schematic description of its evaluation below.

Due to translational invariance the momentum fraction arguments $(x_1,x_2,x_3)$ of the twist-three distributions $\mathfrak{S}^\pm,\mathcal{F}^\pm$ are subject to the constraint
\begin{align}
 x_1+x_2+x_3 &= 0\,.
\end{align}
For definiteness, we will eliminate $x_3$ in favor of $x_1,x_2$ and work with a compact notation:
\begin{align}
 \mathfrak{S}^\pm (x_1,x_2) &\equiv \mathfrak{S}^\pm(x_1,x_2,-x_1-x_2)\,, \label{eq:fx1x2}
\end{align}
and similar for $\mathcal{F}^\pm$.

In general $\mathbb{H}$ is a linear one-dimensional integral operator whose action will be computed numerically.
This is implemented by a discretization of the momentum fraction support.
The choice here will be an equally spaced $N_p\times N_p$ grid, which contains points
\begin{align}
 \bigl(x_{1,2}\bigr)_i &= \frac{2(i+1)}{N_p-1} + \frac{N_p+1}{1-N_p}\,,
\end{align}
where $ i\in\{0,\dots,N_p-1\}$.
Note that the maximum  and minimum values for $i$ correspond to the boundaries $\pm1$ of the support of the correlation functions.
Then, one can evaluate the action of $\mathbb H$ on a function $f$ (in very schematic notation)
\begin{align}
 \mathbb{H}f\bigl(x_i,x_j\bigr) &= \sum_k  w_k\mathbb{K}\bigl(x_i,x_j,x_k\bigr) f\bigl(x_k,x_i+x_j-x_k\bigr)\,,\label{eq:generich}
\end{align}
where $\mathbb{K}$ is the associated integral kernel of $\mathbb{H}$ and $w_k$ are certain weights for the sampling points of the integration, e.g.\ according to Simpson's rule or other Newton-Cotes formulas~\cite{abramowitzstegun}.
The structure appearing in Eq.~\eqref{eq:generich} follows from momentum conservation and the actual support of the integration/summation is determined by the support of $f$ as well as by the form of $\mathbb K$.
We refer to the appendix of~\cite{Braun:2009mi} for all possible kernels.

There is one subtlety specific to this type of evolution equations.
Certain values of momentum fractions, like the zero gluon momentum limit, usually emerging from ``plus-like'' prescriptions, require some caution, see~\cite{Braun:2009mi,Schafer:2012ra,Ma:2012xn,Kang:2012em}.
A possibility to deal with this would be to introduce some regulator, which avoids divisions by zero and ensures that all limits are taken correctly (up to errors of the order of the regulator).
The kernels implemented in this program treat this issue differently: they react on these ``special'' points appropriately with an analytic expression that needs to be hard-coded.

\section{Description of the program\label{sec:description}}

Usually the user wants to customize the program for his or her needs.
In order to do so, there are only a few places where one has to dig into the actual code, while one can widely ignore the low-level routines.
In the following an overview over the different parts of the program is presented and typical user-relevant modifications are pointed out.

\subsection{Source files}

The source of the program consists of $13$ files, which are listed below along with a short description of their contents.
\begin{verbatim}
 D1Function.cpp, D1Function.h
\end{verbatim}
contain a wrapper class for functions of one variable along with some useful routines (like calculating derivatives or saving to a file).
\begin{verbatim}
 HexFunction.cpp, HexFunction.h
\end{verbatim}
contain a wrapper class for functions of two variables with support on a hexagonal simplex. Examples are $\mathfrak{S}^\pm_q$, $\mathcal{F}^\pm$ and functions derived from it, see Sec.~\ref{sec:output}.
\begin{verbatim}
 qfqkernels.h, mixkernels.h,
 fffkernels.h, g2kernels.h
\end{verbatim}
implement the integral operators $\mathbb{H}^\pm_{QQ}$, $\mathbb{H}^\pm_{FQ,QF}$, $\mathbb{H}^\pm_{FF}$ and the integrals necessary for the extraction of $g_2$.
The initial conditions for the evolution can be given in
\begin{verbatim}
 initial.h,
\end{verbatim}
see Sec.~\ref{sec:initial} for details.
\begin{verbatim}
 constants.cpp, constants.h
\end{verbatim}
specify several global parameters for the evolution, see Sec.~\ref{sec:basicparameters}.
\begin{verbatim}
 mathutil.h
\end{verbatim}
contains the implementation of the numerical integration, viz. Eq.~\eqref{eq:generich}.
\begin{verbatim}
 t3evol.cpp
\end{verbatim}
represents the main program.

\subsection{How to compile}

In principle {\tt t3evol} can be compiled with any C++ compiler, for definiteness we assume it to be {\tt g++}.
The basic command to compile the program is
\begin{verbatim}
 g++ t3evol.cpp HexFunction.cpp
 D1Function.cpp constants.cpp
\end{verbatim}
which does not produce the fastest code.
The optimization options
\begin{verbatim}
 -O3 -fno-trapping-math
 -fomit-frame-pointer -funroll-loops
\end{verbatim}
provided by {\tt g++} can be used safely and produced an executable that did the calculation twice as fast compared to the unoptimized case.
A further performance enhancement (up to a few per cent) can be achieved by using profiles.
It is recommended to use all of the above.

In addition the unsafe option
\begin{verbatim}
 -ffast-math
\end{verbatim}
produced the fastest code in the test case.
The result of the computation did not change compared to the ``safe'' option scenario (although it could have).
There is no warranty that this may be true in a general situation, the use of this option is therefore not recommended.

\subsection{Basic parameters\label{sec:basicparameters}}

The files {\tt constants.h} and {\tt constants.cpp} contain some basic parameters for the simulation, which are listed and described below.
\begin{verbatim}
 Np
\end{verbatim}
is the number of discretization points.
Is is identical to $N_p$ from Sec.~\ref{sec:method}.
It should be an odd number, in order to have the central point $(x_1,x_2)=(0,0)$ and the boundary points $(\pm1,.)$, $(.,\pm1)$ lying on the grid.
The central point is characterized by an integer {\verb ZEROPOS }, which is equal to $(N_p-1)/2$ and should not be changed.
\begin{verbatim}
 NumSteps
\end{verbatim}
specifies the number of terms in the series of the formal solution in Eq.~\eqref{eq:solution}. For example, if {\verb NumSteps } is $4$ the exponential will be computed to the accuracy \mbox{$\mathcal{O}\bigl((t-t_0)^4\bigr)$}, neglecting terms of order $\mathcal{O}\bigl((t-t_0)^5\bigr)$.
\begin{verbatim}
 mN
\end{verbatim}
is the nucleon mass $m_N$ in units of $\giga\electronvolt$.
Its default value is set to $\unit{0.938}{\giga\electronvolt}$.
\begin{verbatim}
 Nc, oneoverNc, Cf
\end{verbatim}
are parameters of the special unitary group $\mathit{SU}(N_c)$, namely $N_c$, $1/N_c$ and $C_F=(N_c^2-1)/(2N_c)$, respectively.
Formally a large $N_c$ limit can be taken by setting {\verb oneoverNc } to zero and {\verb Cf } to $N_c/2$.
Of course, the default value is $N_c=3$.
\begin{verbatim}
 Nf
\end{verbatim}
is the number of different flavors participating in the evolution.
Since only up- and down-quark correlations are supported, it should be set to $2.0$.
It is treated as an independent parameter and does not interfere with the number of flavors in $\alpha_s(\mu^2)$, see below.

For the running coupling $\alpha_s(\mu^2)$ a couple of implementations at leading order and next-to-leading order are available.
The default treatment is analogous to GRV98LO, cf.~\cite{Gluck:1998xa} for details and values of the quantities below.
At lowest order the beta function is given in
\begin{verbatim}
 double beta0(double nf),
\end{verbatim}
which depends on the number of active flavors at the scale $\mu$, which in turn is given by
\begin{verbatim}
 int nf(double mu),
\end{verbatim}
which changes its value at the flavor thresholds
\begin{verbatim}
 FTH34, FTH45, FTH56.
\end{verbatim}
The matching parameters are adjusted accordingly in
\begin{verbatim}
 double lambda_LO_GRV98(int nf).
\end{verbatim}
From these parameters the running coupling and the evolution time is computed in
\begin{verbatim}
 double alpha_s_LO_GRV98(double mu)
 double evolutionTime_LO_GRV98(double
 mu0, double mu1, int nf),
\end{verbatim}
\noindent
viz. Eqs.~\eqref{eq:runningalphas}, \eqref{eq:evoltime}.
Other implementations can also be used by replacing each call of the latter function in the main program.

The rest of the constants are used for the numerical integration, cf.\ $w_k$ in Eq.~\eqref{eq:generich}, and can be ignored by the user.

\subsection{Initial distributions\label{sec:initial}}

There are two ways to specify the initial conditions, $\mathfrak{S}^\pm_0 $, $\mathcal{F}^\pm_0$.
The first way is to implement the functions inside the file {\tt initial.h}, namely
\begin{verbatim}
 initialFunction_Splus_u
 initialFunction_Splus_d
 initialFunction_Sminus_u
 initialFunction_Sminus_d
 initialFunction_Fplus
 initialFunction_Fminus
\end{verbatim}
for  $\mathfrak{S}^+_u(t_0)$, $\mathfrak{S}^+_d(t_0)$, $\mathfrak{S}^-_u(t_0)$, $\mathfrak{S}^-_d(t_0)$, $\mathcal{F}^+(t_0)$, $\mathcal{F}^-(t_0)$ respectively.
All of them are real functions of $x_1,x_2,x_3$.
The prescription of Eq.~\eqref{eq:fx1x2} is done by function overloading inside the same file.
By default the program is shipped with the model of~\cite{Braun:2011aw} and $\mathcal{F}^\pm(t_0) = 0$.

The second way is to read data from a file.
The program will check whether inside its directory one (or more) of the following files exist:
\begin{verbatim}
 S+_singlet_initial.txt
 S+_nonsinglet_initial.txt
 S-_singlet_initial.txt
 S-_nonsinglet_initial.txt
 F+_initial.txt
 F-_initial.txt
\end{verbatim}
They correspond to  $\mathfrak{S}^+_u(t_0)+\mathfrak{S}^+_d(t_0)$, $\mathfrak{S}^+_u(t_0)-\mathfrak{S}^+_d(t_0)$, $\mathfrak{S}^-_u(t_0)+\mathfrak{S}^-_d(t_0)$, $\mathfrak{S}^-_u(t_0)-\mathfrak{S}^-_d(t_0)$, $\mathcal{F}^+(t_0)$, $\mathcal{F}^-(t_0)$.
The scope of these input files should be such that it contains numbers in $N_p$ rows and $N_p$ columns, where the number at position $(i,j)$ corresponds to the value of the function at $\bigl((x_1)_i,(x_2)_j\bigr)$.
In fact, the program reads $N_p^2$ tokens, and any data beyond these will be ignored.

Note that the initialization files are automatically used if they are present.
If one does not want to use them, they need to be removed from the directory.
Therefore it is also possible to use a ``hybrid'' mode, giving some functions by file and the others by code.

At this point it is important to note that $\mathfrak{S}^\pm$ and $\mathcal{F}^\pm$ have to obey certain symmetry relations and restrictions on the support, cf.~\cite{Braun:2009mi}.
The user himself is responsible to ensure that these properties are fulfilled in his input.
Apart from that there is no restriction on the functional form for the initial conditions.
In that sense the evolution can be tested on any model, which seems to be necessary given the fact that the available models and estimates differ drastically~\cite{Braun:2011aw,Kanazawa:2010au,Boer:2011fx,Kang:2012xf,Metz:2012ui}.

\subsection{Running the program}

{\tt t3evol} can effectively be called with zero, one or two command line parameters.
If two parameters are given, it will assume that the first one is the initial scale $\mu_0^2$ (in $\giga\electronvolt^2$) and the second one is the final scale $\mu^2$ (in $\giga\electronvolt^2$).
If only one parameter is received, it is assumed that it corresponds to $\mu_0^2$ and the program awaits an input from the user for $\mu^2$.
When no parameter is specified, the program will ask for both $\mu_0^2$ and $\mu^2$.
Note that it is required that $\mu^2>\mu_0^2$ and obviously $\mu_0^2>0$.

The program then continues to do the evolution using the methods of Sec.~\ref{sec:method}, making a separate calculation for each flavor and each $C$-parity sector.
Solutions are written as soon as they are available.

\section{Output\label{sec:output}}

The results of the evolution will be written to several files in the directory of the program, each function to a separate file.
For convenience, this is done for a couple of equivalent representations owing to the variety of conventions in the literature.
One may safely modify the main program to restrict the output to the quantities of the user's interest,
The output filenames are of the form
\begin{verbatim}
 prefix_final.txt,
\end{verbatim}
where {\verb prefix } stands for an identifier for one of the functions widely used in the literature.
The multitude of functions is listed in Tab.~\ref{tab:prefixes}.
\begin{table}
\begin{center}
\begin{tabular}{ccc}
\hline\hline
 {\verb prefix }         & function                             & definition  \\
 \hline
  {\tt DeltaTdFd}        & $\Delta T_{\bar dFd}$                & \cite{Kang:2008ey,Braun:2009mi}\\
  {\tt DeltaTuFu}        & $\Delta T_{\bar uFu}$                & \cite{Kang:2008ey,Braun:2009mi}\\
  {\tt F-}               & $\mathcal{F}^-$                      & \cite{Braun:2009mi}\\
  {\tt F+}               & $\mathcal{F}^+$                      & \cite{Braun:2009mi}\\
  {\tt GFd}              & $G^d_{F}$                            & \cite{Kanazawa:2010au}\\
  {\tt GFu}              & $G^u_{F}$                            & \cite{Kanazawa:2010au}\\
  {\tt N}                & $N$                                  & \cite{Beppu:2010qn}\\
  {\tt O}                & $O$                                  & \cite{Beppu:2010qn}\\
  {\verb S-_d }          & $\mathfrak{S}^-_d$                   & \cite{Braun:2009mi}\\
  {\verb S+_d }          & $\mathfrak{S}^+_d$                   & \cite{Braun:2009mi}\\
  {\verb S-_nonsinglet } & $\mathfrak{S}^-_u -\mathfrak{S}^-_d$ & here\\
  {\verb S+_nonsinglet } & $\mathfrak{S}^+_u -\mathfrak{S}^+_d$ & here\\
  {\verb S-_singlet }    & $\mathfrak{S}^-_u +\mathfrak{S}^-_d$ & here\\
  {\verb S+_singlet }    & $\mathfrak{S}^+_u +\mathfrak{S}^+_d$ & here\\
  {\verb S-_u }          & $\mathfrak{S}^-_u$                   & \cite{Braun:2009mi}\\
  {\verb S+_u }          & $\mathfrak{S}^+_u$                   & \cite{Braun:2009mi}\\
  {\tt T3F- }            & $T_{3F}^-$                           & \cite{Braun:2009mi}\\
  {\tt T3F+ }            & $T_{3F}^+$                           & \cite{Braun:2009mi}\\
  {\tt TdFd}             & $T_{\bar dFd}$                       & \cite{Kang:2008ey,Braun:2009mi}\\
  {\tt TuFu}             & $T_{\bar uFu}$                       & \cite{Kang:2008ey,Braun:2009mi}\\
 \hline
\end{tabular}
\caption{\label{tab:prefixes}Prefix identifier for the generated output functions and reference(s) to their definition.}
\end{center}
\end{table}
These files have a structure of
\begin{align*}
 &x_1\quad x_2 \quad f(x_1,x_2)
\end{align*}
in each line, following the conventions of Eq.~\eqref{eq:fx1x2}.
In addition, for some of the antiquark-gluon-quark type correlations the so-called soft-gluon-pole and soft-fermion-pole configurations are exported in files of the form
\begin{verbatim}
 prefix_final_SGP.txt,
 prefix_final_SFP.txt.
\end{verbatim}
Their content is written in the form
\begin{align*}
 &x \quad f(x,0) \quad\text{(for {\tt SGP})}\,,\\
 &x \quad f(0,x) \quad\text{(for {\tt SFP})}
\end{align*}
in each line, again following the conventions of Eq.\ \eqref{eq:fx1x2}.

As a byproduct one can extract the twist-three contribution to the structure function $g_2$, see~\cite{Braun:2011aw}, via the relation
\begin{align}
 g_2^{\text{tw-3}}(x_B) &=  \sum_q \frac{e_q^2}{2m_N}\int \mathcal{D}x\,\mathfrak{S}_q^+(x)\notag\\
 &\quad \times\biggl(\frac{1 - P_{13}}{x_2} + \frac{d}{dx_3}\biggr)\frac{\theta(x_3 - x_B)}{x_2x_3}\,.
\end{align}
The resulting functions are written to
\begin{verbatim}
 xg2_tw3_proton_final.txt,
 xg2_tw3_neutron_final.txt,
\end{verbatim}
where $g_2$ for the neutron is obtained by a simple isospin rotation from the proton correlators.

\section{Conclusions\label{sec:conclusion}}

To the best of the author's knowledge, {\tt t3evol} represents the first closed solution to the twist-three renormalization group for the single spin asymmetry inducing multi-parton correlators.
It is a first step towards a correct incorporation of QCD evolution effects in phenomenological applications like a global fitting procedure.

The program, especially the implementation of the Hamiltonians, has been thoroughly tested for simple input functions, for which one can calculate the action of $\mathbb H$ analytically.
An example output for a realistic ansatz has already been presented in a model study based on light-cone wave functions~\cite{Braun:2011aw}.
In the particular case of $g_2$, a good agreement with effective approaches~\cite{Braun:2001qx} has been found.
The reproducibility of the anomalous dimensions of $g_2$ is a very strong check and indicates that the program is working correctly.
Further checks may be in order to clarify this claim.

The advantage to choose an arbitrary set of initial distributions is probably the biggest limitation on the performance.
If there was a more or less universal functional form for the correlators, it may be possible to invent more efficient approaches to deal with the evolution in terms of moments analogous to~\cite{Vogt:2004ns}.
Unfortunately the generalization of the Mellin-techniques to the problem at hand does not seem to be straightforward and requires further investigation.

An obvious way to increase the performance of the code in a multi processor environment is to use multithreading.
Since there are essentially four sectors (flavor (non-)singlet, $C$-parity $\pm$) that evolve independently, a parallelization of these calculations suggests itself.
This feature will be reserved for future versions.

\section*{Acknowledgements}

The author is grateful to V.M.~Braun and A.N.~Manashov for collaboration and to M.~Gruber for valuable discussions.
This project is partially supported by the DFG, grant BR2021/5-2 ``Multiparton evolution equations in QCD''.


\section{Introduction}
\label{}

Some games are complex systems which provide information about
social and biological processes \cite{Prost12,Blasius09PRL,Ribeiro13PO,Sigman10FN,Chassy11PO,
Sire09EPJB,Guerra12PA,Petersen08EPL,Bittner09EPJB,BenNaim07EPL,Heuer09EPJB,Ribeiro12PRE}.
They are currently being studied because of the existence of very well documented
registers of games which are useful for statistical analysis.
In particular, recent works have tried to understand the statistics of wins and losses in
baseball teams \cite{Sire09EPJB}, the final standing in basketball
leagues \cite{Guerra12PA}, the distribution of career longevity in baseball \cite{Petersen08EPL},
the football goal distribution \cite{Bittner09EPJB}, and
face to face game rank distributions \cite{BenNaim07EPL}.
The time evolution of the table during a season, for instance,  can be interpreted as a
random walk \cite{Heuer09EPJB}, and long-range correlations have been found in the score
evolution of the game of cricket \cite{Ribeiro12PRE}.
Among them, the  game of chess has a main place in occidental culture
because its intrinsic complexity is viewed as a symbol of intellectual prowess.
Since the skill level of chess players can be correctly identified \cite{Glickman95ACJ}
chess has  contributed to the scientific understanding  of expertise \cite{Chassy11PO}.
In addition, nowadays, there is a big world-wide community of chess players which
makes the game a benchmark for studying, for instance, decision making
processes \cite{Sigman10FN} and population level learning \cite{Ribeiro13PO}.

Exploring a chess database Blasius and T\"{o}njes \cite{Blasius09PRL} observed that the pooled
distribution of all opening weights follows a Zipf law with universal exponent.
This is a remarkable result \cite{Maslov09P,nature09} since the Zipf law is followed in
a range that  comprises  six orders of magnitude. In their work, the authors explain the
results using a multiplicative process with a branching ratio distribution which resembles
the real branching distribution measured in the database. Beyond the very good
agreement between their empirical observations and the proposed model to explain Zipf's law,
they do not provide an explanation of the evolution of the pool of games; disregarding
for instance, the question of possible memory effects between different games.
The development of a given chess game should depend on the expertise
of the players because the first moves (game opening) are usually known in advance
by high level players as a result of their theoretical background.
Other aspects can influence the development of a chess game according to the level of
expertise.
It has been established that skillful players develop outstanding rapid object
recognition abilities that differentiate them from the non-expert
players \cite{Gobet96MC, Gobet96CP}.
In particular, Gobet and Simon \cite{Gobet96CP} suggested that chess players, like experts
in other fields, use long-term memory retrieval structures  or templates in addition to
chunks \cite{Miller56PR,Chase73CP} in short-term memory in order to store information
rapidly.
Then, it can be expected that the differences in the level of expertise of chess
players are reflected in the historical development of chess databases introducing
correlations between games through memory effects.

Currently there is a big corpus of digitized texts which allows
the study of cultural trends tracing memory effects \cite{Diuk12FIN}, opening a new branch
in science known as culturomics \cite{Michel11S}.
In particular, long-range correlations have  been observed in literary corpora \cite{Montemurro02F,Altmann12PNAS}
where Zipf's law  was  also studied using extensive databases \cite{Petersen12SR}.
Testing signatures of memory effects is important because it has been shown that systems
which exhibit Zipf's law need a certain degree of coherence \cite{Cristelli12SR} for its
emergence.
Therefore, the detection of long-range memory effects can be useful to shed
light on the general mechanism behind the Zipf law.

In this work we explored the existence of long-range correlations in game sequences.
To that end, time-series are constructed using  a chronological ordered chess database
similar to the one used by Blasius and T\"{o}njes. In order to support the reliability of
our results, we used two mapping rules to build the discrete time series.
We also used different techniques to detect long-range correlations,namely the rescaled range
analysis and the detrended fluctuation analysis \cite{Beran94}.

\section{Long-range correlation analysis}

Rescaled range analysis ($R/S$) and detrended fluctuation analysis (DFA) are two useful
tools for detecting long-range correlations in discrete time series. The $R/S$ analysis
was introduced by Hurst when studying the regularization problem of the Nile
River \cite{Beran94}, whereas DFA was introduced  by Peng et al. \cite{Peng94PRE} for
detecting long-range (auto-) correlations in time series with non-stationarities.
Both methods use the accumulated time series which can be thought
as the displacement of a one dimensional random walker whose steps
are dictated by the values of the temporal series.

As a general procedure, the discrete time series $X(t)$ is first centered and normalized,

\begin{equation}
\label{eq1}
\tilde{X}(t)= \frac{X(t)-\langle X(t) \rangle}{\hat{\sigma}}
\end{equation}

\noindent where $\hat{\sigma}=\sqrt{\langle X(t)^2 \rangle - (\langle X(t) \rangle})^2$ and
$\langle ... \rangle$ means arithmetic averages over the complete series.
Normalization is useful in order to compare different assignation rules \cite{Montemurro02F}.
Then the accumulated series $Y(t)$ is constructed,

\begin{equation}
\label{eq2}
Y(t)= \sum_{u=1}^{t}\tilde{X}(u).
\end{equation}

Once the series $Y(t)$ is obtained, we applied the rescaled range
analysis \cite{Montemurro02F,Beran94} and the detrended fluctuation analysis  \cite{Peng92N,Kantelhardt01PA,Bryce12SR}.
In the $R/S$ method the time series $\tilde{X}(t)$ is divided into non-overlapping
intervals of size $n$. The range of each interval is computed,

\begin{equation}
\label{eq3}
R(n)= max[Y(1),Y(2), ... ,Y(n)] - min[Y(1),Y(2), ... ,Y(n)],
\end{equation}

\noindent and the corresponding standard deviation,

\begin{equation}
S(n)=\sqrt{\frac{1}{n}\sum_{i=1}^{n} {(X_{i}-m)}^{2}},
\end{equation}

\noindent where $m$ is the average of $X(t)$ in the interval of size $n$.
Then the average of the rescaled  range $E[R(n)/S(n)]$ over all the intervals of size $n$ is calculated.
The Hurst exponent $H$ is obtained varying the size of the
intervals and fitting the data to the expression $Cn^H$ where $H$ is the Hurst exponent.

In the DFA method, the $Y(t)$ series is also segmented into intervals of
size $n$. At each segment the cumulated series is fitted to a polynomial
$y_n(t)$ and the fluctuation $F(n)$ is obtained:

\begin{equation}
\label{eq5}
F(n)=\sqrt{\frac{1}{N}\sum_{i=1}^{N} {[y(i)-y_n(i)]}^{2}},
\end{equation}
\noindent where $N$ is the total number of data points. A log-log plot
of F(n) is expected to be linear and if the slope is less than unity it corresponds
to the Hurst exponent. When $H = 0.5$ the time series resembles a memoryless random
walker, and for $H > 0.5$ ($H < 0.5$) a random walker with persistent (antipersistent)
long-range correlations.

\section{Results}

In our work we use the ChessDB database\footnote{http://chessdb.sourceforge.net/.}
which has $3.5\times 10^6$ registered games; however, all our analysis is based on
a subset of $1.4\times10^6$  obtained by filtering corrupted data. The total time
spanned in the database is about $10$ years, from $1998$ to $2007$.
The first step in our analysis consists in mapping the chronologically ordered
chess games into a discrete time series.
The database contains approximately $380$ games per day uniformly distributed
over $\approx 3650$ days.
The games played on the same day are not chronologically ordered,
hence the minimum time resolution in  the studied temporal series is a day.
The mapping for building the time series is far from being unique, but as a
general rule it should  not introduce spurious
correlations \cite{Montemurro02F, Altmann12PNAS}.
We used two assignation rules.
In the first one, given a subset of consecutive
games we evaluate the similitude of the last of them with the rest.
We call this rule the similitude
assignation rule (SAR). We chose to compare the recent $t$--th game
against the rest using the following mapping expression,

\begin{equation}
\label{eq6}
X(t)=\sum_{t'=t-\tau}^{t-1}S_d(t,t'),
\end{equation}

\noindent where $S_d(t,t') = k \in \{0,1,2,3,...,d\}$ if and only if the $t$--th and $t'$--th games
 are identical up to move $k$ included. Here $\tau +1$  is the total number
of games in the set, and $d$ is the maximum number of moves to be
considered. In this work we considered full length games, i.e. $d$ is the length
of the shortest game to be compared. In our database the mean value of a game length
is $\approx 60$ and the longest game corresponds to $400$ moves.
According to Eq. (\ref{eq6}) when $\tau=1$ then $X(t)=S_d(t,t-1)$, meaning that each point of
the time series uses two consecutive  games. For instance, if the sequence of moves of
two games at $t$ and $t-1$ are --algebraic notation--
$e4,e5;Nf3,Nc6;Bb5,a6;...$ and $e4,e5;Nf3,Nc6;d3,d6;...$, respectively, then $X(t)=4$ because
these two games have four common moves.
In the second rule, $X(t)$ is the popularity at depth $d$ of the t-th game, where the popularity
at depth $d$ of the t-th game is equal to the number
of games in the database that have the same moves up to depth $d$.
We call this rule the popularity assignation rule (PAR) and is similar to the one
introduced by Montemurro and Pury \cite{Montemurro02F} for the study of a literary corpus.
The popularity of a given game is computed using its firsts $d$ movements.
Games that evolve in the same way
up to a depth $d$ determine a subset of games with the same popularity.
The popularity of a given game is equal to the number of elements of the subset to which
it belongs.
Notice that this assignation rule is based on a global characterization of the database;
however, it does not introduce spurious correlations as we will see.

The mapping rules we have chosen are roughly equivalent since games with a popular
sequence of movess are likely to match with other games. Popular games are related
to opening lines which develop to different depths. From the chess opening system
classification it is known that the game lines branch at different depths depending
on the opening system, for instance, Caro-Kann opening has a long sequence of
moves before branching.
However, in a global characerization \cite{Blasius09PRL} the tree of chess moves
has been proved to be self-similar to a certain degree, which means that the topology
 of the tree is nearly independent of the depth,
supporting the equivalence of both mapping rules.

\begin{figure
\includegraphics*[width=12cm,angle=0]{./figure1.eps}
\caption{\label{fig1} The top panels show fragments of the temporal series $X(t)$
(2000 points) obtained using  the assignation by similitude (a) and popularity (b).
The bottom panels show the complete integrated series, $Y(t)$,
similitude (c) and popularity (d).  Insets: distribution of similarities (a)
and popularities (b).  The straight lines correspond to fits using: $f(x)= A e^{\alpha x}$
and $g(x)=B x^{\beta}$, insets (a) and (b) respectively, where $\alpha =-0.49$ and $\beta=-1.53$.}

\end{figure}

The upper panels of Fig. \ref{fig1} show  fragments of the time series obtained using
the two assignation rules described above. On the right, we show the PAR using a game depth
$d=4$ and on the left  the  SAR using Eq.\ref{eq6} with $\tau=1$.
These parameters are appropriated for comparing the two assignation rules.
As we will see the mean value of the SAR series (with $\tau=1$) is small, meaning
that the SAR is testing the first moves, therefore $d=4$ for the PAR is convenient.
On the other hand, we have checked that the Hurst exponent is nearly independent  of $d$
up to $d=20$. For $d>20$ correlations are affected maybe because
the exponent of the power law which controls the distribution of popularities is much greater
than 2, therefore fluctuations in the popularities become relevant \cite{Blasius09PRL}.
In the SAR series the Hurst exponent is nearly independent of $\tau$ for  $\tau \leq 500$.
Beyond that value $H$ drops, and the SAR series do not show long-range correlations.
Note that in the SAR case, increasing $\tau$ implies a sort of time averaging.
The points of the series obtained by the PAR are distributed according to a scale free
distribution as expected.
In the corresponding inset we show this distribution, the exponent obtained by
fitting is in agreement with that obtained by Blasius and T\"{o}njes \cite{Blasius09PRL} (see caption).
SAR series have small values since most of the matches between games are restricted
to the first moves, in fact the  distribution of coincidences is exponential
and has a rapid decay (see inset and caption).
This means that most of the consecutive games are similar in the very first moves,
and a few of them have coincidences beyond the $10^{th}$ moves.
In the bottom panels we show the complete accumulated series for both cases, SAR on
the left and PAR on the right. Persistence can be appreciated from the monotonic
behavior of the curve over long times.

In Fig. \ref{fig2}, we show data corresponding to the $R/S$ analysis
(upper panel)
and  DFA (medium and lower panels). The left side panels show the analysis of the
SAR series  and the right one, that corresponding to the PAR series.
In all cases a log-log plot of the points can be well fitted by a straight line
as can be appreciated. The Hurst exponent obtained from both analyses $R/S$ and DFA
are similar; however, the exponent tends to be larger in the PAR series, $H \simeq 0.75$
as compared to the SAR series $\simeq 0.67$.
According to these values both series exhibit long-range correlations.
Since SAR is a local assignation rule and PAR is based on a global charactarization of the database, SAR series should be noisier than PAR series. According to this the Hurst exponent of SAR is expected to be less than that of the PAR series.
The fact that both methods give nearly the  same value for the
Hurst exponent indicates that  non-stationarities are  negligible in the time series.
When the SAR and PAR time series are randomly shuffled (lower panels), the Hurst exponent
in both analyses becomes close to $0.5$, which corresponds to series without long-range
correlations. This means that the assignation rules we used do not introduce
artificial correlations.

\begin{figure
\includegraphics*[width=12cm,angle=0]{./figure2.eps}
\caption{\label{fig2} $R/S$ analysis (top panel) and DFA analysis
(medium and lower panels) for the similitude (left) and popularity
(right)  assignation rules. The straight dashed lines are linear fits of the data.
The $R/S$ error points  (black crosses) correspond to the maximum and minimum point of
the averaged data.}
\end{figure}

In order to compare the correlations obtained from different subsets  of the total
time series we will first study the dependence of the Hurst exponent on the length
of the series.
In the upper panel of Fig. \ref{fig3} we plotted the  Hurst exponent as
function of the length $N$ of the PAR (full line) and SAR (dashed line) series.
In all this analysis we used the DFA method and we averaged over all
possible sets of consecutive games of length $N$.
A certain correlation can be appreciated between the Hurst exponent of the two series
assignations. The exponent $H$ of the PAR assignation is always greater than that of the
SAR assignation, as in Fig.\ref{fig2}. In both cases $H$ grows as the
length of the series increases reaching a maximum value at nearly $2\times10^5$
and then nearly stabilizes to $H\simeq 0.75$ in the PAR series and $H \simeq 0.65$ in the SAR series.
In the case of the shuffled series (medium panel) we do not observe size effects and
$H$ fluctuates around $0.5$ as expected. It is worth noting that in this analysis, a change
in the  size of the time series is equivalent to a change in the total time spanned.
Finally, in the lower panel of Fig.\ref{fig3} we show the effects of
reducing the series by pruning them. The error bars are obtained from the fitting errors.
In this case we constructed reduced series by  removing  points from complete series but
keeping nearly constant the total time span.
More specifically, to construct the set of reduced series we pick up one point from the
complete series every $s$ points, where $s = 2, 4, 8, ...$.
The result after pruning is similar to that observed in the upper panel.
In both cases, the exponent grows as the size of the series increases.
However, in the pruned series $H$ seems not to be stabilized even in the larger series.
\begin{figure
\includegraphics*[width=12cm,angle=0]{./figure3.eps}
\caption{\label{fig3} Upper panel: Hurst exponent obtained by the DFA method as
function of the length of the time series using the two assignation rules, PAR and
SAR series, full and dashed lines, respectively. Meddle panel: Hurst exponent obtained
from the shuffled series. Lower panel: $H$ versus size in the pruned series.
}
\end{figure}
Since the level of the players can be ranked according to their ELOs,
a rating system introduced by the physicist Arpad Elo \cite{Glickman95ACJ},
we repeated the analysis filtering the database according to different ELO intervals.
We used the following ELO's ranges \cite{Chassy11PO}: $[1,2199]$, $[2200, 2399]$ and above
$2400$. These intervals nearly correspond to: non-experts and master candidates, masters,
and grand masters, respectively. In particular, this partition
also ensures each interval contains nearly the same number of games.
In Fig. \ref{fig4}, we show the Hurst exponent as function of the length of the time series
associated to the ELO's interval defined above.
In the upper panel, we show the distribution of ELO's rank of our database.
As expected it is well described by a  Gaussian function $\exp(\frac{-(x-x_0)^2}{2a^2})$ (continuous line),
with the fitted values, $x_0=2303$ and $a=203$.
On the left we show the  SAR series; and on the right, the PAR series. As in Fig. \ref{fig3},
the behavior of the two assignation
rules is similar and the exponents of PAR series are  slightly greater than in the SAR series.
At short times the Hurst exponent in the three ELO intervals is close to 0.5,
i.e. at this time scale games are uncorrelated and no persistence can be detected.
At intermediate times, the lower ELO intervals show stronger correlations;
and at the longest times, the order reverses and the games of high level players are
the most correlated.

\begin{figure
\includegraphics*[width=12cm,angle=0]{./figure4.eps}
\caption{\label{fig4} Hurst exponents obtained by the DFA method computed in different
ELO intervals as function of the length of the series. SAR left lower
panel and PAR rule right lower panel. Upper panel:  ELO distribution
of the whole database, the vertical lines indicates the ELO's intervals.
The continuous line correspond to a fit using a Gaussian function,
$f(x)=\exp(\frac{-(x-x_0)^2}{2a^2})$, with $x_0=2303$ and $a=203$.}
\end{figure}

When a pruning is implemented in series filtered according to ELO ranks,
the Hurst exponent decreases as the length of the series decreases, as in the
complete database (see Fig.\ref{fig5}). The error bars are obtained from the fitting errors.
At variance with the results of Fig.\ref{fig4}
the series of higher ELOs are always more correlated, suggesting a different
origin of this behavior.

\begin{figure
\includegraphics*[width=12cm,angle=0]{./figure5.eps}
\caption{\label{fig5} Size effects in the Hurst exponent for pruned series in
different ELO ranges. Upper panel PAR series and lower panel SAR series.
}
\end{figure}

\section{Summary and Discussion}
\label{conclu}

The Hurst exponents obtained through $R/S$ and DFA analysis are similar
indicating that the time series analyzed are stationary for both assignation rules,
SAR and PAR.
For the case of PAR and the $R/S$ analysis our results are similar to that reported
in literary corpora \cite{Montemurro02F}.
The value $H\simeq0.5$ that emerges from the shuffled series implies that
the assignation rules we used to build the series do not introduce spurious correlations.
The analysis of the series as a function of the time spanned shows that there
is a threshold in the length of the time series from which long-range correlations
appear. This threshold depends on the assignation rule and on the level of
the players.
When the database is filtered by ELO ranks, our results indicate that long-range
correlations observed at large time scales are related to the presence of high level
players. The games corresponding to intermediate and low level players show correlations
at shorter times  when compared to high level players.
A reduction of the database by pruning indicates that in short series correlations disappear.
At variance with the analysis based on the length of the time spanned, the reduction in
correlations by pruning has a similar behavior irrespective of the ELO ranks analyzed.
Therefore, the reduction  of correlations in this case should have a different origin.

Size effects in the detection of long-range correlations have been analyzed in the work
of Coronado {\it et al.} \cite{Coronado05JBP} using time series generated by the Makse
algorithm \cite{Makse96PRE}. They observed that when DFA is used, size effects are
negligible even in series containing as few as $10^3$ points.
Since we use DFA and in our series size effects were observed
in series having more than $10^3$ elements,
our results must be reflecting intrinsic properties of the system.
However, in our case we have portions of the time series which are randomly shuffled due
to the minimum time resolution of our database, and this is affecting the detection
of persistence in the time series.
In long-memory analysis of persistence if missing data are  replaced using random data drawn from
the empirical distribution of the time series, the Hurst exponent tends to decrease \cite{Wilson03PRE}.
Notice that pruning  does not affect the total time  spanned in the series.
When pruning, the parameter $s$ can go up to $s \approx 380$ and still have on average
one game per day in the pruned series with a time series of length $\approx 3,650$.
According to this, the pruned time series we used should be representative samples of
the temporal evolution of the games.
However, as was pointed out, in pruned series correlations decrease.
In our case, a possible explanation is that pruning affects the representativity of samples because
persistence is a collective effect.
The samples in order to be representative have to include at each time the behavior of a
representative pool of different players.
When series are pruned, the number of samples taken at the minimum resolution time (one day) are
reduced hence the behavior of the pool of players is not correctly probed.

For intermediate and low level players persistence reaches its maximum in two or three years
and then nearly stabilizes. This behavior is independent of the origin of time used
for the analysis. In fact, the points shown in Fig. \ref{fig4} result from the averaging
of equivalent non-overlapping time intervals. For the most specialized players a signature
of persistence appears after one or two years. According to our results, it seems that
very high level players use different strategies as compared to low level players.
It seems they vary the lines of play more often at short times scales, maybe because outstanding players know many opening lines in depth and are less influenced
by the opponents.
In particular, in high level players the correlations at large time scale are the strongest.

This work was partially supported by grants from CONICET, SeCyT--Universidad
Nacional de C\'ordoba (Argentina), and the Academy of Finland, project no 260427.
The authors thank Pedro Pury for useful discussion.

\bibliographystyle{elsart-num}