{'fraction_non_alphanumeric': 0.036716077537058155, 'fraction_numerical': 0.012455883151436174, 'mean_word_length': 4.4585655008891525, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Classical nonlinear response of a chaotic system: collective resonances

We consider the classical response in a chaotic system. In contrast to behavior in integrable or almost integrable systems, the nonlinear classical response in a chaotic system vanishes at long times. The response also reveals certain features of collective resonances which do not correspond to any periodic classical trajectories. The convergence of the response is shown to hold due to the exponential time dependence of the stability matrix. The growing exponentials corresponding to strong instability do not inhibit the convergence. We calculate both linear and second-order response in one of the simplest chaotic systems: free classical motion on a surface of constant negative curvature. We demonstrate the relevance of the model for applications to spectroscopic experiments.

(Dated: October 1, 2018) We consider the classical response in a chaotic system. In contrast to behavior in integrable or almost integrable systems, the nonlinear classical response in a chaotic system vanishes at long times. The response also reveals certain features of collective resonances which do not correspond to any periodic classical trajectories. The convergence of the response is shown to hold due to the exponential time dependence of the stability matrix. The growing exponentials corresponding to strong instability do not inhibit the convergence. We calculate both linear and second-order response in one of the simplest chaotic systems: free classical motion on a surface of constant negative curvature. We demonstrate the relevance of the model for applications to spectroscopic experiments.

I. INTRODUCTION
Time-resolved femtosecond spectroscopy constitutes a powerful tool that probes electronic and vibrational coherent dynamics of complex molecular systems in condensed phase [1,2,3,4,5,6]. Spectroscopic techniques allow for direct measurements of the time-and frequency-domain optical response functions that carry detailed information on the underlying dynamical phenomena. Nonlinear spectroscopy contains additional information on nonequilibrium processes absent in linear measurements.
In complex systems, such as polyatomic molecules, vibrational dynamics of many degrees of freedom can be adequately described in the framework of classical mechanics. In particular, at energies corresponding to room temperatures the complexity often originates from heavy strongly anharmonic vibrational modes that can be treated classically.
A number of studies have been devoted to the nonlinear response in stable (integrable) dynamical systems [7,8,9,10,11]. Classical nonlinear response functions have been shown to diverge with time in the general case [11], although in certain cases the divergence can be removed by thermal averaging using the canonical distribution [7,8]. The latter effect can be viewed as a result of destructive interference of different paths (dephasing), while the divergence is associated with the rephasing.
The divergence would pose a major problem for the whole concept of perturbative response. On the other hand, quantum response functions do not exhibit any divergences. However, application of the fully quantum description is often unmanageable in the systems of interest. It has been demonstrated that systematic expansion cures the divergence when the limit of vanishing * Electronic address: chernyak@chem.wayne.edu is taken after obtaining the asymptotics of long times [9,10]. Alternatively the divergence can be removed by adding noise that models interaction with a bath (or a solvent) [12].
Previous analytical studies of the nonlinear response were mostly limited to fully solvable models, which implies integrable dynamics in the classical limit. Integrable dynamics corresponds to an idealized picture for a realistic physical situation, since even weak perturbations generally destroy integrability. However, the problem is that weak deviations from the integrability do not eliminate the unphysical behavior of the response functions. This follows from the fact that the divergence in integrable systems originates from quasiperiodic motions on invariant tori [9]. According to the Kolmogorov-Arnold-Moser (KAM) theory, most invariant tori are not destroyed by small perturbations that break down integrability [13]. Therefore, the long-time behavior of the response in almost integrable systems is similar to its integrable counterpart.
Stable dynamics is typical for a situation close to equilibrium. At larger energies a generic situation corresponds to dynamical behavior with chaotic features [14]. Moreover, unstable (hyperbolic) dynamics may be more common due to the stability of chaos with respect to perturbations. It has been argued based on results of numerical analysis [15] that chaotic dynamics appears to observe the convergence of the classical response functions. In spite of apparent importance and to the best of our knowledge, the problem of the nonlinear response in strongly chaotic systems has never been addressed using analytical methods. We note that analytical calculations of the response are rarely feasible in nonintegrable systems, whereas numerical simulations of chaotic dynamics are complicated by the exponential divergence of stability matrices [15].
A strongly chaotic (mixing) system is characterized by a special spectrum of the Liouville operator [16,17,18]. The spectrum consists of complex Ruelle-Pollicott (RP) resonances that determine the asymptotic oscillations and decay of the correlations. This yields the linear response function directly related to the two-point correlation functions, in agreement with the fluctuation dissipation theorem (FDT). Nonlinear response, however, turns out to be more involved still with the noticeable effect of the resonances.
In this work we show that (i) the classical linear and nonlinear response of a chaotic system exhibits decay and oscillations as a function of times between the driving pulses, and (ii) the Fourier transform of 2D second-order response function reveals broad and asymmetric peaks that can be viewed as signatures of chaos in underlying dynamics.
We first establish a general qualitative picture of linear and nonlinear response in a classical system with hyperbolic chaotic dynamics and demonstrate that the classical response functions exponentially vanish at large times. To exemplify our arguments we further perform detailed analytical calculations of linear and second-order response functions for a free particle moving along a compact surface of constant negative curvature. The model, that constitutes a well-known example of classical chaos and a prototype for quantum chaos, allows for an exact solution due to strong dynamical symmetry. To the best of our knowledge, we present the first explicit analytical calculation of nonlinear response in a chaotic system.
The manuscript is organized as follows. We begin with reviewing the classical response theory using the Liouville representation of classical mechanics. In Section III we present the qualitative picture for the classical nonlinear response functions and argue that they show exponential decay at long times. We also argue that a generic finite motion of the molecular system with potential interactions is similar to the free motion along an effective compact configuration space of negative curvature. In Section IV we consider the chaotic model of free motion along a Riemann surface of constant negative curvature. We demonstrate how to make use of strong dynamical symmetry, and reduce the model to a much simpler problem of 1D dynamics on a circle. In Section V we present our calculation of the linear and second-order response functions. This calculation constitutes the main technical result of the paper. The results of numerical calculations of 2D spectroscopic signals are presented in Section VI. All necessary technical details are summarized in the Appendices.

II. CLASSICAL RESPONSE IN LIOUVILLE SPACE
In this section we review a general formalism of the classical response theory, using the language of probability distributions (e.g. see Ref. 15). We adopt the Liouville picture of classical mechanics where the system state that evolves in time is described by a distribution ρ(η, t), whereas the observables are represented by func-tions in phase space. This language is advantageous in a chaotic case, since mixing (e.g., hyperbolic) dynamics is characterized by strong instability in trajectory space, and individual trajectories do not provide meaningful information on the system relaxation.
Consider classical Hamilton dynamics that occurs in phase space M equipped with a Poisson bracket. Evolution of an externally driven system is described by the classical Liouville equation with a time-dependent Hamiltonian H T (t) ∂ρ(η, t) ∂t = −L T (t)ρ(η, t) , (1) where H(η) is the time independent Hamiltonian of the undriven system, E(t) is the time dependent uniform driving field, and dipole moment (polarization) f (η) describes the coupling of the system to the driving field. This sign choice of in the Poisson bracket {f, g} corresponds to a convention that the action of the Liouville operatorL T determines the phase space velocity aṡ η =L T (t)η.
In many cases f (η) is also an observable, e.g. polarization that creates spectroscopic signals, and the measured value is given by the time dependent average The undriven system with Hamiltonian H is characterized by equilibrium distribution ρ 0 (η) that depends only on the conserved energy H of the system, i.e. {H, ρ 0 } = 0. The average of f for the equilibrium distribution ρ 0 (η) normally vanishes (in spectroscopic literature this is often referred to as the absence of the permanent dipole): The response theory is usually based on a perturbative expansion of the measured signal in powers of the driving field. The response functions can be naturally obtained via calculating the correction to the equilibrium density distribution due to the driving field. The full Liouville equation is repeatedly solved, this results in an expansion ρ = ρ 0 + ρ 1 + . . . where ρ n represents the contribution of n-th order in the field E. Therefore, we start with solving the equation to get the linear in E correction to the density distribution: HereL is the Liouville operator of the unperturbed dynamics that corresponds to the Hamiltonian H. The corrections to ρ 0 in all orders of E can be obtained by solving equations similar to Eq. (5). This leads to the natural representation of the observable quantity in the form of an expansion which defines response functions S n (t; . . . , τ 1 ). The response function of order n depends on n + 1 time variables, it can be reduced to n time segments t n = t − τ n , . . . , t 1 = τ 2 − τ 1 in the case when the unperturbed dynamics has no explicit time dependence. For time segments representation, the linear and second-order response functions read: These expressions for the first-and second-order response will be the starting point of our general qualitative analysis of response functions in hyperbolic systems. They will be also utilized in our analytical calculations for a particular chaotic system where strong dynamical symmetry of the phase space enables us to solve the problem completely.

A. Integrable and almost integrable dynamics
Dynamics of a classical Hamiltonian system in the vicinity of an equilibrium (stationary) point can be adequately described by a system of uncoupled harmonic oscillators. This is achieved by expanding the classical Hamiltonian H up to second order in dynamical variables (the first-order terms vanish for an expansion around a stationary point). A harmonic system represents the simplest example of integrable dynamics.
The evolution of integrable systems with n degrees of freedom is naturally described in terms of n canonical pairs of action and angle variables. Action variables c j with j = 1, . . . , n are first integrals of motion {H, c j } = 0 in involution, i.e. {c i , c j } = 0, with the Hamiltonian H(η) = H(c 1 , . . . , c n ) depending on the integrals of motion only. The corresponding phase space vector fieldsĉ j defined byĉ j f = {c j , f }, therefore, commute and describe motions along the tori with frequencies ω j = ∂H(c 1 , . . . , c n )/∂c j that depend parametrically on the integrals of motion. Since for fixed values of the integrals of motion the angular velocities are constant, such dynamics is referred to as is conditionally periodic (or quasi-periodic) motion. For given values of action variables determined by initial conditions, the trajectories lie on the corresponding n-dimensional torus which is a subspace of 2n-dimensional phase space. The motion on the torus is strictly periodic if the ratios of corresponding angular velocities are rational (resonant torus); if the angular velocities are incommensurate, the trajectories densely cover the torus (non-resonant torus). Harmonic dynamics represents the simplest particular case of integrable dynamics when the frequencies ω i , i = 1, . . . , n are independent of the action variables c j .
The nonlinear response function for classical integrable dynamics have been shown to have power-like divergence at long times which can be eliminated by invoking a quantum description [7,8,9,10,11]. Divergence of nonlinear response in integrable systems does not imply unphysical behavior by itself, since integrable dynamics is an idealization. In real systems quantum effects, interactions with the bath or irregular dynamics may provide a necessary regularization. While the first two phenomena have been discussed in the literature, we concentrate on the latter. A small, yet generic, perturbation of the system Hamiltonian destroys integrability. However such a perturbation does not break down the linear divergence in the response functions. This follows from the qualitative picture of almost integrable dynamics established by the celebrated Kolmogorov-Arnold-Moser (KAM) perturbation theory [13]. The KAM theory states that a sufficiently small perturbation does not destroy most nonresonant tori, which means that in the invariant subspace of the entire phase space represented by the remaining distorted tori the dynamics preserves its quasiperiodic nature. Although motion inside the instability zones represented by destroyed tori becomes chaotic, their relative measure in phase space is small for small perturbations to the integrable dynamics. In the simplest n = 2 case the zones of instability are confined between remaining invariant tori.
It has been pointed out [8,9] that the divergence or the classical response functions in integrable systems originates from quasiperiodic nature of the underlying motion. Combined with the picture of almost integrable dynamics established by KAM theory, this demonstrates the unphysical divergence of the response functions is stable with respect to at least weak deviations from integrability. We will see that chaotic regions (instability zones) do not contribute to the divergence, therefore the diverging terms will decrease with increasing the deviation from the integrable situation due to decrease of the amount of surviving invariant tori.

B. Chaotic dynamics
Integrable and almost integrable Hamiltonian dynamics considered in subsection III A is typical for low energies when the system is moving in the neighborhood of a stable stationary point (equilibrium). Islands of stable dynamics completely vanish for larger deviations from the integrability caused by higher characteristic energies. The motion will be typically chaotic for higher temperatures or in essentially nonequilibrium processes such as photoinduced molecular dynamics on the excited electronic adiabatic surfaces.
In this subsection we present a qualitative picture of nonlinear response in chaotic systems. Due to FDT, the stability matrices were found to affect the response functions starting with the second order [15,19]. Therefore, the divergence of the nonlinear response functions in the integrable case can be attributed to the growth with time of certain stability matrix components. In the case of strong chaos there are components of the stability matrices that grow exponentially ∼ e λt with time, λ being the Lyapunov exponent. Numerical simulations have demonstrated that at least for some examples of chaotic dynamical systems the classical nonlinear response functions are free of unphysical divergences [15].
In this subsection we present qualitative arguments that rationalize physical exponentially decaying at long times behavior of classical nonlinear response functions in systems that exhibit strong enough chaotic behavior. More precisely we will consider mixing systems also known as A-flows, Smale or uniformly hyperbolic dynamical systems (for a nice overview without too many details see Ref. 20). In these systems at all points of the (2n − 1)-dimensional energy shell, or all relevant points (the relevant points belong to the so-called nonwandering subset) in the Smale's case, one can define n + unstable and n − stable tangent directions with n + + n − = 2n − 2, so that a small deviation from a trajectory along the stable (unstable) directions decays exponentially in time for forward (backward) dynamics. The stable (unstable) directions can be locally integrated to obtain stable (unstable) manifolds of dimension n − (n + ). Note that n + (n − ) is the number of positive (negative) Lyapunov exponents. Also note that in the Hamiltonian dynamics case the volume in phase space is conserved so that all Lyapunov exponents sum up to zero We will consider the simplest case of the lowest dimension coordinate space dimension n = 2 that allows for chaotic dynamics. The isoenergetic shell is 3dimensional, we have n + = n − = 1, and two nontrivial Lyapunov exponents ±λ.
A schematic picture of first-and second-order response formation in a hyperbolic chaotic system is shown in Fig. 1. Since chaotic dynamics described in terms of phase-space trajectories is extremely unstable, using the dual representation via evolution of distributions (Liouville representation) will be advantageous. We start with the linear response function S (1) (t). Since the equilibrium distribution depends on the system energy only, we can recast the expression for the linear response [Eq. (10)] in a form that actually represents the FDT Since dxf = 0 (no permanent dipole), we can consider the simplest set-up when the function f is represented by two small separate regions in the phase space where it adopts positive and negative values (an extension of our arguments to the general case is straightforward). By our assumption the linear sizes of the regions a ≪ l are small compared to the linear size l of the compact phase space. Evolution during time t changes the shapes of the regions. In the case of hyperbolic dynamics for t ≫ 1, the shape becomes similar to ribbon-like fettuccine: elongated along the unstable direction by a factor e λt , narrowed along the stable one by e −λt (we reiterate that λ is the positive the Lyapunov exponent) and unchanged along the flow. The long ribbon becomes evenly folded in the entire accessible energy shell. The overlap of the distributions f and e −Lt f ∂ E ρ 0 in Eq. (13) is represented by a large number N 1 (t) of disconnected regions. The dimensions of each region are typically a, ae −λt and a along the flow, in the stable and unstable direction, respectively. This gives the volume of the single 3D region in the energy shell v 1 (t) ∼ a 3 e −λt . A typical distance d(t) between the fettuccine in the overlap regions can be easily estimated as d(t) ∼ l 3 a −2 e −λt . This estimates the number of disconnected regions of the overlap as Since f assumes opposite signs in two initial regions, the distribution e −Lt f ∂ E ρ 0 consists of two positively and negatively "charged" fettuccine. The cancellations result in a signal determined by a typical fluctuation proportional to √ N 1 , and the overlap integral attains a factor √ N 1 v 1 . This results in the decaying long-time asymptotic of the linear response function S (1) (t) ∼ e −λt/2 .
We are now in a position to consider the second-order response function S (2) (t 1 , t 2 ). Propagating the observable f in Eq. (11) backward in time we interpret the second-order response as an overlap of the distributions f − (η) = exp(Lt 2 )f and ξ j ∂ j f + (η) with f + = exp(−Lt 1 )f ∂ E ρ 0 and the vector field ξ j ∂ j = {f, ·}. Similar to the case of linear response for t ≫ 1, the shapes of both distributions f ± (η) become similar to ribbon-like fettuccine, unchanged along the flow. Since f − results from reverse dynamics, the distributions f + and f − are elongated along the unstable and stable directions of the positive time evolution, by the factors e λt1 and e λt2 , respectively. Each of distributions f − and f + consists of two positively and negatively "charged" fettuccine originating from two oppositely signed separate regions of dipole moment distribution f . Since the vector field ξ introduced earlier is zero outside the support of f , the integration in Eq. (15) is restricted to the overlap of three distributions: f and f ± . The overlap of each of distributions f + and f − with f is represented by N 1 (t 1 ) and N 1 (t 2 ) disconnected fettuccine pieces as found in Eq. (14). Therefore the overlap of all three distributions is represented by N 2 (t 1 , t 2 ) ∼ N 1 (t 1 )N 1 (t 2 ) ∼ (l/a) 6 e −λ(t1+t2) disconnected regions. The linear dimensions along the isoenergetic shell of each disconnected region are estimated as ae −λt1 , ae −λt2 , and a, which corresponds to the stable, unstable, and flow directions, respectively. Therefore the volume of each disconnected region is v 2 (t) ∼ a 3 e −λ(t1+t2) .
If instead of the second-order response function, we were dealing with a three-point correlation function the rest of the story would be straightforward. A three-point correlation function can be represented in a form of the overlap integral I (2) (t 1 , t 2 ) with the differential operator ξ j ∂ j = {f, ·} replaced by the operator of multiplying by f (η). In full analogy with the linear response case, the typical value of the three-point correlation results from cancellations of oppositely signed contributions of similar magnitudes and therefore attains a factor √ N 2 v 2 ∼ e −λ(t1+t2)/2 that describes the long-time decaying behavior of the correlation function. The situation with the nonlinear response function is apparently more complicated since generally the vector field ξ has a component along the stable direction. The derivative ξ j ∂ j of a sharp feature in f + along the stable direction can create exponentially large e λt1 factors. This is the Liouville space signature of the exponentially growing components of the stability matrix, which affects the response starting with second order due to FDT [15,19]. The exponentially growing components of the stability matrix may seem to be a reason for an exponential divergence of the nonlinear response, since interaction with the driving field can be considered as a kick leading to an infinitesimal deviation that grows exponentially with time. This would actually happen if the initial distribution was δ-functional concentrated at some point in phase space, and the signal was measured by a deviation of the perturbed trajectory from its unperturbed counterpart. However, the dipole f that describes the system coupling to the driving field is represented by a smooth function. As shown below due to this smoothness the exponentially diverging terms cancel out completely. This demonstrates that response in chaotic systems should be treated using the Liouville (distribution-based) representation of classical dynamics, while its dual trajectory-based counterpart may lead to misleading naive picture.
To prove the harmlessness of the derivative in the second-order response function, we decompose the vector field ξ = ξ E + ξ 0 + ξ + + ξ − into the direction along the energy, flow, unstable, and stable components. This leads to a natural decomposition of the overlap integral, according to Eq. (15): The term I + involves a derivative along the stable direction that provides an additional exponentially small e −λt1 factor, since the distribution f + is elongated along the unstable direction. The dangerous term that involves derivatives along the sharp feature is represented by I (2) − . The aforementioned cancellation can be seen after the integration by parts: The first term in Eq. (17) includes time-independent distribution divξ − . The second term contains the derivative in the stable direction of the distribution f − . Therefore, it acquires an additional exponentially small factor e −λt2 and becomes negligible at long times. Finally, we have at long times All terms in Eq. (18) do not have derivatives in the stable and unstable directions that provide exponentially growing or decaying factors. Therefore, according to the arguments presented above, the second-order response function has a long-time asymptotic S (2) (t 1 , t 2 ) ∼ e −λ(t1+t2)/2 similar to the three-point correlation function. Stated differently, the exponentially growing components of the stability matrices that in principle enter the expressions for nonlinear response functions [15,19], do not affect the long-time behavior of the nonlinear response, due to the smoothness of the dipole function f (η) that describes the system coupling to the driving field. The qualitative arguments developed in this subsections for the lowest-dimensional case n = 2 can be extended to a general case n ≥ 2 in a straightforward way. This leads to the same asymptotic expressions for magnitudes of response functions S (1) (t) ∼ e −λt/2 , and S (2) (t 1 , t 2 ) ∼ e −λ(t1+t2)/2 , where λ should be interpreted in the sense of Eq. (12) as the sum of positive Lyapunov exponents.

C. Effective negative curvature of configuration space
In subsection III B we have established a qualitative picture of response in classical systems with hyperbolic and mixing dynamics. In the forthcoming sections we support our qualitative arguments by performing explicit analytical calculations of the first-and second-order response functions in a model system of a free particle moving along a Riemann surface of constant negative curvature. Due to strong dynamical symmetry in this system (see section IV for the details) the response functions can be calculated explicitly. In this subsection we present a rationale why this simple model system can be viewed as a representative example that bears the basic qualitative features of chaotic dynamics in a molecule or collection of molecules for high enough energies.
First of all classical motion of a molecular system can be represented as classical dynamics of a multidimensional particle moving in a potential U (r), where r stands for a complete set of nuclear coordinates. Trajectories of a particle with energy E in arbitrary potential U (r) are known to be the same as for a free motion in a curved space with the metric g ik = (1 − U (r)/E)δ ik [13,14]. Although the potential generates nonuniform motion along the trajectories, the dominant feature of the dynamics is the exponentially growing separation between the trajectories in the regions where the metric curvature is negative. The regions of the configuration space that correspond to negative curvature work as defocusing lenses, causing instability, i.e. divergence of close trajectories. For example, close trajectories diverge every time they approach the boundary of the classically inaccessible island or pass a region of a potential local maximum that belongs to the accessible region. Except for occasional symmetric integrable cases, passing through the stable regions with positive curvature cannot compensate for the instability, since stability reflects oscillatory, rather than converging features in the dynamics of the trajectory deviation. In particular, existence of unstable regions combined with ergodicity ensures exponential divergence of trajectories over long enough times.
In addition, when the motion is finite, the accessible part of the configuration space at a given energy can be multiply connected. In the simplest n = 2 case of two coordinates the motion occurs inside a disk-like region punctured by g forbidden islands (see Fig. 2). Boundaries of the accessible area are represented by curved lines where the potential energy U (r) coincides with the total energy E. Some fraction of trajectories approaches the boundaries so close that this can be qualified as reflection. Utilizing the original argument of Sinai [21] reflection can be interpreted as continuing motion on the antipode replica of the accessible region glued to its original via the boundaries (see Ref. 13). The resulting compact surface has topology of a sphere with g handles (Riemann surface of genus g). According to the Gauss-Bonnet theorem, in the g > 1 case the average Gaussian curvature is negative, which implies regions of instability and results in unstable (hyperbolic) dynamics.
We also note that loose distinction between reflected and deflected trajectories that approach the boundaries should not be a matter of concern since classical dynamics constitutes an approximation for a quantum mechanical problem. The approximation is not valid near the boundaries of the classically accessible region, where the quantum uncertainty takes over. More formally, one can consider a trajectory reflected (or, equivalently, contin-ued on the antipode replica of the accessible region) if it penetrates the near-boundary region where quantum effects become important. The classical picture of reflection becomes exact for the hard wall potential concentrated on the boundaries, e.g. for Sinai billiards [21] and Lorentz gas.
Although in a generic system the effective curvature depends on position r in configuration space, we will rationalize the described above qualitative picture of response in a chaotic system by calculating the linear and second-order classical response functions for free motion in a Riemann surface M 2 of constant negative (Gaussian) curvature. This model that allows an exact solution has been serving as a prototype of classical chaos, as well as an example of semiclassical quantization [13,14,22].

A. Free particle on a Riemann surface
In this section we describe classical dynamics of a free particle moving along a Riemann surface M 2 of constant negative curvature (Gaussian curvature). This will be done by making use of strong dynamical symmetry. This system is one of the best studied models of chaotic dynamics. Although detailed reviews of its symmetry exist in the literature (see e.g. Ref. 22), we summarize the necessary results for the sake of completeness and further extension to the case of Langevin dynamics. Noisy dynamics associated with geodesic flows on surfaces with constant negative curvature is considered in detail elsewhere [28].
The Lagrangian of a free particle contains only kinetic energy: L = mg ikṙ iṙk /2. The corresponding unperturbed classical Hamiltonian does not depend on x if expressed in terms of the absolute value of the momentum ζ. Here we introduced covariant and contravariant metric tensors g ik and g ik . The curvature (also referred to as the Gaussian curvature) of configuration space M 2 is expressed in terms of derivatives of the metric tensor. Since unperturbed classical dynamics conserves energy, the value of ζ is preserved, and evolution actually occurs in the reduced phase space (energy shell) M 3 . If we are not interested in a trivial case of zero energy, distributions can be thought of as functions ρ(x, ζ), with x ∈ M 3 and ζ > 0 is the absolute value of the particle momentum. A compact 3-dimensional smooth manifold M 3 is the subspace of the phase space that consists of points with unit length of the momentum vector. Points x ∈ M 3 are specified by two coordinates r and momentum direction angle θ. The dynamical symmetry for the free-particle dynamics originates from the fact that there is a smooth action of the group G ∼ = SO(2, 1) in the reduced phase space of the system, consistent with the dynamics. The details are presented in Appendix A. The dynamical symmetry has a simple and clear interpretation in infinitesimal terms, i.e. action of the corresponding Lie algebra so(2, 1), where the algebra generators are implemented as vector fields in M 3 .
As discussed in Section II, the evolution in the phase space is determined by equationη =Lη. In this equa-tionLη is a vector field that belongs to the tangent space of phase space point η. Hereafter we adopt an agreement used in differential geometry by identifying a first-order differential operator of differentiating along the vector field with the vector field itself. Due to energy conservation, vector fieldL is tangent to the reduced phase space M 3 . We introduce a vector field σ 1 that determines the geodesic flow and generates the translation along trajectories, so that the Liouville operator of free motion iŝ Another natural vector field in M 3 is σ z which corresponds to the generator of the momentum rotation and keeps the coordinates unchanged. This operator can be represented as σ z = ∂/∂θ in terms of a derivative with respect to the particle momentum direction θ.
Operators σ 1 and σ z can be easily expressed in terms of r and p. Performing the variables transformation we obtain: Here E ij is the asymmetric tensor, E 12 = −E 21 = (det g ik ) −1/2 , E 11 = E 21 = 0. Both operators σ 1 and σ z obviously conserve energy: σ 1 ζ = σ z ζ = 0. Finally, we introduce the third vector field in M 3 as σ 2 = [σ 1 , σ z ] where the commutator of vector fields is understood as the commutator of the corresponding differential operators. A simple calculation shows that in the case of constant negative curvature the vector fields σ z , σ 1 , and σ 2 form the Lie algebra so(2, 1) with respect to the vector field commutator (see Appendix B for more details): The group action of SO(2, 1) in the reduced phase space M 3 is obtained by integration of the algebra action.
Any compact surface of constant negative curvature may be represented as a result of factorization of the hyperbolic plane with respect to translations that constitute a finitely generated infinite discrete subgroup of SO(2, 1) (see Appendix A). The entire hyperbolic plane possesses the complete SO(2, 1) symmetry, which makes the infinite motion there integrable due to the presence of the integral of motion similar to the angular momentum in the case of SO(3). The factorization destroys global SO(2, 1) symmetry and folds back the trajectories, which renders the dynamics strongly chaotic.
The Poisson bracket can be computed in a standard way, making use of the fact that locally it coincides with the Poisson bracket for a free particle moving in the entire hyperbolic plane H (see Appendix B). The canonical Poisson bracket may be rewritten in terms of the linear differential operators introduced above acting on two functions f (x, ζ) and g(x, ζ) as The Poisson bracket in Eq. (21) is expressed as a bilinear form of generators of a Lie algebra with constant coefficients. This is is another important manifestation of the dynamical symmetry in the problem of free motion on the surface of constant curvature. The form of the Poisson bracket as well as the commutation relations suggest that the Lyapunov exponent is equal to λ = √ −K. This can be clearly seen from the Jacobi equationÿ + Ky = 0 for the magnitude of the normal component of the deviation from a given geodesic [13]. The stable and unstable directions are given by linear combinations of σ 2 and σ z . The displacements along the flow σ 1 are conserved. The deviation of the trajectory at t = 0 can be written as δη(0) = σ(η), where σ is a first order differential operator corresponding to a certain direction in the phase space. The evolution of the phase point η+δη(0) is given by e σ1t (η + ση) = η(t) + e σ1t σe −σ1t η(t). We see that after time t the deviation from the phase point η(t) is determined by the vector field σ(t) = e σ1t σe −σ1t . Commutation relations (20) allow to find this Heisenberg representation of any operator decomposed over the basis elements of so(2, 1). In particular, we find that and thus conclude that the local stable and unstable directions are given by σ 2 − σ z and σ 2 + σ z , respectively. The fact that the form of these vector fields is conserved by the dynamics is an important manifestation of the dynamical symmetry. Dynamical symmetry also leads to the following general relation: for l = 1, 2 or z that is valid for any smooth function f (x) where dx in M 3 is an invariant integration measure with respect to the SO(2, 1) action (see Appendix A for details). Eq. (23) implies an integration-by-part rule that will be an important ingredient of our analytical calculations.
In conclusion we emphasize that the dynamical symmetry with respect to the action of the group G ∼ = SO(2, 1) does not mean symmetry in a usual sense, i.e. that the system dynamics commutes with the group action, but rather reflects the fact that the vector fieldL that determines the classical dynamics is represented by an element of the corresponding Lie algebra so (2,1). This allows to apply representation theory as described below.
B. Decomposition of the free-particle dynamics using irreducible representations The smooth action of G in M 3 can be interpreted as that the space H of smooth functions in M 3 constitutes a representation of G, which turns out to be a unitary representation (see Refs. 23, 24, 25 and Appendix A), and therefore can be decomposed into a direct sum of irreducible representations of G. The spectrum Spec 0 (M 2 ) ⊂Ĝ of the Riemann surface is defined as a discrete subset of the spaceĜ of irreducible unitary representations of the principal series that participate in the decomposition of functions in M 3 into irreducible representations: The spectrum Spec 0 (M 2 ) consists of imaginary numbers s that characterize representations of the principal series. The effects of the inclusion of other (complementary and discrete) series are discussed after their review in Appendices A and C. The evolution in the space of distributions is decomposed into a set of uncoupled evolutions that correspond to relevant irreducible representations. The component dynamics is determined by the reduced Liouville operators: whereL(s) corresponds to an irreducible representation labeled by s. The distribution components can be further decomposed as using the eigenstates of the angular momentum operator σ z . These satisfy the following properties: where we introduced the raising and lowering operators σ ± = σ 1 ± iσ 2 that are anti-Hermitian conjugated, i.e.
The description of unitary representations of G, reviewed in detail in Appendix C (see also Refs. 24,25), allows to represent the eigenstates ψ k (x; s) by the functions on the circle Ψ k (u) for any given s ∈ Spec(M 2 ). In the case of imaginary s (principal series), when Ψ m (u) constitutes an orthogonal normalized set with the natural scalar product, this leads to normalized functions ψ k (x; s). Identification of the functions ψ k (x) with Ψ k (u) establishes an isomorphism between two irreducible representations, the first being H s that participates in the decomposition of Eq. (A4), the second being its standard representation in function in a circle described in Appendix C. According to the Shur lemma (see, e.g. [23]) the identification (isomorphism) is determined up to a factor. Its absolute value can be fixed by requiring that the function ψ 0 (x; s), identified with Ψ 0 (u) = 1 is normalized. To fix also its phase we note that ψ 0 (x; s) does not depend on θ and turns out to be an eigenfunction of the Laplacian operator ∇ 2 in M 2 as described by Eq. (C5). Since the Laplacian is a real operator, the eigenfunction ψ 0 (x; s) can be chosen to be real. Hereafter, we implement an agreement that ψ 0 (x; s) identified with Ψ 0 (u) = 1 is real. This determines the functions up to a sign, the latter being of no importance: The form of σ 1 in the angular representation is fixed up to a phase [24]: Thus, the original dynamics of a distribution g s (x; s) in Eq. (27) is mapped onto an effective classical dynamics of a distribution G(u) = ∞ k=−∞ ρ s,k (ζ)Ψ k (u; s) defined on a circle with the Liouville operator (29). The resulting effective problem is one-dimensional and can be solved exactly.
The term that describe the interaction with the driving field can be also decomposed in irreducible representations. The coupling (polarization) f is represented by a function of the particle position only (a function in M 2 ) or, stated differently, a phase-space function that does not depend on the particle momentum. It can be equivalently interpreted as a function in M 3 independent of θ, hence σ z f = 0. Thus, f can be expanded as The sum in Eq. (30) runs over the spectrum of the compact surface corresponding to the principal series representations characterized by imaginary s (see Appendix C). We emphasize that, as noted earlier, the functions ψ 0 (x; s), being actually functions of the particle position only, are the eigenfunctions of the Laplacian operator ∇ 2 on the Riemann surface. Typically the dipole f is a slow function of coordinates, and therefore only few terms provide substantial contributions to the expansion of Eq. (30). We also note that, although the Laplacian di-agonalization on an arbitrary Riemann surface with constant curvature is a complex problem, its eigenfunctions are tangible and intuitively simple objects. Combined with the previous results, the expansion (30) leads to closed analytical expressions for the response functions of the original problem.

C. Effective dynamics on the circle
We have identified the dynamical symmetry that allows to map the original chaotic dynamics on M 2 onto a tractable dynamical problem on the circle. In this subsection we consider reduced dynamics that corresponds to a principal series representation labeled by s with Ims > 0. In the following calculations of the response functions, we will need the expansion of e −Lt ψ 0 (x; s) over basis vectors ψ k (x; s), In particular, the linear response is expressed in terms of the coefficient A 0 . These coefficients that are actually matrix elements of the evolution operator between the angular harmonics have been calculated earlier [18,22]. The result immediately follows from the description of irreducible representations using a construction of an induced representation (see, e.g. Ref. 24). The derivation presented below allows an extension to the Langevin dynamics, associated with the original classical dynamical problem. The picture based on the Fokker-Planck equation will be presented in detail elsewhere [28]. The coefficients A k (t; s) are calculated by implementing the representation on the circle, introduced above. Since e −L(s)t Ψ 0 (u) represents a function obtained as the result of the evolution operator action on Ψ 0 (u) ≡ 1, it can be found by solving the dynamical equation ∂ t g(t, u) +Lg(t, u) = 0. For the principal series representation labeled by s, the equation adopts the following form: The solution of the first order partial differential equation supplemented with an initial condition g(0, u) = 1 can be found using the method of characteristics: At long times t > 0 the solution is concentrated near the stable stationary point u = π. This reflects the collapse of the reduced phase space distribution function along the stable direction. Indeed, the width of the region where e −L(s)t Ψ 0 (u) is not exponentially small vanishes as ∝ e −t , according to the fact that for the chaotic dynamics on M 2 with Gaussian curvature K the Lyapunov exponent is equal to √ −K. The expansion coefficients in Eq. (31) are given by The integral can be reduced to a standard integral representation of the Gauss hypergeometric function Γ(z) being the gamma function [26]. An exact expression for k ≥ 0 reads: Various representations of these coefficients using other special functions have been derived in the context of twopoint correlations [18,22]. Since solution g is even, g(t, −u) = g(t, u), coefficients A k are symmetric, A −k = A k . The time reversal property g(−t, u) = g(t, u + π) .
(37) that follows directly from Eq. (32), immediately implies A k (−t) = (−1) k A k (t). We can represent A k (t; s) in an alternative form suitable for studying its long-time behavior [26]: Its expansion in e −2t using Eq. (35) corresponds to the spectral decomposition of the regularized evolution operator for the irreducible representation labeled by s. The properly interpreted eigenvalues of the Liouville operator are known as Ruelle-Pollicott resonances for chaotic systems. They can be obtained using a physical approach. Adding weak Langevin noise results in adding an infinitesimal second-order differential diffusion operator to the classical Liouville operator. This results in a Fokker-Planck operatorL = −κ∇ 2 +L that represents a regularized version of the Liouville operator [27].
The resonances are more often treated as mathematical objects in a form of generalized eigenfunctions in a properly chosen rigged Hilbert space (see e.g. Ref. 18). The physical approach to spectral decomposition of linear and nonlinear response functions for a free particle moving on a surface of constant negative curvature has been developed in Ref. 28. We will be interested in long-time behavior of the response functions and expect that the asymptotics at t → ∞ are relevant. However, we will see that, since the approximation breaks down for higher harmonics k e t and does not vanish at k → ∞, it cannot be immediately utilized for nonlinear response calculations where the result is given by the infinite series over k. In the region t 1 and 1 ≪ k ≪ e 2t one can use an approximation for A k (t) in terms of the MacDonald function (modified Bessel function of the second kind):

V. CALCULATION OF RESPONSE FUNCTIONS
In this section we substantiate the qualitative picture of response presented in Section III which is based on general expressions (10, 11) for a particular model of free motion on a compact surface with constant negative curvature. Heretofore we have not made any assumptions about the shape of the initial distribution ρ 0 in the general expressions for response functions. In the case of free motion, the equilibrium phase space density ρ 0 depends on the momentum absolute value only. The fluctuation-dissipation theorem (FDT) relates the linear response function to the correlation function in a system at thermal equilibrium. To reduce bulky calculations, in both linear and second-order response functions, we employ an analog of FDT and notice that the action of the evolution operator on the most inner Poisson bracket can be reduced to the time derivative. It allows recasting Eqs. (10,11) in a form: The general expressions given by Eqs. (41) and (42) apply to any equilibrium distribution, including the canonical ρ 0 ∝ e −βE , microcanonical ρ 0 ∝ δ(E − E 0 ), and cumulative microcanonical ∂ E ρ 0 ∝ δ(E − E 0 ) distributions. To obtain the final expressions in the most transparent form we will use the cumulative microcanonical distribution with the cut-off energy set to E 0 = 1/2.
We are able to approach an apparently intractable problem of calculating the response in a chaotic system because of strong dynamical symmetry present.
The preliminary steps of our calculation have been described in the preceding section. Dynamical symmetry leads to decomposition of the original dynamical problem into a set of uncoupled evolutions that correspond to irreducible representations of the dynamical symmetry group SO(2, 1). The dipole moment f is expanded over momentum-independent basis functions ψ 0 (x; s) according to Eq. (30), where s labels relevant irreducible representations. The relevant matrix elements A k (t; s) of the evolution (Perron-Frobenius) operator are given by Eq. (36). The coefficients A k (t; s) constitute the dynamical part in the calculation of the response functions. The geometrical part of the calculation is represented by the integrals over the reduced phase space, hereafter referred to as (geometrical) matrix elements. The geometrical part turns out to be trivial in the linear response case while the calculations for the second-order response are more involved. Computation of the matrix elements involved in the nonlinear response is one of the main technical results of this manuscript. The secondorder response function will be obtained in a form of an infinite series. A nontrivial character of the series convergence reflects a nontrivial way how the exponentially growing terms cancel out in the nonlinear response, as described in Section III B.
In this section for the sake of simplicity the calculations are performed in the case when the dipole moment f has a single component f = ψ 0 (x; s) in the expansion over irreducible representations (30). Stated differently, the dipole function f is represented by a single eigenmode of the Laplace operator. Some details for the case of two components, f = B s1 ψ 0 (x; s 1 ) + B s2 ψ 0 (x; s 2 ) (superposition of two Laplacian eigenmodes), are presented in Appendix D. Numerical results for the second-order response function in this simplest case that involves different resonances are given in Section VI.

A. Linear response
Since the basis functions ψ k (x; s) with different k are mutually orthogonal, the linear response function is determined by the dynamical part alone and can be represented in a form The second-order response function for a free particle on a Riemann surface with constant negative curvature is obtained by substituting the specific expression for the Poisson bracket [Eq. (21)] into the generic expression given by Eq. (42). We further simplify the calculation by making use of the evolution operator unitarity e −Lt2 † = eL t2 with respect to the natural scalar defined by Eq. (C2). Performing the integration over ζ by parts and making use of ∂ρ 0 /∂ζ| ζ=0,∞ = 0 results in: The transformation to Eq. (45) is an important step towards expressing the result in terms of the quantities calculated in the previous sections. This allows to avoid propagating the vector fields σ 2 and σ z , and instead deal only with the action of the evolution operators eL t2 and e −Lt1 on the dipole moment f . The latter is decomposed according to Eq. (30) in zero harmonics ψ 0 (x; s) of irreducible representations labeled by s, i.e. the Laplacian eigenmodes. Then we use the expansion (31) to express e −Lt ψ 0 (x; s) in terms of the angular harmonics ψ k (x; s) which constitutes a basis set in the representation. Now we approach the most challenging part of the calculation represented by the integration over the reduced phase space. Generally, according to the decomposition of f , the integration may involve functions from different irreducible representations. In this section we consider the simplest case f = ψ 0 (x; s) of the dipole function represented by a single eigenmode of the Laplacian. A more complicated case of two representations is treated in Appendix D.
One needs to find a way to calculate the integrals of triple products like ψ * n (x; s)(σ 1,2 ψ 0 (x; s))ψ m (x; s) over the reduced phase space M 3 . This matrix element in the second-order response function is essentially less trivial part than its counterpart in the linear response function, where the integral of the double product of ψ n (x; s) has been done based on orthogonality and normalization of the functions. The problem arises from the fact that integrals of relevant triple products over the reduced phase space do not have natural representation in terms of the effective system on the circle. Nevertheless, we will show that this geometrical part of the second-order response can be worked out using relatively simple tools, and can be expressed in terms of few quantities related to the Laplacian eigenmodes.
The integration in reduced phase space x = (r, θ) is performed over the particle coordinates r ∈ M 2 and the locally defined angle θ that represents the 2D momentum direction. The local (for fixed r) dependence of ψ n (x; s) on θ is given by ψ n (x; s) ∝ e inθ due to the definition ∂ θ ψ n (x; s) = inψ n (x; s) of the basis functions. The operators σ 1 and σ 2 applied to ψ 0 (x; s) in the middle of the integrand are expressed via the ladder operators σ ± . Therefore, integration over θ in Eq. (45) results in zero unless the integrand is locally independent of θ, i.e.
We now introduce the matrix element a k and b k as As established in Section IV, the differential operators σ ± can be moved from one part of the integrand to another one, and thus one can integrate by parts [Eq. (24)]. Integrating by parts in Eq. (48) and making use of the properties (28) of σ operators we derive the following implicit recurrence relations: We further make use of Eq. (49) to express a k+1 via a k and b k+1 , and applies Eq. (50) to find b k+1 in terms of a k and b k . This allows to recast the recurrence relations between the matrix elements in an explicit form: The matrix elements b k can be excluded to connect three adjacent matrix elements a k−1 , a k and a k+1 : As mentioned earlier ψ 0 can be chosen real, which leads to the relations a 1 = a −1 = a 0 /2 (see Appendix D for the details). This allows to solve the recurrence relations and express all matrix elements in terms of a single real number a 0 . In addition, in Appendix D we consider a more complicated case of matrix elements composed of functions that involve different representations.
Putting all contributions to the second-order response function together, we can write it down in the following form: Here we used the short notation A n (t) for A n (t; s) and the fact that the product (2n+1−2s)A * n+1 (t 2 ; s)A n (t 1 ; s) is real. We used the symmetries A −n (t) = A n (t) and a −n = a n discussed in Section IV C and Appendix D to restrict the consideration to non-negative n.
We note that the coefficients A * n (−t 2 ; s) at negative time −t 2 < 0 enter the expansion for the backward evolution of ψ 0 (x; s), which appears to be crucial for the calculation below. The time reversal symmetry A n (−t; s) = (−1) −n A n (t; s) was also used in Eq. (53).
The result for the second-order response function (53) is written in the form of an ordinary series of terms whose expressions are either completely known analytically (see Eq. (36)) or recurrently expressed via a single number a 0 that characterizes function ψ 0 (x; s). The functions ψ 0 do not depend on the angle θ, i.e. ψ 0 (x) = φ(r; λ), where, according to Eq. (C6) φ(r; λ) is an eigenmode of the Laplace operator ∇ 2 with the eigenvalue λ on our Riemann surface. Therefore, Eqs. (A3) and (C6) allow to express the only undetermined parameters s and a 0 in Eq. (53) in terms of the relevant eigenmode φ(r; λ) of the Laplacian: Finding the eigenmodes of the Laplacian analytically is a complex problem. Fortunately this is not necessary for our purposes. It is known that for a given constant negative curvature K < 0 there is a family of non-equivalent Riemann surfaces of genus g that can be parametrized by a space M g with the dimension dimM g = 4g − 2 known as the moduli space. In particular, a 0 and s can be viewed as some nontrivial functions on the moduli spaces M g , and computation of these functions is a complex problem. Yet, we can note that for any given values of a 0 and s we can find a Riemann surface that implements them according to Eq. (54). Therefore, we can treat a 0 and s as free parameters that characterize some Riemann surface of constant negative curvature.
In a more general case, when the dipole f is represented by a finite superposition of N f Laplacian eigenmodes, the second-order response can be expressed in terms of a larger set of parameters s j and a sisj s k 0 , with i, j, k = 1, . . . , N f as described in Appendix D. We can apply the same argument to treat them as free parameters for a Riemann surface of a high enough genus g.

C. Second-order response: summation of the series
We have expressed the second-order response in terms of an ordinary converging series for the cases of a single Laplacian mode N f = 1 [Eq. (53)] and two modes N f = 2 [Eq. (D7)], respectively, in the expansion of the dipole (30). A general expression for a finite number N f of modes has a similar structure. Each term of the series is known analytically in the form that involves special functions and solutions of recurrence relations. Or goals are to derive the long-time asymptotics of the secondorder response S (2) , and develop a computationally efficient procedure for S (2) (t 1 , t 2 ) at finite times.
We start with the asymptotic behavior of the secondorder response function for large t 1 and t 2 . Naively, one would plug the long-time asymptotic of A n (t; s) from Eq. (39) into Eq. (53). However, this asymptotic of A n (t; s) is valid only for fixed n, if ne −t ≪ 1. This can be confirmed by comparing the consecutive terms in the expansion of the hypergeometric function in Eq. (36) in powers of e −t . Moreover, although A n (t; s) represent the Fourier expansion coefficients of a smooth function (33), their asymptotic (39) does seem to vanish as n → ∞. Counting the powers of n in the second part of the summand in Eq. (53) estimated using these asymptotic expressions reveals that the resulting series fails to converge.
Indeed, the series (53) can be represented in the form where the dependence of F (n; t 1 , t 2 ) on n is slow for n ≫ 1. We may find the long-time behavior of the terms as The double expansion originates from the decomposition of A n (t; s) in powers of e −2t . The first term in the latter decomposition, specified by Eq. (39), turns out to be independent of n in the limit of large n, apart from irrelevant slow quasi-oscillations ∝ n s that cannot affect the series convergence. A naive long-time asymptotic of the series is given by n (−1) n F 00 (n; t 1 , t 2 ), where the n-th term is estimated as F 00 (n; t 1 , t 2 ) = n 2 (a n − a n+1 )e − t 1 +t 2 where quasi-oscillatory dependence of a n , C n (s) and D n (s) on n ≫ 1 does not affect the convergence. Taking into account the asymptotic form of a n [see Eq. (E1)], we conclude that the series is obviously divergent as In what follows we show how to overcome the apparent divergence and calculate the asymptotic of the second-order response function analytically.
The series in Eq. (55) that represents the response function must converge for fixed values of t 1 and t 2 . Indeed, the terms F (n; t 1 , t 2 ) ∝ exp[−2n(e −t1 + e −t2 )] vanish exponentially, although only for extremely large n ≫ e t1 , e t2 ≫ 1. The decay rate is determined by the fact that g(t, u) ≡ e −Lt Ψ 0 (u) in Eq. (33) is smooth on scales that do not exceed the smallest fettuccine size e −t . The ultimate convergence allows to safely regroup the terms of the series: After the terms are regrouped the initial naive approach based on approximating F (n; t 1 , t 2 ) by its longtime asymptotic F 00 (n; t 1 , t 2 ) results in a converging series, since the linear combination in the summand represents the discrete counterpart of the second derivative d 2 F 00 (n)/dn 2 which decays as ∝ n −3/2 . According to Eq. (57) all terms in the resulting series have the same time dependence. Therefore, the long-time asymptotic of the response function S (2) (t 1 , t 2 ) is represented by a superposition of four components e −(t1+t2)/2±s(t1±t2) with the coefficients in the form of convergent series. Coefficients F kl (n; t 1 , t 2 ) at higher orders in the expansion (56) grow faster with increasing n as F kl (n; t 1 , t 2 ) ∼ n 2(n+k) F 00 (n; t 1 , t 2 ). The series re-grouping approach of Eq. (58) can be generalized to eliminate the apparent divergences for higher-order terms in Eq. (56), and obtain the higher-order terms in the asymptotic expansion of the second-order response in powers of e −2t1,2 . Instead of that we suggest an alternative procedure that allows (i) to demonstrate the existence of a long-time asymptotic expansion of S (2) in powers of e −2t1 and e −2t2 , (ii) derive relatively simple expressions for the expansion co-efficients in any order, and (iii) develop an efficient numerical scheme for computation of the response at finite times.
We start with deriving the asymptotic expansion. To that end, provided t 1 , t 2 ≫ 1, we introduce an intermediate N ≫ 1, so that for n > N the terms F (n; t 1 , t 2 ) are represented by smooth functions of n. We further partition the sum S of the series into the finite sum S N = n≤N (−1) n F (n) and the remainder R N = n>N (−1) n F (n), followed by evaluating the remainder. Implementing the definition [26] of the Euler polynomials E n (y), we derive the following identity, based on the Taylor expansion of (1 + e −x ) −1 in x = d/dz: which is used to calculate the sum of consecutive terms pairwise, so that the remainder R N = n>N (−1) n F (n) of the almost alternating series becomes related to the first term F (N ) that is not included in R N , and its derivatives: We further introduce the following "improved" partial sums: The reason for introducing the improved partial sums P can be viewed as a set of counter-terms that eliminate the divergence in the first term. In the limit N → ∞ Eq. (61) leads to the expansion The expansion of Eq. (63) is represented by an asymptotic rather than converging series. It can be viewed as a double expansion of the second-order response function S (2) (t 1 , t 2 ) in Ruelle-Pollicott resonances which originates from the uncoupled evolution of the chaotic system during time intervals t 1 and t 2 . This will be demonstrated explicitly [28] by analyzing the noiseless limit of the corresponding Langevin dynamics.
At this point we should note that our derivation of the asymptotic expansion has been somewhat frivolous. First of all, since the terms F (n; t 1 , t 2 ) of our original series involve the coefficients a m , we do not actually know a function F (z; t 1 , t 2 ) that represents the series, only its values for positive integer integer z = n are available. In particular the derivatives ∂ m N F m1m2 (N ) in Eq. (62) have been not defined yet. Second, in obtaining the asymptotic expansion we were not controlling the neglected terms. This is especially dangerous in our case since we derive the complete asymptotic, which involves computing the terms that are exponentially small compared to more senior terms in the expansion. We need to make sure that the terms that are neglected in computing a certain expansion coefficient are small also compared to the higher terms that are kept in the asymptotic expansion.
A sketch of the appropriate derivation is presented in Appendix E. In particular, we demonstrate that the coefficients a n that enter the expressions for F (n; t 1 , t 2 ) can be represented by some analytical function a(z), and therefore the terms F (n; t 1 , t 2 ) of the original series are represented by some analytical function F (z; t 1 , t 2 ). We will further demonstrate how to compute derivatives F (m) (N ) without explicit knowledge of the function F (z). We also show how a proper choice of N and n in P (M) N allows to control terms neglected in the asymptotic expansion.
In the remainder of this subsection we will implement the summation procedure in the form of an efficient numerical scheme for computing the second-order response S (2) (t 1 , t 2 ) at arbitrary times. We reiterate that all terms in the series for S (2) are known in their analytical form (via recurrence relations) apart from few parameters determined by the geometry of the surface, as discussed at the end of the previous subsection and in Appendix D.
In the simplest case N f = 1 the single parameter is given by Eq. (54). The series for S (2) is absolutely convergent, which becomes remarkable only at n e t1 or n e t2 , and therefore a straightforward computation involves an exponentially growing with time number of terms with relatively complex structure. At given times t 1 and t 2 the terms F (n; t 1 , t 2 ) of the series constitute a sequence of numbers. Our numerical scheme uses the procedure described above and reduces the problem to the summation of a relatively small and independent of time number of terms. This allows to compute S (2) (t 1 , t 2 ) with minimal numerical effort.
We can justify the numerical procedure separately for two overlapping intervals of t 1 and t 2 . In the case min(t 1 , t 2 ) 1 the series can be truncated already at N ∼ 10. At n 10 the decay of the terms F (n; t 1 , t 2 ) is exponential, and the remainder is negligible. Therefore, the numerical value of S (2) is obtained by the summation of several terms.
In the case min(t 1 , t 2 ) 1 we are not far from the asymptotic region t 1 , t 2 ≫ 1, and S (2) may be found by a simple numerical implementation of the remainder calculation (59). Since we use few terms in the expansion (59) and approximate them by finite differences, the numerical precision is determined by how smoothly the numbers F (n; t 1 , t 2 ) behave as n increases. This is essentially determined by the value of n. The smoothness is only slightly influenced by the values of t 1 , t 2 1, which can be rationalized by considering the principal asymptotic of S (2) at t 1 , t 2 ≫ 1. The series (55) is not purely alternating because of the presence of quasi-oscillations ∝ n ±s = e ±i|s| ln n in F (n; t 1 , t 2 ). To ensure a sufficiently smooth behavior of F (n; t 1 , t 2 ) at n > N one must choose larger N for larger values of |s|.
In practice, a reasonable relative error 10 −3 is achieved if merely N ∼ 10 2 or less terms are retained in the partial sum S N for s ∼ 5i. In the expansion (59) the third term identically vanishes because E 2 (1) = 0. We use two first terms in the expansion, and represent the first derivative by the discrete difference with the third-order accuracy. Therefore, the neglected contributions scale at best as O(F 00 (N ; t 1 , t 2 )e −4 min(t1,t2) ). The efficiency of the procedure is remarkable. We need to calculate only around N ∼ 10 2 terms, which should be compared with a much bigger number e 10 > 10 4 at which the summand F (n; t 1 , t 2 ) starts to decay exponentially if min(t 1 , t 2 ) ∼ 10. At the first sight, summation of 10 2 terms versus 10 4 does not seem to be a crucial improvement, given today's computational capacities. However, the bottleneck here is computing the individual terms that are represented via the hypergeometric functions. The computational effort grows dramatically for the terms F (n; t 1 , t 2 ) where the series starts to converge in the absolute sense. Besides, straightforward summation of alternating series is not free from accuracy issues. Overall, generating a 2D spectroscopic signal, with the developed approach, requires several minutes of CPU time using such a simple and high-level tool as Mathematica. This qualifies for an "almost analytical" calculation of the complete signal.
In conclusion, we emphasize that only a finite number of the series terms is used in the summation procedure. Thus, the asymptotic expansion of the second-order response function turns out to be determined by the expansions of A n (t; s) at t → ∞. The latter may be viewed as spectral decompositions of the evolution operator matrix elements. In fact, the strongly chaotic system we consider is characterized by Ruelle-Pollicott resonances ω ν = ν + 1/2 ± s. Only even resonances with ν = 2k, where k = 0, 1, . . ., contribute to the spectral decomposition. This result for the second-order response is nontrivial because the expansion can be made only after the summation of the series over angular harmonics. This means that the expansion is only asymptotic. The convergence issues can be completely avoided by introducing infinitesimal noise. Nonzero noise provides the convergence of the series over angular harmonics for a given pair of resonances [28]. In the limit of the vanishing diffusion coefficient, the terms of the series turn into their noiseless counterparts introduced in Eq. (56). The application of the remainder summation procedure, described above, requires only smoothness of the series terms. Therefore, the sum of the series is determined by a finite number of terms, and the asymptotic expansion appears independent of the vanishing diffusion coefficient. This noise regularization is meaningful, since the way the convergence is enforced is irrelevant for the asymptotic expansion.

VI. NUMERICAL RESULTS
Experimental data on 2D time-domain spectroscopy that probes the response function S (2) (t 1 , t 2 ) is usually presented using the so-called 2D spectrum which is obtained by a numerical 2D Fourier transform of the response function with respect to t 1 and t 2 : To provide a spectroscopic view of chaotic vibrational dynamics, in this section we present some numerical results on 2D spectra for a particular strongly chaotic system studied in this paper.
As we have shown above, the 2D response in our model can be expressed in terms of the properties of the Laplacian eigenmodes that participate in the expansion of the dipole function f (r). Since the dipole is generally a smooth function of the system coordinates, the number N f of relevant eigenmodes is typically small. To study the important features of the signals we consider the second-order response function in the simplest cases N f = 1 and N f = 2. The latter represents the simplest situation that involves diagonal and cross peaks in the 2D spectrum. Analytical results for this case which generalize Section V are derived in Appendix D.
In the case N f = 1 the second-order response function depends on a spectral parameter s and an additional pa-rameter a 0 . Both parameters are expressed in terms of the only Laplacian eigenmode that represents f (r) [see Eqs. (43), (54) and (C6)].
For N f = 2 we have two spectral parameters s 1 and s 2 , and four additional parameters. Similarly to the N f = 1 case, all parameters can be expressed in terms of the two relevant Laplacian eigenmodes. As explained at the end of Section V C, the parameters can be considered as independent attributes of a particular Riemann surface. For the demonstration of qualitative features, we use a particular choice of the parameters.
The absolute value of the 2D spectroscopic signal |S (2) (ω 1 , ω 2 )| is presented in Fig. 3. Panel (a) corresponds to the N f = 1 case with the only spectral parameter s = 5i. Panel (b) shows the N f = 2 case with two spectral parameters s 1 = 5i and s 2 = 3i. The real and imaginary components of the 2D spectra are presented in Fig. 4. In the first case we see a diagonal peak with a pronounced stretched feature along ω 1 direction. In the second case we also see cross peaks accompanied by similar stretched features.
Diagonal and off-diagonal peaks are also observed in spectroscopic signals for harmonic and almost harmonic vibrational dynamics. In this case the positions of the peaks are given by the frequencies of the underlying periodic motions, whereas the width is determined by the system-bath interactions and has about the same value in both frequency directions. In the chaotic case the peak positions are determined by Ruelle-Pollicott resonances, rather than by the frequencies of some specific periodic orbits. The stretching in spectroscopic signals originates from time-domain damped "breathing" oscillations with a variable period caused by strong nonlinearity of the underlying vibrational dynamics and can be interpreted as a signature of chaos.

VII. CONCLUSION
In the present paper we studied classical response in a strongly chaotic system. Using the Liouville space representation for classical dynamics we found that the response functions exhibit damped oscillations with time. The decay is attributed to the mixing property of the strongly chaotic dynamics which leads to the efficient equilibration in phase space.
Previous studies of the integrable systems have revealed the appearance of unphysical divergences in the classical nonlinear response functions at long times. The divergence originates from the linear time growth of stability matrix elements in integrable systems. Chaotic dynamics is characterized by the exponential growth of certain stability matrix elements.
We have described a general qualitative picture of response in strongly mixing (hyperbolic) systems. The exponential growth of the stability matrices that can potentially lead to exponential divergence in response functions has been interpreted using the Liouville space representation of classical mechanics: As a result of evolution a smooth initial distribution becomes extremely sharp along the stable directions. The second interaction with the driving field involves derivatives of the evolved distribution along the stable directions. This could result in exponentially growing terms in nonlinear response functions. We have demonstrated that due to smoothness of the initial distribution the dangerous terms completely cancel out, and the nonlinear response functions exhibit exponential decay at long times. Stated differently, in mixing systems stability matrices have exponentially growing and exponentially decaying components. Due to the smooth character of the dipole function that describes the coupling to the driving field, the growing components of the stability matrices are simply eliminated from the game, and the physical behavior of the nonlinear response functions is determined by the exponentially decaying components.
To confirm the established qualitative picture we performed detailed calculation of the linear and second-order response for a chaotic model of a free particle moving along a compact surface with constant negative curvature. The model possesses dynamical symmetry that allows for an exact solution by applying the group representation theory. We found the long-time asymptotic behavior of the linear and second-order response function that has a form of exponentially damped oscillations. Complete asymptotic series obtained in this paper can be viewed as expansions in Ruelle-Pollicott resonances. The expansion of the linear response in RP resonances is in agreement with the earlier results for the two-time correlation functions [18], since the latter is directly related to the linear response via the FD theorem.
The analogue of the RP expansion for the nonlinear response is of more nontrivial nature. The eigenmodes that represent the RP resonances are generalized, rather than smooth functions. In the linear case the initial smooth distribution should be decomposed in the RP modes. The signal is computed by convoluting the RP modes with the smooth dipole function. Both operations are welldefined for generalized functions. In the case of nonlinear response the second interaction with the driving field involves acting on a generalized function with a differential operator followed by projecting it onto a generalized function. The legitimacy of the latter operation is less obvious, and is related to aforementioned cancellations of the dangerous terms. The RP decomposition for the linear response is represented by a converging series, whereas the nonlinear response is given by an asymptotic series. Computation of the expansion coefficients in the nonlinear case requires a delicate summation of almost sign alternating series. These are other signatures of the nontrivial character of the RP decomposition in the nonlinear case.
In our forthcoming publication [28] we consider Langevin dynamics associated with the model considered here. Spectral decomposition in the corresponding eigenmodes of the Fokker-Planck operator is a stable legitimate procedure. The spectral decomposition in the presence of noise is a converging series. We show that in the limit of vanishing noise the converging series becomes asymptotic and reproduces the expansion derived in this paper. The asymptotic expansions of the linear and nonlinear response functions can be interpreted as decompositions in RP resonances.
We suggest to apply our results for the interpretation of spectroscopic data. Classical chaos is quite generic for large molecules. Even for smaller number of nuclear degrees of freedom, e.g. in small systems of hydrogen bonds [29], the shape of the effective potential energy can make the dynamics chaotic. We have argued that restricted motion in the regions with negative curvature of the molecular potential can be qualitatively described as free motion along a compact Riemann manifold with negative curvature. Moreover, the chaos may be caused even by the positive but inhomogeneous curvature [30] of the configuration space.
We have considered in detail the second-order response which is absent in the spectroscopy of the bulk materials or in isotropic environments. Spectroscopy on the surface [31], in the absence of the inversion symmetry, can measure nonvanishing response functions of even order.
A dynamical system of the general type has the mixed phase space which includes both chaotic regions and stability islands. Even in such systems, when the true longtime asymptotic decay of correlations is rather powerlike than exponential, RP resonances may by noticeable, leading to the intermediate asymptotic, exponential decay persisting for a very long time [32].
When 2D spectroscopic data is interpreted in terms of the underlying dynamics the peaks are usually attributed to some periodic motions in the system. We have shown that chaotic dynamics can result in similar diagonal and off-diagonal peaks in spectroscopic signals that should be attributed to RP resonances, and are nor related to any specific periodic orbits. Note that the frequencies and decay rates of the RP resonances can be retrieved from the dynamical ζ-function that can be represented as a product over all periodic orbits [17]. Stated differently, the RP resonances are related to periodic orbits, yet in an extremely collective way, and may not be attributed to any particular periodic motions. Therefore, they can be referred to as collective chaotic resonances. Our numerical results show pronounced stretched features associated with diagonal and off-diagonal peaks. These features result from breathing damped oscillations that originate from contributions of multiple RP resonances with similar oscillation frequencies (imaginary parts) and different damping rates (real parts). Breathing oscillations are known to be typical for strongly nonlinear dynamical systems. On the other hand, 2D spectra in harmonic or almost harmonic systems coupled to a multi-mode harmonic bath usually show similar peak patterns in both frequency directions. We suggest that the stretched versus non-stretched peak shape could have a capacity of distinguishing between the collective versus individual nature of resonances in spectroscopic signals. This Appendix contains some basic aspects on the geometry of Riemann surfaces with constant negative curvature that are necessary for employing the SO(2, 1) dynamical symmetry. This allows to decompose the freeparticle dynamics using irreducible representations and map the original dynamical problem onto some effective 1-dimensional dynamics in a circle.
The fundamental group Γ g = π 1 (M 2 g ) of a Riemann surface M 2 g of genus g which describes noncontractible closed paths is generated by 2g elements a j , b j , where j = 1, . . . , g, with the only relation g j=1 a j b j a −1 j b −1 j = 1. In the case g > 1 the compact surface is covered H → M 2 g by the hyperbolic plane that can be implemented as a pseudosphere determined by the equation y 2 1 + y 2 2 − y 2 0 = −1 with y 0 > 0 embedded in the 3D Minkowski space with a metric dℓ 2 = −dy 2 0 + dy 2 1 + dy 2 2 (other well-known equivalent implementations include the Poincare disc and the upper complex half-plane). The fundamental group acts freely in H, so that M 2 g ∼ = Γ g \H. The reduced phase space M 3 g , considered as a bundle S 1 → M 3 → M 2 , whose fibers are unit velocity or momentum vectors, being pulled back to H forms a bundle S 1 → G → H, where G ∼ = SO(2, 1) can be represented as a pseudo-orthogonal group. This can be visualized as follows: G ∼ = SO(2, 1) consists of points in H (positions) that can be considered as 3D vectors with norm −1 in the Minkowski metric, together with unit velocity vectors that can be interpreted as norm 1 vectors orthogonal to the position vectors (with respect to the Minkowski metric). Extending these two pseudo-orthonormal vectors to a pseudo-orthonormal basis set, we interpret G as the space of pseudo-orthonormal basis sets, the latter can be thought of as the pseudoorthogonal group SO(2, 1). Factorizing G ∼ = SO(2, 1) with respect to the right action of its maximal compact subgroup K ∼ = SO(2) we arrive at H ∼ = G/K. The action of Γ g in H can be naturally extended to a left action of Γ g in G, which determines the embedding Γ g ⊂ G. This leads to convenient for our purposes representations M 2 g ∼ = Γ g \G/K and M 3 g ∼ = Γ g \G. In particular this interprets the action of G in M 3 g as originating from right action of G in itself.
The presented picture has a very transparent interpretation. We start with a Riemann surface M 2 g of genus g > 1 with constant negative scalar curvature K = −1.
As stated in section IV we have three canonical vector fields in M 3 g , i.e. σ z = ∂/∂θ, σ 1 that determines the geodesic flow (classical dynamics of a free particle), and σ 2 = [σ 1 , σ z ]. In the case of constant curvature they form the Lie algebra so(2, 1) with respect to the vector field commutator. Explicit expressions for the vector fields are presented in Section IV and in Appendix B. This defines an action of so(2, 1) in M 3 g that can be natu-rally extended to G → H, considered as the pull-back of M 3 g → M 2 g to H. The algebra action in G can be integrated to a group action, which implies that locally G has a structure of the universal Lie group associated with the algebra so(2, 1). Careful studying of the global properties of G shows that in fact G ∼ = SO(2, 1). This immediately implies that H → M 2 g is equivalent to SO(2, 1)/SO(2) as a Riemann space, i.e. H is the hyperbolic space. This has a very important implication that the local structure of any Riemann surface M 2 g with constant scalar curvature K = −1 is locally equivalent to the hyperbolic space H, whereas M 3 g is locally equivalent to G. In particular, this implies that all local quantities such as the Laplacian operator ∇ 2 in M 2 g , the Poisson bracket ω and the geodesic flow in M 3 g can be actually computed in H and G, followed by expressing them using the canonical vector fields σ l , l = 1, 2, z. This is how, e.g. Eq. (21) can be derived in a straightforward way.
Finally we note that for a Riemann surface with g > 1 and not necessarily constant curvature the metric determines a complex-analytical structure, and the universal cover H → M 2 g is equivalent to the hyperbolic plane and preserves the complex-analytical structure. This implies M 2 g ∼ = Γ g \H. Since the group of conformal diffeomorphisms of H coincides with its isometry group SO(2, 1), we have Γ g ⊂ G. This defines a metric in M 2 g that has constant curvature K = −1 and is conformallyequivalent to the original one. Since that, denoting by M g the moduli space of complex analytical structures for genus g we can describe a Riemann surface with constant curvature by a pair (η, K) with η ∈ M g and K < 0 being the scalar curvature.
The representations of M 3 and M 2 in terms of the groups K, Γ ⊂ G M 3 ∼ = Γ\G , M 2 ∼ = Γ\G/K ∼ = Γ\H ∼ = M 3 /K , (A1) H ∼ = G/K are very convenient for describing invariant integration measures. Since G, K, and Γ are unimodular groups, according to the theorem on invariant measures in homogeneous spaces [23,24,25], there is a unique invariant measure in Γ\G that is locally equivalent (in the sense of covering) to the Haar invariant measure in G. The latter generates an invariant measure in G/K, which is locally equivalent to the measure in M 2 generated by the metric. Combined with the main property of invariant measures in homogeneous spaces, this implies for any integrable function g(x) This measure provides an invariant scalar product in the space H of functions in M 3 . Besides, implementing functionsg(r) in M 2 as functions g(x) in M 3 that do not depend on θ, i.e. σ z g = 0, we have and, therefore H can be decomposed into a direct sum of irreducible representations of G In Eq. (A4) H (0) denotes the one-dimensional unit representation that represents constant functions, H s denote the principal series representations (with purely imaginary values of s ∈ Spec 0 (M 2 ) where Ims > 0) and complimentary series representations (with real values of 0 ≤ s < 1/2), whereas H (n) are the discrete series representations with the integer factors m n ≥ 0 describing how many times a representation participates in the decomposition.
APPENDIX B: ALGEBRA so(2, 1) AND POISSON BRACKET In this Appendix we derive the Poisson bracket in the form of Eq. (21) with the commutation relations between the operators σ l (l = 1, 2, z) given by Eq. (20).
We start with the canonical form of the Poisson bracket in terms of coordinates r i and conjugated momenta p i = ∂L/∂ṙ i : where we imply summation over repeated indices. Since free motion conserves the absolute value of the momentum ζ, it is convenient to make another local choice of variables in the tangent bundle T M 3 . We represent the particle momentum in terms of ζ and θ, the polar angle which determines the momentum direction. The algebra element σ z which generates momentum rotations in T M 3 reads σ z = ∂/∂θ.
The third vector field in the reduced phase space is described by the differential operator Next, we show that σ 1 , σ 2 , σ z form the algebra so(2, 1). The coefficients of the first order linear differential operator [σ 2 , σ z ] include the metric tensor and its first derivatives. A straightforward yet tedious calculation performed in the basis set of four vectors related to the canonical variables p i , q i results in the relation [σ 2 , σ z ] = −σ 1 . A similar calculation yields the commutator [σ 1 , σ 2 ] = −K(r)σ z which contains configuration space curvature K(r) and is indeed expressed only via the remaining operator σ z . In the case of the constant negative curvature, rescaling σ 1 and σ 2 leads to the so(2, 1) commutation relations: [σ 1 , σ 2 ] = σ z , [σ 1 , σ z ] = σ 2 , [σ 2 , σ z ] = −σ 1 .
A number of transformations allow to express the Poisson bracket (B1) in terms of the operators ∂ ζ , σ 1 , σ 2 , σ z . The simple form reflects the presence of the dynamical symmetry. The first two terms of the Poisson bracket determine the phase space velocity. The constant coefficient in front of the second term can be alternatively deduced from the condition S (2) (t 1 , 0) = 0. In this Appendix we describe a convenient implementation of unitary irreducible representation of SO(2, 1) in terms of functions in a circle S 1 . This provides an aforementioned mapping of the free-particle dynamics onto a 1-dimensional problem. Since the group G ∼ = SO(2, 1) has a double covering SL 2 (R) → SO(2, 1), irreducible representations of G are provided by even irreducible representations of SL 2 (R), the latter being well-known (see e.g., [24]). This situation is conceptually close the the case of irreducible representations of SO(3) given by even, i.e. integer-spin, representations of a double cover SU (2) → SO(3) of SO(3). We follow the approach of [24] for SL 2 (R) and translate it to the language convenient for G ∼ = SO(2, 1).
The space H s of an irreducible representation of G that belongs to the principal or complimentary series has a convenient basis set Ψ k of angular harmonics, so that σ z Ψ k = ikΨ k . The anti-Hermitian conjugated raising and lowering operators σ ± = σ 1 ±iσ 2 , Using the allow for the the following relations between the basis functions: In the case of imaginary s (principal series) Ψ m (u) are normalized, and correspond to normalized functions ψ k (x; s), i.e.
where dx is the invariant measure in M 3 defined up to a constant [23,24,25]. We can alternatively view Eq. (C1) as a definition of a so(2, 1) representation parametrized by s. In Appendix B we directly verify that the operators σ ± , σ z satisfy the so(2, 1) commutation relations given by Eq. (20). The fourth group of matrix elements includes a s2s2s1 k , a s2s1s2 k and a s1s2s2 k which obey the relations similar to those for the third group but with s 1 and s 2 interchanged. Thus they expressed in terms of a single real number a s2s2s1 0 . We assume that four numbers a s1s1s1 0 , a s2s2s2 0 , a s1s1s2 0 and a s1s2s2 0 to be independent. This is rationalized in a detailed discussion at the end of subsection V B.

APPENDIX E: ASYMPTOTIC EXPANSION
In this appendix we present some details involved in a derivation of the asymptotic expansion given by Eq. (63).
First of all the derivation is based on the assumption that the the terms F (n) in our series can be represented by some smooth function F (z). This is true for the coefficients A n that enter Eq. (53) due to their explicit form given by Eq. (38). However, the coefficients a n that are also a part of Eq. (53) are given as a discrete set according to Eq. (52). To show that a k can be represented by a smooth function a(z) we note that since a k is a solution of a recurrence relation, its k → ∞ behavior can be readily analyzed, which results in an expansion a k = k − So far we have demonstrated that the terms F (k) in our series are represented by an analytical function F (z), and therefore an approach for calculating the asymptotic expansion introduced in Section V C is in principle legitimate. However, this is not enough yet, since the asymptotic expansion coefficients S (2) m1m2 involve the derivatives ∂ m z F m1m2 (z) at integer z = k values of the argument, as prescribed by Eq. (62). Apparently, it requires explicit knowledge on the analytical function F (z). However, what we really need as a good enough approximation for the derivatives in Eq. (62) so that the neglected terms vanish in the limit N → ∞. This can be achieved by implementing a very simple scheme. For an analytical function F (z) denote by F  In the remainder of this appendix we show how to control the neglected terms in a derivation of the asymptotic series. We start with noting that the asymptotic expansion of Eq. (63) should be understood in a way that for any pair (M 1 , M 2 ) of positive integers we have Full control over the neglected terms can be achieved by a proper choice of the intermediate N (t 1 , t 2 ) as a function of t 1 , t 2 → ∞ and n in Eq. (61). Note that both N (t 1 , t 2 ) and n also depend on M 1 and M 2 . We choose N (t 1 , t 2 ) = exp (γ min(t 1 , t 2 )) , (E5) so that for small enough γ we have 1 ≪ N ≪ e t1 , e t2 and we can expand F (n; t 1 , t 2 ) for n ≤ N in powers of e −2t1 and e −2t2 according to Eq. (56). If we keep the terms up to the degree (M 1 , M 2 ), the neglected terms can be estimated as ∼ N 2(M1+M2)+r exp (−2 min(t 1 , t 2 )) e −2(M1t1+M2t2) , where r = 3 corresponds to a very loose estimate. By choosing γ to be small enough, specifically (2(M 1 + M 2 ) + r) γ < 2, the neglected terms are o e −2(M1t1+M2t2) and can be safely omitted. Finally, no matter how small γ is we can always find M to be large enough (and independent of t 1 and t 2 ), so that N −M (t 1 , t 2 ) = o e −2(M1t1+M2t2) , and the terms omitted in Eq. (62) can be also safely neglected.