Correlation Functions of Quantum Artin System

Abstract

It was conjectured by Maldacena, Shenker and Stanford that the classical chaos can be diagnosed in thermal quantum systems by using an out-of-time-order correlation function. The Artin dynamical system defined on the fundamental region of the modular group SL(2,Z) represents a well defined example of a highly chaotic dynamical system in its classical regime. We investigated the influence of the classical chaotic behaviour on the quantum–mechanical properties of the Artin system calculating the corresponding out-of-time-order thermal quantum–mechanical correlation functions. We demonstrated that the two- and four-point correlation functions of the Liouiville-like operators decay exponentially with temperature dependent exponents and that the square of the commutator of the Liouiville-like operators separated in time grows exponentially.

1. Introduction

The hyperbolic Anosov C-systems have exponential instability of their trajectories and as such represent the most natural chaotic dynamical systems [1]. Of special interest are C-systems that are defined on closed surfaces of the Lobachevsky plane of constant negative curvature1. An example of such a system has been introduced in a brilliant article published in 1924 by the mathematician Emil Artin [6]. The dynamical system is defined on the fundamental region of the Lobachevsky plane which is obtained by the identification of points congruent with respect to the modular group $SL(2,Z)$, a discrete subgroup of the Lobachevsky plane isometries [7,8,9]. The fundamental region $\mathcal{F}$ in this case is a hyperbolic triangle in Figure 1. The geodesic trajectories are bounded to propagate on the fundamental hyperbolic triangle. The geodesic flow in this fundamental region represents a highly chaotic dynamical system with exponential instability of its trajectories, has mixing of all orders, Lebesgue spectrum and non-zero Kolmogorov entropy [1,2,3,4,5,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28].

In the classical regime the correlation functions are defined as an integral over a pair of functions/observables A and B in which the first one is stationary and the second one evolves with the geodesic flow $g_t$:

$$\begin{array}{ccc}\hfill {\mathcal{D}}_t(A,B)& =& {\int }_{\mathcal{M}}A\left(g\right)\overline{B\left(g{g}_t\right)}d\mu \left(g\right).\hfill \end{array}$$

(1)

Earlier investigation of the classical correlation functions of the geodesic flows was performed in [17,18,19,20,29] by using different approaches including Fourier series for the $SL(2,R)$ group, zeta functions for the geodesic flows, relating the poles of the Fourier transform of the correlation functions to the spectrum of an associated Ruelle operator, the methods of unitary representation theory, spectral properties of the corresponding Laplacian and other approaches. In recent articles [30,31] the authors demonstrated exponential decay of the correlation functions with time on the classical phase space. The result was derived using differential geometry, the group-theoretical methods of Gelfand and Fomin, the time evolution equations, and the properties of automorphic functions on $\mathcal{F}$. The exponential decay rate was expressed in terms of the entropy $h\left(\mathcal{F}\right)$ of the system:

$$\begin{array}{ccc}\hfill |{\mathcal{D}}_t(A,B)|& \le & M{e}^{-Kh\left|t\right|},\hfill \end{array}$$

(2)

where M and K are constants depending on the smoothness of the functions. In the classical regime the exponential divergence of the geodesic trajectories results in the universal exponential decay of its classical correlation functions [30,31].

In the present article, we are interested to study the quantum–mechanical correlation functions of highly chaotic Artin system. The interest in studying the quantum–mechanical behaviour of highly chaotic systems is rooted in the fact that the fundamental interactions, electroweek and strong interactions are described by the Yang–Mills gauge theory which is highly chaotic in its classical regime [21,22,27,32,33,34,35]. In recent years the quantum–mechanical properties of highly chaotic systems were widely discussed in series of publications [26,27,33,34,36]. The investigation was initiated by the analysis of black holes evaporation thermodynamics [26,27,36]. Therefore, there is a great interest in the quantisation of the highly chaotic dynamical systems and investigation of their quantum–mechanical properties. It is important to learn how and to what extent the classical chaos influences the quantum–mechanical behaviour. Our main focus in this article is to study the behaviour of the quantum–mechanical correlation functions of the Artin system because on the one hand it is highly chaotic and on the other hand it has a deep mathematical origin. It was conjectured by Maldacena, Shenker and Stanford [26] that the classical chaos influences the thermal properties of the corresponding quantum–mechanical system and can be traced by investigation of out-of-time-order correlation functions, as well as of the square of the commutator of the operators separated in time.

In order to investigate the behaviour of the correlation functions in the quantum–mechanical regime it is necessary to know the spectrum of the system and the corresponding wave functions. In the case of the modular group the energy spectrum has a continuous part, which originates from the asymptotically free motion inside an infinitely long vertical “y-channel” extended in the vertical direction of the fundamental region, as well as infinitely many discrete energy states corresponding to bounded motion at the “bottom” of the fundamental triangle. The spectral problem has a deep number-theoretical origin and was partially solved in a series of pioneering articles [37,38,39,40]. It was only partially solved because the discrete spectrum and the corresponding wave functions are not known analytically. The general properties of the discrete spectrum have been derived using the Selberg trace formula [39,40,41,42,43,44]. Numerical calculations of the discrete energy levels were performed for many energy states [45,46,47].

In the next section we shall describe the geometry of Lobachevsky hyperbolic plane and of the fundamental region which corresponds to the modular group $SL(2,Z)$, the geodesic flow on that region and the quantisation of the system. The derivation of the Maass wave functions [37] for the continuous spectrum will be reviewed in detail. We shall use the Poincaré representation for the Maass non-holomorphic automorphic wave functions. We introduce a natural physical variable $\tilde{y}$ for the distance in the vertical direction2 on the fundamental triangle ${\int }_{dx=0}ds=\int dy/y=lny=\tilde{y}$ and the corresponding momentum $p_y$ in order to represent the Maass wave functions (54) in a physically insightful manner

$$\begin{array}{c}\hfill {\psi }_{p_y}(x,\tilde{y})=e^{-i{p}_y\tilde{y}}+\frac{\theta (\frac{1}{2}+i{p}_y)}{\theta (\frac{1}{2}-i{p}_y)}{e}^{i{p}_y\tilde{y}}+\frac{4}{\theta (\frac{1}{2}-i{p}_y)}\sum _{l=1}^{\infty }{\tau }_{i{p}_y}\left(l\right)K_{i{p}_y}\left(2\pi l{e}^{\tilde{y}}\right)cos\left(2\pi lx\right),\end{array}$$

(3)

where the functions $\theta \left(s\right)$, ${\tau }_{iu}\left(l\right)$ and $K_{iu}\left(y\right)$ will be defined precisely in Section 3 and are given by the Formulas (48), (44) and (41). Indeed, the first two terms describe the incoming and outgoing plane waves. The plane wave $e^{-i{p}_y\tilde{y}}$ incoming from infinity of the y axis in Figure 2 (the vertex $\mathcal{D}$) elastically scatters on the boundary $ACB$ of the fundamental triangle $\mathcal{F}$ shown in Figure 1. The reflection amplitude is a pure phase and is given by the expression in front of the outgoing plane wave $e^{i{p}_y\tilde{y}}$:

$$\frac{\theta (\frac{1}{2}+i{p}_y)}{\theta (\frac{1}{2}-i{p}_y)}=exp\left[i\phi \left(p_y\right)\right].$$

(4)

The rest of the wave function describes the standing waves $cos\left(2\pi lx\right)$ in the x direction between boundaries $x=\pm 1/2$ with the amplitudes $K_{i{p}_y}\left(2\pi ly\right)$, which are exponentially decreasing with index l Figure 2. The continuous energy spectrum is given by the formula

$$E=p_y^2+\frac{1}{4}.$$

(5)

The wave functions of the discrete spectrum have the following form [37,38,39,40,45,46,47]:

$$\begin{array}{ccc}\hfill {\psi }_n\left(z\right)& =& \sum _{l=1}^{\infty }{c}_l\left(n\right)K_{i{u}_n}\left(2\pi l{e}^{\tilde{y}}\right)\left\{\begin{array}{c}cos\left(2\pi lx\right)\hfill \\ sin\left(2\pi lx\right)\hfill \end{array},\right.\hfill \end{array}$$

(6)

where the spectrum $E_n=\frac{1}{4}+u_n^2$ and the coefficients $c_l\left(n\right)$ are not known analytically, but were computed numerically for many values of n [45,46,47].

Having in hand the explicit expression of the wave function one can analyze a quantum–mechanical behaviour of the correlation functions defined in [26,48]:

$$\begin{array}{cc}\hfill & {\mathcal{D}}_2(\beta ,t)=\langle A\left(t\right)B\left(0\right)e^{-\beta H}\rangle ,{\mathcal{D}}_4(\beta ,t)=\langle A\left(t\right)B\left(0\right)A\left(t\right)B\left(0\right)e^{-\beta H}\rangle \end{array}$$

(7)

$$\begin{array}{cc}\hfill & C(\beta ,t)=-\langle {[A\left(t\right),B\left(0\right)]}^2{e}^{-\beta H}\rangle ,\end{array}$$

(8)

where in our case the operators A and B are chosen to be of the Liouiville type:

$$A\left(N\right)=e^{-2N\tilde{y}},N=1,2,\dots \dots$$

(9)

Analyzing the basic matrix elements of the Liouiville-like operators (9) we shall demonstrate that all two- and four-point correlation functions (7) decay exponentially with time:

$${\mathcal{D}}_2(\beta ,t)\sim exp\left(-\frac{t}{t_2\left(\beta \right)}\right)\mathrm{and}{\mathcal{D}}_4(\beta ,t)\sim exp\left(-\frac{t}{t_4\left(\beta \right)}\right),$$

(10)

with exponents $t_2\left(\beta \right)$ and $t_4\left(\beta \right)$ that depend on temperature. These exponents define the dissipation time and reflect the fact that with increasing temperature the correlation functions will universally decay. Their behaviour is shown in Figure 3 and Figure 4.

Alternatively to the exponential decay of the correlation functions (7), the square of the commutator of the Liouiville-like operators separated in time (8) grows exponentially Figure 5. This growth is reminiscent to the local exponential divergence of trajectories in the Artin system when it is considered in the classical regime [30,31]. At short time scales the exponential growth of the commutator (8) shown in Figure 5 has the following form:

$$C(\beta ,t)\sim exp\left(\frac{2\pi }{\chi \left(\beta \right)}t\right),$$

(11)

where the behaviour of the exponent $\chi \left(\beta \right)$ is presented in Figure 5. The temperature dependence of the exponent $\chi \left(\beta \right)$ is smooth and there are no signs of any kind of phase transitions.

In our calculation of the quantum–mechanical correlation functions we shall use a perturbative expansion in which the high-mode Bessel function in (54) and (59) will be considered as perturbations. We found that our calculations are stable with respect to these perturbations and do not influence the final results. The reason is that in the integration region of the matrix elements (62) the high-mode Bessel functions are exponentially small.

2. Quantisation of Artin System

Let us start with Poincare model of the Lobachevsky plane, i.e., the upper half of the complex plane: H = $\{z\in \mathbb{C}$, $\Im z>0\}$ supplied with the metric (we set $z=x+iy$)

$$\begin{array}{c}\hfill d{l}^2=\frac{d{x}^2+d{y}^2}{y^2},\end{array}$$

(12)

which has Ricci scalar $R=-2$. Isometries of this space are given by $SL(2,\mathbb{R})$ transformations. The $SL(2,\mathbb{R})$ matrix (a, b, c, d are real and $ad-bc=1$)

$$g=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$$

acts on a point z by linear fractional substitutions:

$$\begin{array}{c}\hfill z\to \frac{az+b}{cz+d}.\end{array}$$

(13)

Note also that g and $-g$ give the same transformation, hence the effective group is $SL(2,\mathbb{R})/{\mathbb{Z}}_2$. We’ll be interested in the space of orbits of a discrete subgroup $G\subset SL(2,\mathbb{R})$ in H. Our main example will be the modular group $G=SL(2,\mathbb{Z})$. A nice choice of the fundamental region $\mathcal{F}$ of $SL(2,\mathbb{Z})$ is displayed in Figure 1. The fundamental region $\mathcal{F}$ of the modular group consists of those points between the lines $x=-\frac{1}{2}$ and $x=+\frac{1}{2}$ that lie outside the unit circle in Figure 1. The modular triangle $\mathcal{F}$ has two equal angles $\alpha =\beta =\frac{\pi }{3}$ and the third one at $\mathcal{D}$ equal to zero, $\gamma =0$, thus $\alpha +\beta +\gamma =2\pi /3<\pi$. The area of the fundamental region is finite and equal to $\frac{\pi }{3}$ and has the topology of sphere due to “gluing” the opposite edges of the triangle. The invariant area element on the Lobachevsky plane is proportional to the square root of the determinant of the metric (12): $d\mu \left(z\right)=\frac{dxdy}{y^2},$, thus $\mathrm{Area}\left(\mathcal{F}\right)={\int }_{-\frac{1}{2}}^{\frac{1}{2}}dx{\int }_{\sqrt{1-x^2}}^{\infty }\frac{dy}{y^2}=\frac{\pi }{3}.$

Consider geodesic flow on $\mathcal{F}$, which is conveniently described by the least action principle $\delta S=0$, where (cf. (12)):

$$\begin{array}{c}\hfill S=\int Ldt=\int \frac{\sqrt{{\dot{x}}^2+{\dot{y}}^2}}{y}dt.\end{array}$$

(14)

By varying the action, we get the equations of motion

$$\begin{array}{c}\frac{d}{dt}\frac{\dot{x}}{y\sqrt{{\dot{x}}^2+{\dot{y}}^2}}=0,\frac{d}{dt}\frac{\dot{y}}{y\sqrt{{\dot{x}}^2+{\dot{y}}^2}}+\frac{\sqrt{{\dot{x}}^2+{\dot{y}}^2}}{y^2}=0.\hfill \end{array}$$

(15)

Notice the invariance of the action and of the equations of motion under time reparametrizations $t\to t\left(\tau \right)$. Presence of a local (“gauge”) symmetry indicates that we have a constrained dynamical system. One particularly convenient choice of gauge fixing specifying the time parameter t proportional to the proper time, is archived by imposing the condition

$$\begin{array}{c}\hfill \frac{{\dot{x}}^2+{\dot{y}}^2}{y^2}=2H,\end{array}$$

(16)

where H is a constant. In this gauge the Equation (15) will take the following form [49]:

$$\begin{array}{c}\frac{d}{dt}\left(\frac{\dot{x}}{y^2}\right)=0,\frac{d}{dt}\left(\frac{\dot{y}}{y^2}\right)+\frac{2H}{y}=0.\hfill \end{array}$$

(17)

Defining the canonical momenta $p_x$, $p_y$ conjugate to the coordinates x, y as

$$\begin{array}{c}\hfill p_x=\frac{\dot{x}}{y^2},p_y=\frac{\dot{y}}{y^2},\end{array}$$

(18)

we shall get the geodesic Equation (17) in the Hamiltonian form:

$$\begin{array}{c}\hfill \dot{p_x}=0,\dot{p_y}=-\frac{2H}{y}.\end{array}$$

(19)

Indeed, after defining the Hamiltonian as

$$\begin{array}{c}\hfill H=\frac{1}{2}{y}^2(p_x^2+p_y^2)\end{array}$$

(20)

the corresponding equations will take the following form:

$$\begin{array}{c}\dot{x}=\frac{\partial H}{\partial p_x}=y^2{p}_x,\dot{y}=\frac{\partial H}{\partial p_y}=y^2{p}_y\hfill \\ \dot{p_x}=-\frac{\partial H}{\partial x}=0,\dot{p_y}=-\frac{\partial H}{\partial y}=-y(p_x^2+p_y^2)=-\frac{2H}{y}\hfill \end{array}$$

(21)

and coincide with (18) and (19). The advantage of the gauge (16) is that the Hamiltonian (20) coincides with the constraint. Now it is fairly standard to quantize this Hamiltonian system. We simply replace in (20)

$$p_x=-i\frac{\partial }{\partial x},p_y=-i\frac{\partial }{\partial y}$$

and consider the (time independent) Schrödinger equation

$$H\psi =E\psi .$$

The resulting equation explicitly reads:

$$\begin{array}{c}\hfill -y^2({\partial }_x^2+{\partial }_y^2)\psi =E\psi .\end{array}$$

(22)

On the left hand side one recognises the Laplace operator (with an extra minus sign) in Poincare metric (12). It is easy to see that the Hamiltonian is a positive semi-definite Hermitian operator. Indeed, for any quadratically integrable function $\psi (x,y)$ that fulfills the periodic boundary condition (26)

$$\begin{array}{c}\hfill -\int {\psi }^{*}(x,y)y^2({\partial }_x^2+{\partial }_y^2)\psi (x,y)\frac{dxdy}{y^2}=\int \left(\right|{\partial }_x\psi (x,y){|}^2+|{\partial }_y\psi (x,y){|}^2)dxdy\ge 0.\end{array}$$

(23)

It is convenient to introduce parametrization of the energy $E=s(1-s)$ and to rewrite this equation as

$$\begin{array}{c}\hfill -y^2({\partial }_x^2+{\partial }_y^2)\psi (x,y)=s(1-s)\psi (x,y).\end{array}$$

(24)

As far as E is real and semi-positive and parametrisation is symmetric with respect to $s\leftrightarrow 1-s$ it follows that the parameter s should be chosen within the range

$$s\in [1/2,1]\mathrm{or}s=1/2+iu,u\in [0,\infty ].$$

(25)

One should impose the “periodic” boundary condition on the wave function with respect to the modular group

$$\begin{array}{c}\hfill \psi \left(\frac{az+b}{cz+d}\right)=\psi \left(z\right),\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\in SL(2,Z)\end{array}$$

(26)

in order to have the wave function which is properly defined on the fundamental region $\overline{\mathcal{F}}$ shown in Figure 1. Taking into account that the transformation $T:z\to z+1$ belongs to $SL(2,Z)$, one has to impose the periodicity condition $\psi \left(z\right)=\psi (z+1)$. Thus we have a Fourier expansion

$$\begin{array}{c}\hfill \psi (x,y)=\sum _{n=-\infty }^{\infty }{f}_n\left(y\right)exp\left(2\pi inx\right).\end{array}$$

(27)

Inserting this into Equation (24), for the Fourier component $f_n\left(y\right)$ we get

$$\begin{array}{c}\hfill \frac{d^2{f}_n\left(y\right)}{d{y}^2}+(s(1-s)-4{\pi }^2{n}^2)f_n\left(y\right)=0.\end{array}$$

(28)

For the case $n\ne 0$ the solution which exponentially decays at large y reads

$$\begin{array}{c}\hfill f_n\left(y\right)=\sqrt{y}{K}_{s-\frac{1}{2}}\left(2\pi n\right|y\left|\right),\end{array}$$

(29)

and for $n=0$ one simply gets

$$\begin{array}{c}\hfill f_0\left(y\right)=c_0{y}^s+c_0^{{}^{\prime }}{y}^{1-s}.\end{array}$$

(30)

Thus the solution can be represented in the following form:

$$\begin{array}{c}\hfill \psi (x,y)=c_0{y}^s+c_0^{{}^{\prime }}{y}^{1-s}+\sqrt{y}\sum _{\genfrac{}{}{0pt}{}{n=-\infty }{n\ne 0}}^{\infty }{c}_n{K}_{s-\frac{1}{2}}\left(2\pi n\right|y\left|\right)exp\left(2\pi inx\right),\end{array}$$

(31)

where the coefficients $c_0,c_0^{{}^{\prime }},c_n$ should be defined so that the wave function fulfills the boundary conditions (26). Thus one should impose also the invariance with respect to the second generator of the modular group $SL(2,Z)$, that is, with respect to the transformation $S:z\to -1/z$:

$$\begin{array}{c}\hfill \psi \left(z\right)=\psi (-1/z).\end{array}$$

(32)

This functional equation defines the coefficients $c_0,c_0^{{}^{\prime }},c_n$. We found that it is much easier to resolve it by using the full group of $SL(2,Z)$ transformations acting on a particular solution (30) in order to derive the eigenfunctions obtained by Maass [37]. The wave function generated in this way will be invariant with respect to the $SL(2,Z)$ transformations. We shall follow this approach in the next section.

3. Continuous Spectrum and the Reflection Amplitude

As we just mentioned above in order to get $SL(2,Z)$ invariant solutions, one should define the coefficients $c_0,c_0^{{}^{\prime }}$ and $c_n$ in (31). Another option is to start from a particular solution and perform summation over all nonequivalent shifts of the argument by the elements of $SL(2,Z)$, that is, by using the Poincaré series representation [7,8]. We will demonstrate this strategy by using the simplest solution (30), (31) with $c_0=1,c_0^{{}^{\prime }}=0$:

$$\psi \left(z\right)=y^s={(\Im z)}^s.$$

Let us denote by ${\Gamma }_{\infty }$ the subgroup of $\Gamma =SL(2,Z)$, generating shifts $z\to z+n$, $n\in Z$. Explicitly the elements of ${\Gamma }_{\infty }$ are given by $2\times 2$ matrices:

$$\begin{array}{c}\hfill g_n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right).\end{array}$$

(33)

Since $y^s$ is already invariant with respect to ${\Gamma }_{\infty }$, we should perform summation over the conjugacy classes ${\Gamma }_{\infty }\setminus \Gamma$. Let us define these conjugacy classes. If two $SL(2,Z)$ matrices

$$\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\mathrm{and}\left(\begin{array}{cc}{a}^{\prime }& b^{\prime }\\ c^{\prime }& d^{\prime }\end{array}\right)$$

belong to the same class, then by definition for some $n\in Z$

$$\left(\begin{array}{cc}{a}^{\prime }& b^{\prime }\\ c^{\prime }& d^{\prime }\end{array}\right)=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$$

so that $c^{\prime }=c$, $d^{\prime }=d$, $a^{\prime }-a=nc$ and $b^{\prime }-b=nd$. Since $ad-bc=1$, it follows that a and c do not have a common devisor. In fact, the opposite is also true. Given a pair of mutually prime integers $(c,d)$ it is always possible to find a pair of integers $(a,b)$ such that $ad-bc=1$. For any other pair $(a^{\prime },b^{\prime })$ satisfying the same condition $a^{\prime }d-b^{\prime }c=1$, the relations $a^{\prime }-a=nc$ and $b^{\prime }-b=nd$ are satisfied for some integer n. Thus we established a bijection between the set of mutually prime pairs $(c,d)$ with $(c,d)\ne (0,0)$ and the set of conjugacy classes ${\Gamma }_{\infty }\setminus \Gamma$. The fact that the integers $(c,d)$ are mutually prime integers means that their greatest common divisor (gcd) is equal to one: $gcd(c,d)=1$. As a result, it is defined by the classical Poincaré series representation [7,8] and for the sum of our interest we get3

$$\begin{array}{c}\hfill {\psi }_s\left(z\right)\equiv \frac{1}{2}\sum _{\gamma \in {\Gamma }_{\infty }\setminus \Gamma }{(\Im \left(\gamma z\right))}^s=\frac{1}{2}\sum _{\genfrac{}{}{0pt}{}{(c,d)\in {\mathbb{Z}}^2}{gcd(c,d)=1}}\frac{y^s}{{({(cx+d)}^2+c^2{y}^2)}^s},\end{array}$$

(34)

where, as explained above, the sum on the r.h.s. is taken over all mutually prime pairs $(c,d)$. The series (34) is convergent when $\Re s>1$. We used also the simple relation

$$\Im \gamma z\equiv \Im \frac{az+b}{cz+d}=\frac{y}{{(cx+d)}^2+c^2{y}^2}.$$

To simplify further the sum let us multiply both sides of the Equation (34) by [37]

$$\sum _{n=1}^{\infty }\frac{1}{n^{2s}}\equiv \zeta \left(2s\right)$$

so that we shall get

$$\begin{array}{c}\hfill \zeta \left(2s\right){\psi }_s\left(z\right)=\frac{1}{2}\sum _{n=1}^{\infty }\sum _{\genfrac{}{}{0pt}{}{(c,d)\in {\mathbb{Z}}^2}{gcd(c,d)=1}}\frac{y^s}{{({(ncx+nd)}^2+{\left(nc\right)}^2{y}^2)}^s}.\end{array}$$

(35)

It is easy to be convinced now that the set of all pairs $(nc,nd)$ with n a positive integer and $(c,d)$ mutually prime coincides with the set of all pairs of integers $(m,k)$ which are not simultaneously zero. Indeed, given a pair $(m,k)$ we can factor out the greatest common divisor n and represent it as $(nc,nd)$ with mutually prime $(c,d)$. Thus we arrive at the Eisenstein series representation of the wave function:

$$\begin{array}{c}\hfill \zeta \left(2s\right){\psi }_s\left(z\right)=\frac{1}{2}\sum _{\genfrac{}{}{0pt}{}{(m,k)\in {\mathbb{Z}}^2}{(m,k)\ne (0,0)}}\frac{y^s}{{({(mx+k)}^2+m^2{y}^2)}^s}.\end{array}$$

(36)

Since the r.h.s. of this equation is periodic in x with period 1, we can expand it in Fourier series. Our next goal is to find the coefficients of this expansion:

$$\begin{array}{c}\hfill c_l\left(y\right)=\frac{1}{2}\sum _{\genfrac{}{}{0pt}{}{(m,k)\in {\mathbb{Z}}^2}{(m,k)\ne (0,0)}}{\int }_0^1\frac{y^s{e}^{-2\pi ilx}dx}{{({(mx+k)}^2+m^2{y}^2)}^s}.\end{array}$$

(37)

First let us handle the term with $m=0$:

$$\begin{array}{c}\hfill \frac{1}{2}\sum _{\genfrac{}{}{0pt}{}{k\in \mathbb{Z}}{k\ne 0}}{\int }_0^1\frac{y^s{e}^{-2\pi ilx}dx}{k^{2s}}={\delta }_{l,0}\zeta \left(2s\right)y^s.\end{array}$$

(38)

For the sum over non-zero m’s let us note that we may drop the factor $1/2$ and sum over $m\ge 1$. Indeed, the sum over negative m’s can be reverted to a sum over positive ones through redefinition $k\to -k$. For fixed positive m it is instructive to represent k as $k=nm+s$ thus splitting the initial sum over $k\in \mathbb{Z}$ into double sum over $r=0,1,\dots ,m-1$ and $n\in \mathbb{Z}$. In this way after few simple manipulations we get

$$\begin{array}{c}\sum _{m=1}^{\infty }\sum _{r=0}^{m-1}\sum _{n\in \mathbb{Z}}{\int }_0^1\frac{y^s{e}^{-2\pi ilx}dx}{{({(m(x+n)+r)}^2+m^2{y}^2)}^s}=\sum _{m=1}^{\infty }\sum _{r=0}^{m-1}{\int }_{-\infty }^{\infty }\frac{y^s{e}^{-2\pi ilx}dx}{{({(mx+r)}^2+m^2{y}^2)}^s}\hfill \\ =\sum _{m=1}^{\infty }\sum _{r=0}^{m-1}{m}^{-2s}{y}^{1-s}{e}^{\frac{2\pi ilr}{m}}{\int }_{-\infty }^{\infty }\frac{cos\left(2\pi \right|l\left|yx\right)dx}{{(x^2+1)}^s}.\hfill \end{array}$$

(39)

The last integral is expressed in terms of modified Bessel K function:

$$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }\frac{cos\left(2\pi \right|l\left|yx\right)dx}{{(x^2+1)}^s}=\left\{\begin{array}{c}\frac{2{\pi }^s}{\Gamma \left(s\right)}{\left|ly\right|}^{s-\frac{1}{2}}{K}_{s-\frac{1}{2}}\left(2\pi \right|l\left|y\right),\mathrm{i}fl\ne 0\hfill \\ \frac{\sqrt{\pi }\Gamma \left(s-\frac{1}{2}\right)}{\Gamma \left(s\right)},\mathrm{if}l=0.\hfill \end{array}\right.\end{array}$$

(40)

A useful alternative representation of modified Bessel K function which makes its properties more transparent is given by

$$K_{iu}\left(y\right)=\frac{1}{2}{\int }_{-\infty }^{\infty }{e}^{-ycosht}{e}^{iut}dt.$$

(41)

This expression allows analytical continuation of the wave function from the region $\Re s>1$ in (34) into the whole complex plane s because the Bessel $K_s\left(y\right)$ functions are well defined for any s. Besides, an easy examination shows that the finite sum is:

$$\begin{array}{c}\hfill \sum _{r=0}^{m-1}{e}^{\frac{2\pi ilr}{m}}=\left\{\begin{array}{c}m\mathrm{if}m\mathrm{divides}l\hfill \\ 0\mathrm{otherwise}\hfill \end{array}.\right.\end{array}$$

(42)

To summarise, for the Fourier coefficients (37) we shall get

$$\begin{array}{c}\hfill c_l\left(y\right)=\frac{2{\pi }^s}{\Gamma \left(s\right)}{\tau }_{s-\frac{1}{2}}\left(\right|l\left|\right)\sqrt{y}{K}_{s-\frac{1}{2}}\left(2\pi \right|l\left|y\right),\mathrm{if}l\ne 0,\end{array}$$

(43)

where

$$\begin{array}{c}\hfill {\tau }_{\nu }\left(n\right)=\sum _{a·b=n}{\left(\frac{a}{b}\right)}^{\nu },\end{array}$$

(44)

while for $l=0$:

$$\begin{array}{c}\hfill c_0\left(y\right)=\frac{\sqrt{\pi }\Gamma (s-\frac{1}{2})\zeta (2s-1)}{\Gamma \left(s\right)}{y}^{1-s}.\end{array}$$

(45)

Thus we recovered the second solution $y^{1-s}$ in (31) and calculated the coefficient $c_0^{{}^{\prime }}$ in front of it. Thus the invariant solution (36) takes the following form:

$$\begin{array}{lll}\zeta \left(2s\right){\psi }_s(x,y)& =& \zeta \left(2s\right)y^s+\frac{\sqrt{\pi }\Gamma (s-\frac{1}{2})\zeta (2s-1)}{\Gamma \left(s\right)}{y}^{1-s}+\\ & +& \sqrt{y}\frac{4{\pi }^s}{\Gamma \left(s\right)}\sum _{l=1}^{\infty }{\tau }_{s-\frac{1}{2}}\left(l\right)K_{s-\frac{1}{2}}\left(2\pi ly\right)cos\left(2\pi lx\right).\end{array}$$

(46)

Using Riemann’s reflection relation

$$\begin{array}{c}\hfill \zeta \left(s\right)=\frac{{\pi }^{s-\frac{1}{2}}\Gamma \left(\frac{1-s}{2}\right)}{\Gamma \left(\frac{s}{2}\right)}\zeta (1-s)\end{array}$$

(47)

and introducing the notation

$$\begin{array}{c}\hfill \theta \left(s\right)={\pi }^{-s}\zeta \left(2s\right)\Gamma \left(s\right)\end{array}$$

(48)

we arrive at the elegant final expression for the energy eigenfunctions obtained by Maas [37]:

$$\begin{array}{ccc}\hfill \theta \left(s\right){\psi }_s\left(z\right)& =& \theta \left(s\right)y^s+\theta (1-s)y^{1-s}+4\sqrt{y}\sum _{l=1}^{\infty }{\tau }_{s-\frac{1}{2}}\left(l\right)K_{s-\frac{1}{2}}\left(2\pi ly\right)cos\left(2\pi lx\right),\hfill \end{array}$$

(49)

This wave function is well defined in the complex s plane and has a simple pole at $s=1$. The physical continuous spectrum was defined in (25), where $s=\frac{1}{2}+iu$, $u\in [0,\infty ]$ so that

$$E=s(1-s)=\frac{1}{4}+u^2.$$

(50)

The continuous spectrum wave functions ${\psi }_s(x,y)$ are delta function normalisable [37,38,39,40,41,42]. The wave function (49) can be conveniently represented also in the form

$$\begin{array}{ccc}\hfill {\psi }_{\frac{1}{2}+iu}\left(z\right)& =& y^{\frac{1}{2}+iu}+\frac{\theta (\frac{1}{2}-iu)}{\theta (\frac{1}{2}+iu)}{y}^{\frac{1}{2}-iu}+\frac{4\sqrt{y}}{\theta (\frac{1}{2}+iu)}\sum _{l=1}^{\infty }{\tau }_{iu}\left(l\right)K_{iu}\left(2\pi ly\right)cos\left(2\pi lx\right),\hfill \end{array}$$

(51)

where

$$\begin{array}{c}\hfill K_{-iu}\left(y\right)=K_{iu}\left(y\right),{\tau }_{-iu}\left(l\right)={\tau }_{iu}\left(l\right).\end{array}$$

(52)

The physical interpretation of the wave function becomes more transparent when we introduce the new variables

$$\tilde{y}=lny,p=-u,E=p^2+\frac{1}{4}$$

(53)

as well as an alternative normalisation of the wave function:

$$\begin{array}{c}{\psi }_p(x,\tilde{y})\equiv y^{-\frac{1}{2}}{\psi }_{\frac{1}{2}+iu}\left(z\right)=\\ =e^{-ip\tilde{y}}+\frac{\theta (\frac{1}{2}+ip)}{\theta (\frac{1}{2}-ip)}{e}^{+ip\tilde{y}}+\frac{4}{\theta (\frac{1}{2}-ip)}{\sum }_{l=1}^{\infty }{\tau }_{ip}\left(l\right)K_{ip}\left(2\pi l{e}^{\tilde{y}}\right)cos\left(2\pi lx\right).\end{array}$$

(54)

Indeed, the first two terms describe the incoming and outgoing plane waves. The plane wave $e^{-ip\tilde{y}}$ incoming from infinity of the y axis in Figure 1 and Figure 2 (the vertex $\mathcal{D}$) elastically scatters on the boundary $ACB$ of the fundamental region $\mathcal{F}$ in Figure 1. The reflection amplitude is a pure phase and is given by the expression in front of the outgoing plane wave $e^{ip\tilde{y}}$:

$$\frac{\theta (\frac{1}{2}+ip)}{\theta (\frac{1}{2}-ip)}=exp\left[i\phi \right(p\left)\right].$$

(55)

The rest of the wave function describes the standing waves $cos\left(2\pi lx\right)$ in the x direction between boundaries $x=\pm 1/2$ with the amplitudes $K_{ip}\left(2\pi ly\right)$, which are exponentially decreasing with index l.

In addition to the continuous spectrum, the system (24) has a discrete spectrum [37,38,39,40,41,42]. The number of discrete states is infinite: $E_0=0<E_1<E_2<\dots \to \infty$, the spectrum is extended to infinity—unbounded from above—and lacks any accumulation point except infinity. Let us denote the wave functions of the discrete spectrum by ${\psi }_n\left(z\right)$ so that the expansion into the full set of basis vectors will take the form

$$f(x,\tilde{y})=\sum _{n\ge 0}{a}_n{\psi }_n(x,\tilde{y})+\frac{1}{2\pi }{\int }_0^{\infty }{a}_p{\psi }_p(x,\tilde{y})dp$$

(56)

and the Parseval identity will be

$$\begin{array}{c}{\left|\right|f\left(z\right)\left|\right|}^2={\sum }_{n\ge 0}|a_n{|}^2+\frac{1}{2\pi }{\int }_0^{\infty }{\left|a_p\right|}^{2}dp,\end{array}$$

(57)

and

$$\begin{array}{c}{\sum }_{n\ge 0}{\psi }_n(x,\tilde{y}){\psi }_n^{*}(x_1,{\tilde{y}}_1)+\frac{1}{2\pi }{\int }_0^{\infty }{\psi }_p(x,\tilde{y}){\psi }_{-p}(x_1,{\tilde{y}}_1)dp={\delta }^{\left(2\right)}(z-z_1).\end{array}$$

(58)

The wave functions of the discrete spectrum have the following form [37,38,39,40,45,46,47]:

$$\begin{array}{ccc}\hfill {\psi }_n\left(z\right)& =& \sum _{l=1}^{\infty }{c}_l\left(n\right)\sqrt{y}{K}_{i{u}_n}\left(2\pi ly\right)\left\{\begin{array}{c}cos\left(2\pi lx\right)\hfill \\ sin\left(2\pi lx\right)\hfill \end{array},\right.\hfill \end{array}$$

(59)

where the spectrum $E_n=\frac{1}{4}+u_n^2$ and the coefficients $c_l\left(n\right)$ are not known analytically, but were computed numerically for many values of n [45,46,47]. Having explicit expressions of the wave functions one can analyze the quantum–mechanical behaviour of the correlation functions, which we shall investigate in the next sections.

4. Correlation Functions

First let us calculate the two-point correlation function:

$$\begin{array}{c}\hfill {\mathcal{D}}_2(\beta ,t)=\langle A\left(t\right)B\left(0\right)e^{-\beta H}\rangle =\sum _n\langle n|e^{iHt}A\left(0\right)e^{-iHt}B\left(0\right)e^{-\beta H}|n\rangle =\\ \hfill =\sum _{n,m}{e}^{i(E_n-E_m)t-\beta E_n}\langle n\left|A\left(0\right)\right|m\rangle \langle m\left|B\left(0\right)\right|n\rangle .\end{array}$$

(60)

The energy eigenvalues (50) are parametrised by $n\to \frac{1}{2}+iu$, $E_n=\frac{1}{4}+u^2$ and $m\to \frac{1}{2}+iv$, $E_m=\frac{1}{4}+v^2$, thus

$$\begin{array}{c}\hfill {\mathcal{D}}_2(\beta ,t)={\int }_0^{+\infty }{\int }_0^{+\infty }dudv{e}^{i(u^2-v^2)t-\beta (\frac{1}{4}+u^2)}\\ \hfill {\int }_{\mathcal{F}}{\psi }_{\frac{1}{2}-iu}\left(z\right)A{\psi }_{\frac{1}{2}+iv}\left(z\right)d\mu \left(z\right){\int }_{\mathcal{F}}{\psi }_{\frac{1}{2}-iv}\left(w\right)B{\psi }_{\frac{1}{2}+iu}\left(w\right)d\mu \left(w\right),\end{array}$$

(61)

where the complex conjugate function is ${\psi }_{\frac{1}{2}+iu}^{*}\left(z\right)={\psi }_{\frac{1}{2}-iu}\left(z\right)$. Defining the basic matrix element as

$$\begin{array}{c}\hfill A_{uv}={\int }_{\mathcal{F}}{\psi }_{\frac{1}{2}-iu}\left(z\right)A{\psi }_{\frac{1}{2}+iv}\left(z\right)d\mu \left(z\right)={\int }_{-1/2}^{1/2}dx{\int }_{\sqrt{1-x^2}}^{\infty }\frac{dy}{y^2}{\psi }_{\frac{1}{2}-iu}\left(z\right)A{\psi }_{\frac{1}{2}+iv}\left(z\right)\end{array}$$

(62)

for the two-point correlation function we shall get

$$\begin{array}{c}\hfill {\mathcal{D}}_2(\beta ,t)={\int }_{-\infty }^{+\infty }{e}^{i(u^2-v^2)t-\beta (\frac{1}{4}+u^2)}{A}_{uv}{B}_{vu}dudv.\end{array}$$

(63)

In terms of the new variables (53) the basic matrix element (62) will take the form

$$\begin{array}{c}\hfill A_{pq}={\int }_{-1/2}^{1/2}dx{\int }_{\frac{1}{2}log(1-x^2)}^{\infty }dy{\psi }_p^{*}(x,y)\left(e^{-\frac{1}{2}y}A{e}^{\frac{1}{2}y}\right){\psi }_q(x,y).\end{array}$$

(64)

The matrix elements (62) and (64) play a fundamental role in the investigation of the correlation functions because all correlations can be expressed through it. One should also choose appropriate observables A and B. The operator $y^{-2}$ seems very appropriate for two reasons. Firstly, the convergence of the integrals over the fundamental region $\mathcal{F}$ will be well defined. Secondly, this operator is reminiscent of the exponentiated Liouiville field since $y^{-2}=e^{-2\tilde{y}}$. Thus we are interested in calculating the matrix element (64) for the observables in the form of the Liouiville-like operators 4:

$$A\left(N\right)=e^{-2Ny}.$$

(65)

We shall get

$$\begin{array}{c}\hfill A_{pq}\left(N\right)={\int }_{-1/2}^{1/2}dx{\int }_{\frac{1}{2}log(1-x^2)}^{\infty }dy{\psi }_p^{*}(x,y)e^{-2Ny}{\psi }_q(x,y),N=1,2,\dots \end{array}$$

(66)

To calculate the above matrix elements we shall use a perturbative expansion in which the part of the wave function (54) containing the Bessel functions and the contribution of the discrete spectrum (59) will be considered as perturbation. As we shall demonstrate below, these terms of the perturbative expansion are small and do not influence our results. The reason behind this fact is that in the integration region $\Im z\gg 1$ of the matrix element (62) the Bessel function decays exponentially, as one can be convinced by inspecting (41). Therefore, the contribution of these high modes is small (analogues to the so called mini-superspace approximation in the Liouville theory). Thus we shall consider the perturbative expansion over high frequency modes $l=1,2,...$ in (54) and (59). In the first approximation of the wave function (54) for the matrix element we shall get

$$\begin{array}{lll}{A}_{pq}\left(N\right)& =& {\int }_{-1/2}^{1/2}dx\left(\frac{{(1-x^2)}^{-N+\frac{p-q}{2i}}}{2N+i(p-q)}+\frac{{(1-x^2)}^{-N+\frac{p+q}{2i}}{e}^{-i\phi \left(q\right)}}{2N+i(p+q)}\right.\\ & & \left.+\frac{{(1-x^2)}^{-N-\frac{p+q}{2i}}{e}^{i\phi \left(p\right)}}{2N-i(p+q)}+\frac{{(1-x^2)}^{-N-\frac{p-q}{2i}}{e}^{i\left(\phi \right(p)-\phi (q\left)\right)}}{2N-i(p-q)}\right),\end{array}$$

(67)

where we used (55). Integration over x can be performed exactly with the result for the basic matrix element of the following form:

$$\begin{array}{c}\hfill A_{pq}\left(N\right)=\frac{{}_2{F}_1\left(\frac{1}{2},N+i\frac{p-q}{2};\frac{3}{2};\frac{1}{4}\right)}{2N+i(p-q)}+\frac{{}_2{F}_1\left(\frac{1}{2},N+i\frac{p+q}{2};\frac{3}{2};\frac{1}{4}\right)}{2N+i(p+q)}{e}^{-i\phi \left(q\right)}+\\ \hfill +\frac{{}_2{F}_1\left(\frac{1}{2},N-i\frac{p+q}{2};\frac{3}{2};\frac{1}{4}\right)}{2N-i(p+q)}{e}^{i\phi \left(p\right)}+\frac{{}_2{F}_1\left(\frac{1}{2},N-i\frac{p-q}{2};\frac{3}{2};\frac{1}{4}\right)}{2N+i(p-q)}{e}^{i\left(\phi \right(p)-\phi (q\left)\right)},\end{array}$$

(68)

where the reflation phase $\phi \left(p\right)$ was defined in (55). Thus for the two-point correlation function we shall get

$$\begin{array}{c}\hfill {\mathcal{D}}_2(\beta ,t)={\int }_{-\infty }^{+\infty }{e}^{i(p^2-q^2)t-\beta (\frac{1}{4}+p^2)}{A}_{pq}\left(N\right)A_{qp}\left(M\right)dpdq.\end{array}$$

(69)

The correlation function is between two Liouiville-like fields in the power N and M respectively. This expression is very convenient for analytical and numerical analyses. It is expected that the two-point correlation function should decay exponentially [26], therefore we fitted the points by the curve

$${\mathcal{D}}_2(\beta ,t)\sim K_2\left(\beta \right)e^{-\frac{t}{t_2\left(\beta \right)}}.$$

(70)

The $t_2\left(\beta \right)$ is the dissipation time and defines one of the characteristic time scales in the quantum–mechanical system5. The exponential decay of the two-point correlation function with time and the dependence of the exponent $t_2\left(\beta \right)$ and of the prefactor $K_2\left(\beta \right)$ as a function of temperature are presented in Figure 3. As one can see, at high and low temperatures the dissipation time tends to the fixed values. These limiting values we have calculated by using the expressions (71) and (72). At large $\beta \to \infty$ the main contribution came from the zero momentum region $p=0$:

$$\begin{array}{c}\hfill {\mathcal{D}}_2(\infty ,t)={\int }_{-\infty }^{+\infty }{e}^{-i{q}^{2}t}{A}_{0q}\left(N\right)A_{q0}\left(M\right)dq\sim K_2(\infty )e^{-\frac{t}{t_2(\infty )}},\end{array}$$

(71)

and at $\beta \to 0$ from the region $p=q$:

$$\begin{array}{c}\hfill {\mathcal{D}}_2(0,t)={\int }_{-\infty }^{+\infty }{e}^{i(p^2-q^2)t}{A}_{pq}\left(N\right)A_{qp}\left(M\right)dpdq\sim K_2\left(0\right)e^{-\frac{t}{t_2\left(0\right)}}.\end{array}$$

(72)

The corresponding limiting values in dimensionless units are shown on the l.h.s. of Figure 3.

Four-Point Correlation Function

It was conjectured by Maldacena, Shenker and Stanford [26] that the classical chaos can be diagnosed in thermal quantum systems by using an out-of-time-order correlation functions as well as by the square of the commutator of the operators which are separated in time. The out-of-time four-point correlation function of interest was defined in [26,48] as follows:

$$\begin{array}{c}\hfill {\mathcal{D}}_4(\beta ,t)=\langle A\left(t\right)B\left(0\right)A\left(t\right)B\left(0\right)e^{-\beta H}\rangle =\\ \hfill =\sum _{n,m,l,r}{e}^{i(E_n-E_m+E_l-E_r)t-\beta E_n}\langle n\left|A\left(0\right)\right|m\rangle \langle m\left|B\left(0\right)\right|l\rangle \langle l\left|A\left(0\right)\right|r\rangle \langle r\left|B\left(0\right)\right|n\rangle .\end{array}$$

(73)

The other important observable which we shall consider here is the square of the commutator of the Liouiville-like operators separated in time [26]

$$C(\beta ,t)=\langle {[A\left(t\right),B\left(0\right)]}^2{e}^{-\beta H}\rangle .$$

(74)

The energy eigenvalues we shall parametrise as $n\to \frac{1}{2}+iu$, $m\to \frac{1}{2}+iv$, $l\to \frac{1}{2}+il$ and $r\to \frac{1}{2}+ir$, thus from (73) we shall get

$$\begin{array}{ccc}\hfill {\mathcal{D}}_4(\beta ,t)& =& {\int }_{-\infty }^{+\infty }{e}^{i(u^2-v^2+l^2-r^2)t-\beta (\frac{1}{4}+u^2)}{A}_{uv}{B}_{vl}{A}_{lr}{B}_{ru}dudvdldr.\hfill \end{array}$$

(75)

In terms of the variables (53) the four-point correlation function (75) will take the following form:

$$\begin{array}{c}\hfill {\mathcal{D}}_4(\beta ,t)={\int }_{-\infty }^{+\infty }{e}^{i(p^2-q^2+l^2-r^2)t-\beta (\frac{1}{4}+p^2)}{A}_{pq}\left(N\right)A_{ql}\left(M\right)A_{lr}\left(N\right)A_{rp}\left(M\right)dpdqdldr.\end{array}$$

(76)

In order to compute the square of the commutator (74), one should also consider the following four-point correlation function:

$$\begin{array}{c}\hfill {\mathcal{D}}_4^{\prime }(\beta ,t)=\langle A\left(t\right)B\left(0\right)B\left(0\right)A\left(t\right)e^{-\beta H}\rangle =\\ \hfill =\sum _{n,m,l,r}{e}^{i(-E_m+E_r)t-\beta E_n}\langle n\left|A\left(0\right)\right|m\rangle \langle m\left|B\left(0\right)\right|l\rangle \langle l\left|B\left(0\right)\right|r\rangle \langle r\left|A\left(0\right)\right|n\rangle ,\end{array}$$

(77)

$$\begin{array}{c}\hfill {\mathcal{D}}_4^{″}(\beta ,t)=\langle B\left(0\right)A\left(t\right)A\left(t\right)B\left(0\right)e^{-\beta H}\rangle =\\ \hfill =\sum _{n,m,l,r}{e}^{i(E_m-E_r)t-\beta E_n}\langle n\left|B\left(0\right)\right|m\rangle \langle m\left|A\left(0\right)\right|l\rangle \langle l\left|A\left(0\right)\right|r\rangle \langle r\left|B\left(0\right)\right|n\rangle ,\end{array}$$

(78)

$$\begin{array}{c}\hfill {\mathcal{D}}_4^{‴}(\beta ,t)=\langle B\left(0\right)A\left(t\right)B\left(0\right)A\left(t\right)e^{-\beta H}>=\\ \hfill =\sum _{n,m,l,r}{e}^{i(-E_n+E_m-E_l+E_r)t-\beta E_n}\langle n\left|B\left(0\right)\right|m\rangle \langle m\left|A\left(0\right)\right|l\rangle \langle l\left|B\left(0\right)\right|r\rangle \langle r\left|A\left(0\right)\right|n\rangle .\end{array}$$

(79)

As it was suggested in [26], the most important correlation function indicating the traces of the classical chaotic dynamics in quantum regime is (74)

$$\begin{array}{c}\hfill C(\beta ,t)=-{\mathcal{D}}_4(\beta ,t)+{\mathcal{D}}_4^{\prime }(\beta ,t)+{\mathcal{D}}_4^{″}(\beta ,t)-{\mathcal{D}}_4^{‴}(\beta ,t),\end{array}$$

(80)

where in the case of the Artin system these correlation functions in the momentum space representation have the form:

$$\begin{array}{c}{\mathcal{D}}_4^{\prime }(\beta ,t)+{\mathcal{D}}_4^{″}(\beta ,t)=\\ 2{\int }_{-\infty }^{+\infty }{e}^{-\beta (\frac{1}{4}+p^2)}cos(q^2-r^2)t{A}_{pq}\left(N\right)A_{ql}\left(M\right)A_{lr}\left(N\right)A_{rp}\left(M\right)dpdqdldr,\end{array}$$

(81)

and

$$\begin{array}{c}{\mathcal{D}}_4(\beta ,t)+{\mathcal{D}}_4^{‴}(\beta ,t)=\\ 2{\int }_{-\infty }^{+\infty }{e}^{-\beta (\frac{1}{4}+p^2)}cos(p^2-q^2+l^2-r^2)t{A}_{pq}\left(N\right)A_{ql}\left(M\right)A_{lr}\left(N\right)A_{rp}\left(M\right)dpdqdldr.\end{array}$$

(82)

Figure 4 shows the behaviour of the four-point correlations as a function of temperature and time. All of the four-point correlation functions decay exponentially:

$${\mathcal{D}}_4(\beta ,t)\sim K_4\left(\beta \right)e^{-\frac{t}{t_4\left(\beta \right)}}.$$

(83)

Now we can turn to the investigation of the commutator (74), (80) which can be represented in the following form:

$$\begin{array}{c}\hfill C(\beta ,t)=2{\int }_{-\infty }^{+\infty }{e}^{-\beta (\frac{1}{4}+p^2)}\{cos(q^2-r^2)t-cos(p^2-q^2+l^2-r^2)t\}\\ \hfill A_{pq}\left(N\right)A_{ql}\left(M\right)A_{lr}\left(N\right)A_{rp}\left(M\right)dpdqdldr,\\ \hfill ={\int }_{-\infty }^{+\infty }{e}^{-\beta (\frac{1}{4}+p^2)}sin\frac{1}{2}(p^2+l^2-2r^2)t·sin\frac{1}{2}(p^2+l^2-2q^2)t\\ \hfill A_{pq}\left(N\right)A_{ql}\left(M\right)A_{lr}\left(N\right)A_{rp}\left(M\right)dpdqdldr,\end{array}$$

(84)

where we used (81) and (82). Due to the sinus functions the $C(\beta ,t)$ vanishes at $t=0$ and is proportional to $t^2$ at the very short time scale: $(p^2+l^2)t\ll 1,r^{2}t\ll 1,q^{2}t\ll 1$, $p^2\ll 1/\beta$. Generally the commutator $\left[A\right(0),A(t\left)\right]$ between operators at different moments of time is not vanishing because the operator $A\left(0\right)$ and the time dependent operator $A\left(t\right)=e^{iHt}A\left(0\right)e^{-iHt}$ do not coincide, for example in the case of a simple oscillator, where $[x\left(0\right),x\left(t\right)]=\frac{i\hslash }{m\omega }sin\left(\omega t\right)$. For larger time intervals the value of the matrix elements $A_{pq}$ (66) plays a crucial role and essentially influences the behaviour of the $C(\beta ,t)$.

In contrast to the exponential decay of the correlation functions the square of the commutator of the operators $e^{-Ny\left(t\right)}$ separated in time grows exponentially, Figure 5. This growth is reminiscent to the local exponential divergence of trajectories in the Artin system when it is considered in the classical [30,31] and in quasi-classical limits. In the quasi-classical limit the commutator $C(\beta ,t)$ changes proportionally to the square of trajectory deviation $\langle {\left(\frac{\delta y\left(t\right)}{\delta y\left(0\right)}\right)}^2\rangle$ which grows exponentially [30,31]. The growth of the commutator has been parameterised in the following form:

$$C(\beta ,t)\sim C\left(\beta \right)e^{\frac{2\pi }{\chi \left(\beta \right)}t},C{(\beta ,t)}_{max}/C{(\beta ,t)}_{Artin}\sim R\left(\beta \right)e^{(\frac{1}{\beta }-\frac{1}{\chi \left(\beta \right)})2\pi t},$$

(85)

where the behaviour of the exponent $\chi \left(\beta \right)$ and ratio with respect to the maximal growth exponent is presented in Figure 5. This beautifully confirms the fact that the correlation function $C(\beta ,t)$ indeed grows exponentially with time, as it takes place in its classical counterpart. As one can see, the exponent $\chi \left(\beta \right)$ defining the behaviour of the commutator $C(\beta ,t)$ in (84) and (85) slowly decreases with $\beta$.

5. Discussion

In order to check if the results are sensitive to the truncation of the high modes of the Maass wave function (49), as well as to the contribution of the discrete energy levels (59), we included the corresponding terms into the integration of the basic matrix element $A_{uv}$ in (62). The first discrete energy level is at $u_1\approx 9.53$ and has the wave function given in (59) [46]. The value of the corresponding matrix element (62) at $u=v=u_1$ is of order ${10}^{-15}$. In comparison, the contribution from continuous energy spectrum to the same matrix element is varying in the interval $[0.01-1]$. We also included two consecutive terms of the discrete sum of the Maass wave function (54) with $l=1,2$ in our calculation of the basic matrix element $A_{uv}$ and found that the changes are of the same order. Thus we found that their influence on the behaviour of the correlation functions is negligible mainly because the Bessel function (41) is exponentially small in the fundamental region. The numerical values of the exponents are changing in the range of a few percentage points and do not influence the results. In summary, all two- and four-point correlation functions decay exponentially. The commutator $C(\beta ,t)$ in (84) and (85) grows exponentially with an exponent that is almost constant, Figure 5.