Correlation Functions of Quantum Artin System

: 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 deﬁned on the fundamental region of the modular group SL(2,Z) represents a well deﬁned example of a highly chaotic dynamical system in its classical regime. We investigated the inﬂuence 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.


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 curvature 1 . 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 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]. 1 The definition of the Anosov C-systems can be found in his excellent work [1]. In order to provide this definition one should use such mathematical concepts as tangent vector bundle, derivative mapping, contracting and expanding linear spaces and others. The properties of the C-systems [1], of the Kolmogorov entropy [2][3][4] and the description of their periodic trajectories can be found in the recent review article [5]. The non-compact fundamental region F of a finite area is represented by the hyperbolic triangle ABD. The vertex D is at infinity of the y axis and corresponds to a cusp. The edges of the triangle are the arc AB, the rays AD and BD. The points on the edges AD and BD and the points of the arcs AC with CB should be identified by the transformations w = z + 1 and w = −1/z in order to form a closed non-compact surfaceF by "gluing" the opposite edges of the modular triangle together. The hyperbolic triangle OAB can be considered equally well as the fundamental region. The modular transformations of the fundamental region F create a regular tessellation of the whole Lobachevsky plane by congruent hyperbolic triangles. The trajectory K is geodesic passing through the point (x, y) of F in the v direction.
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 : D t (A, B) = M A(g)B(gg t )dµ(g). (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 F . The exponential decay rate was expressed in terms of the entropy h(F ) of the system: 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ỹ for the distance in the vertical direction 2 on the fundamental triangle dx=0 ds = dy/y = ln y =ỹ and the corresponding momentum p y in order to represent the Maass wave functions (54) in a physically insightful manner where the functions θ(s), τ iu (l) and K iu (y) 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 −ip yỹ incoming from infinity of the y axis in Figure 2 (the vertex D) elastically scatters on the boundary ACB of the fundamental triangle 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 ip yỹ : The rest of the wave function describes the standing waves cos(2πlx) in the x direction between boundaries x = ±1/2 with the amplitudes K ip y (2πly), which are exponentially decreasing with index l Figure 2. The continuous energy spectrum is given by the formula 2 Physical distances are obtained by integration of the metric (12) along the trajectories.
e ipỹ e −ipỹ D Figure 2. The incoming and outgoing plane waves (3). The plane wave e −ipỹ incoming from infinity of the y axis (the vertex D) elastically scatters on the boundary ACB of the fundamental triangle ABD which defines the fundamental region F (see also Figure 1). The reflection amplitude is a pure phase and is given by the expression in front of the outgoing plane wave e ipỹ . The rest of the wave function describes the standing waves in the x direction between boundaries x = ±1/2 with the amplitudes, which are exponentially decreasing.
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]: where in our case the operators A and B are chosen to be of the Liouiville type: lAnalyzing 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: with exponents t 2 (β) and t 4 (β) 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 Figures 3 and 4.  Figure 3. The exponential decay of the two-point correlation function D 2 (β, t) as a function of time at temperature β = 1. The points are fitted by the curve K 2 (β) exp (−t/t 2 (β)). The exponent t 2 (β) has a well defined high and low temperature limits. The limiting values in dimensionless units are t 2 (0) ≈ 0.276 and t 2 (∞) ≈ 0.749. The temperature dependence of K 2 (β) is shown at the far right.
). The temperature dependence of the exponent t 4 (β) have a well defined high and low temperature limits and is shown at far right. The corresponding limiting values of the function t 4 (β) in dimensionless units are t 4 (0) = 0.112 and t 4 (∞) = 0.163. The behaviour of the exponent t 2 (β) of the two-point correlation function is shown in Figure 3.
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: where the behaviour of the exponent χ(β) is presented in Figure 5. The temperature dependence of the exponent χ(β) 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.

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 ∈ C, z > 0} supplied with the metric (we set z = x + iy) which has Ricci scalar R = −2. Isometries of this space are given by SL(2, R) transformations.
The SL(2, R) matrix (a,b,c,d are real and ad − bc = 1) acts on a point z by linear fractional substitutions: Note also that g and −g give the same transformation, hence the effective group is SL(2, R)/Z 2 . We'll be interested in the space of orbits of a discrete subgroup G ⊂ SL(2, R) in H. Our main example will be the modular group G = SL(2, Z). A nice choice of the fundamental region F of SL(2, Z) is displayed in Figure 1. The fundamental region F of the modular group consists of those points between the lines x = − 1 2 and x = + 1 2 that lie outside the unit circle in Figure 1. The modular triangle F has two equal angles α = β = π 3 and the third one at D equal to zero, γ = 0, thus α + β + γ = 2π/3 < π. The area of the fundamental region is finite and equal to π 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µ(z) = dxdy y 2 , thus Area(F ) = Consider geodesic flow on F , which is conveniently described by the least action principle δS = 0, where (cf. (12)): By varying the action, we get the equations of motion Notice the invariance of the action and of the equations of motion under time reparametrizations t → t(τ). 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 conditioṅ where H is a constant. In this gauge the Equation (15) will take the following form [49]: Defining the canonical momenta p x , p y conjugate to the coordinates x, y as we shall get the geodesic Equation (17) in the Hamiltonian form: Indeed, after defining the Hamiltonian as the corresponding equations will take the following form: 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) and consider the (time independent) Schrödinger equation The resulting equation explicitly reads: 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 ψ(x, y) that fulfills the periodic boundary condition (26) It is convenient to introduce parametrization of the energy E = s(1 − s) and to rewrite this equation as As far as E is real and semi-positive and parametrisation is symmetric with respect to s ↔ 1 − s it follows that the parameter s should be chosen within the range One should impose the "periodic" boundary condition on the wave function with respect to the modular group in order to have the wave function which is properly defined on the fundamental regionF shown in Figure 1. Taking into account that the transformation T : z → z + 1 belongs to SL(2, Z), one has to impose the periodicity condition ψ(z) = ψ(z + 1). Thus we have a Fourier expansion Inserting this into Equation (24), for the Fourier component f n (y) we get For the case n = 0 the solution which exponentially decays at large y reads and for n = 0 one simply gets Thus the solution can be represented in the following form: where the coefficients c 0 , c 0 , 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 → −1/z : This functional equation defines the coefficients c 0 , c 0 , 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.

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 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 = 0: Let us denote by Γ ∞ the subgroup of Γ = SL(2, Z), generating shifts z → z + n, n ∈ Z. Explicitly the elements of Γ ∞ are given by 2 × 2 matrices: Since y s is already invariant with respect to Γ ∞ , we should perform summation over the conjugacy classes Γ ∞ \Γ. Let us define these conjugacy classes. If two SL(2, Z) matrices 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 s > 1. We used also the simple relation To simplify further the sum let us multiply both sides of the Equation (34) by [37] ∞ ∑ n=1 1 n 2s ≡ ζ(2s) so that we shall get 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. 3 The factor 1/2 below is introduced for removing the double degeneracy due to the fact that SL(2, Z) elements ±γ both act on z in the same way. 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: 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: First let us handle the term with m = 0: For the sum over non-zero m's let us note that we may drop the factor 1/2 and sum over m ≥ 1. Indeed, the sum over negative m's can be reverted to a sum over positive ones through redefinition k → −k. For fixed positive m it is instructive to represent k as k = nm + s thus splitting the initial sum over k ∈ Z into double sum over r = 0, 1, . . . , m − 1 and n ∈ Z. In this way after few simple manipulations we get The last integral is expressed in terms of modified Bessel K function: A useful alternative representation of modified Bessel K function which makes its properties more transparent is given by This expression allows analytical continuation of the wave function from the region s > 1 in (34) into the whole complex plane s because the Bessel K s (y) functions are well defined for any s. Besides, an easy examination shows that the finite sum is: To summarise, for the Fourier coefficients (37) we shall get where while for l = 0: Thus we recovered the second solution y 1−s in (31) and calculated the coefficient c 0 in front of it. Thus the invariant solution (36) takes the following form: Using Riemann's reflection relation and introducing the notation we arrive at the elegant final expression for the energy eigenfunctions obtained by Maas [37]: 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 = 1 2 + iu, u ∈ [0, ∞] so that The continuous spectrum wave functions ψ 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 where The physical interpretation of the wave function becomes more transparent when we introduce the new variablesỹ = ln y, p = −u, E = p 2 + 1 4 (53) as well as an alternative normalisation of the wave function: Indeed, the first two terms describe the incoming and outgoing plane waves. The plane wave e −ipỹ incoming from infinity of the y axis in Figures 1 and 2 (the vertex D) elastically scatters on the boundary ACB of the fundamental region 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ỹ : The rest of the wave function describes the standing waves cos(2πlx) in the x direction between boundaries x = ±1/2 with the amplitudes K ip (2πly), which are exponentially decreasing with index l.
The wave functions of the discrete spectrum have the following form [37][38][39][40][45][46][47]: where the spectrum E n = 1 4 + u 2 n and the coefficients c l (n) 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.

Correlation Functions
First let us calculate the two-point correlation function: The energy eigenvalues (50) are parametrised by n → 1 2 + iu, E n = 1 4 + u 2 and m → 1 2 + iv, where the complex conjugate function is ψ * 1 2 +iu (z) = ψ 1 2 −iu (z). Defining the basic matrix element as for the two-point correlation function we shall get In terms of the new variables (53) the basic matrix element (62) will take the form 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 F will be well defined. Secondly, this operator is reminiscent of the exponentiated Liouiville field since y −2 = e −2ỹ . Thus we are interested in calculating the matrix element (64) for the observables in the form of the Liouiville-like operators 4 : We shall get 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 z 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 where we used (55). Integration over x can be performed exactly with the result for the basic matrix element of the following form: where the reflation phase ϕ(p) was defined in (55). Thus for the two-point correlation function we shall get 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 The t 2 (β) is the dissipation time and defines one of the characteristic time scales in the quantum-mechanical system 5 . The exponential decay of the two-point correlation function with time and the dependence of the exponent t 2 (β) and of the prefactor K 2 (β) 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 β → ∞ the main contribution came from the zero momentum region p = 0: and at β → 0 from the region p = q: 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 5 We would like to thank Gabriel Poghosyan for numerical calculation of the two-point correlation function presented in The other important observable which we shall consider here is the square of the commutator of the Liouiville-like operators separated in time [26] The energy eigenvalues we shall parametrise as n → 1 2 + iu, m → 1 2 + iv, l → 1 2 + il and r → 1 2 + ir, thus from (73) we shall get In terms of the variables (53) the four-point correlation function (75) will take the following form: In order to compute the square of the commutator (74), one should also consider the following four-point correlation function: As it was suggested in [26], the most important correlation function indicating the traces of the classical chaotic dynamics in quantum regime is (74) where in the case of the Artin system these correlation functions in the momentum space representation have the form: 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: Now we can turn to the investigation of the commutator (74), (80) which can be represented in the following form: where we used (81) and (82). Due to the sinus functions the C(β, t) vanishes at t = 0 and is proportional to t 2 at the very short time scale: (p 2 + l 2 )t 1, r 2 t 1, q 2 t 1, p 2 1/β. Generally the commutator [A(0), A(t)] between operators at different moments of time is not vanishing because the operator A(0) and the time dependent operator A(t) = e iHt A(0)e −iHt do not coincide, for example in the case of a simple oscillator, where [x(0), x(t)] = ih mω sin(ωt). 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(β, t).
In contrast to the exponential decay of the correlation functions the square of the commutator of the operators e −Ny(t) 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(β, t) changes proportionally to the square of trajectory deviation ( δy(t) δy(0) ) 2 which grows exponentially [30,31]. The growth of the commutator has been parameterised in the following form: C(β, t) ∼ C(β) e 2π χ(β) t , C(β, t) max /C(β, t) Artin ∼ R(β) e where the behaviour of the exponent χ(β) and ratio with respect to the maximal growth exponent is presented in Figure 5. This beautifully confirms the fact that the correlation function C(β, t) indeed grows exponentially with time, as it takes place in its classical counterpart. As one can see, the exponent χ(β) defining the behaviour of the commutator C(β, t) in (84) and (85) slowly decreases with β.

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 ≈ 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(β, t) in (84) and (85) grows exponentially with an exponent that is almost constant, Figure 5.