Energy returns in global AdS_4

Recent studies of the weakly nonlinear dynamics of probe fields in global AdS$_4$ (and of the nonrelativistic limit of AdS resulting in the Gross-Pitaevskii equation) have revealed a number of cases with exact dynamical returns for two-mode initial data. In this paper, we address the question whether similar exact returns are present in the weakly nonlinear dynamics of gravitationally backreacting perturbations in global AdS$_4$. In the literature, approximate returns were first pointed out numerically and with limited precision. We first provide a thorough numerical analysis and discover returns that are so accurate that it would be tantalizing to sign off the small imperfections as an artifact of numerics. To clarify the situation, we introduce a systematic analytic approach by focusing on solutions with spectra localized around one of the two lowest modes. This allows us to demonstrate that in the gravitational case the returns are not exact. Furthermore, our analysis predicts and explains specific integer numbers of direct-reverse cascade sequences that result in particularly accurate energy returns (elaborate hierarchies of more and less precise returns arise if one waits for appropriate longer multiple periods in this manner). In addition, we explain, at least in this regime, the ubiquitous appearance of direct-reverse cascades in the weakly nonlinear dynamics of AdS-like systems.


Introduction
Investigations of weakly nonlinear dynamics of global Anti-de Sitter (AdS) spacetime, initiated by the numerical evidence for its nonlinear instability presented in [1], have resulted in a number of surprises. This includes evidence for turbulent cascades leading to black hole formation starting from small initial data [1], as well as non-collapsing initial configurations [2][3][4][5][6][7], and a number of intriguing results within the resonant approximation of the AdS dynamics accurately discribing slow weakly nonlinear transfer of energy between the linearized modes [8][9][10]: selection rules prohibiting certain interaction channels between the linearized modes [9][10][11][12], extra conservation laws [10,13,14], dual energy cascades [14], and strong numerical evidence for turbulence within the resonant system [15] (for a review, see [16]). While the bulk of these considerations have focused on the case of spherically symmetric perturbations of gravity-scalar field systems, extensions to pure gravity within squashed sphere ansätze [17] and more general perturbations outside spherical symmetry (starting with [18]) have also appeared.
Among the many surprises of the sort mentioned above is the observation of Fermi-Pasta-Ulam (FPU)-like returns of energy configurations in [8]. Conventionally, the Fermi-Pasta-Ulam paradox refers to surprisingly close returns of the energy distribution between linearized modes to its initial configuration observed in the pioneering numerical study [19] of weakly nonlinear oscillator chains (for a review, see [20]). In the context of weakly nonlinear dynamics of AdS, the energy transfer between the linearized modes is accurately captured by the approximate resonant system, also called the Two-Time Framework (TTF) in [8] and the renormalization flow equation or the time-averaged system in [9,10]. One remarkable observation of [8] is that, for initial data with all of the energy in the two lowest modes, this resonant system shows a close return (within a few percent) of the energy distribution to the initial configuration after a few direct and reverse cascades of energy to shorter wavelength modes and back. This is reminiscent of the Fermi-Pasta-Ulam phenomenon (though we emphasize that there are also some significant differences: from the standpoint of the effective resonant system, the dynamics is strongly nonlinear, as the weak nonlinearity parameter has been completely scaled out, while the conventional Fermi-Pasta-Ulam paradox is essentially weakly nonlinear).
A number of developments that have occured since the publication of [8] call us to reexamine the issue of energy returns in global AdS 4 . The complexity of the gravitational dynamics in AdS, even when treated within the resonant approximation, has encouraged examination in [21][22][23][24] of a number of related systems for nonlinear probe fields in AdS, without gravitational interactions. These systems, while possessing much simpler nonlinearities, share many structural features of the gravitational problem, including the patterns in the effective resonant system. It turned out that such probe field resonant systems derived in [21][22][23] possess remarkable analytic features, including exact energy returns for two-mode initial data. We also mention the closely related problems arising from nonlinear Schrödinger equations in isotropic harmonic traps [25][26][27][28][29][30][31][32]. These systems likewise display perfect energy returns, and their close relation to the AdS systems is explained by the nonlinear Schrödinger equation in a harmonic trap emerging as a nonrelativistic limit of nonlinear wave equations in AdS [23].
In view of the exact returns for two-mode initial data which can be analytically demon-strated for a number of systems closely related to weakly nonlinear gravitational dynamics in AdS, a question naturally comes up: could it be that the returns are in fact perfect and the imperfections observed in the simulations of [8] are a numerical artifact? Indeed, the simulations of [8] were performed at a modest numerical precision, which typically involved truncating the system to 30 or 47 lowest modes out of the infinite tower of modes. While a few subsequent studies, including [15] and [37], were performed at higher numerical precision, they did not target the question of energy returns, and made no comments on this issue. We therefore revisit the problem, and examine the returns of two-mode initial data in global AdS 4 in hope of elucidating the question whether such returns are exact or approximate. The course of our investigation involves both numerical and analytic parts. After reviewing in section 2 the basic setup of the resonant approximation to nonlinear gravitational dynamics in global AdS 4 , we dedicate section 3 to a numerical implementation of this resonant dynamical system at a much higher precision than in [8], specifically targeting the question of dynamical returns for two-mode initial data. In this way, we discover strikingly accurate returns that are visually indistinguishable from the perfect returns for related systems observed in [21,23,28,[30][31][32]. Nonetheless, the tiny deviation from perfect returns does not appear to go away completely with increased numerical precision. We highlight the subtleties of numerical simulations of the type of problems we are considering, and are forced to look for analytic clarification of the paradoxical situation we observe. With this goal in mind we turn, in section 4, to regimes in which our two-mode initial data are close to one-mode initial data. In this situation, the problem becomes analytically tractable and one can see that the returns are inexact. This should be contrasted with the scenario of [21,23,28,[30][31][32], where, for related systems, generic two-mode initial data display exact returning behaviors. We thus confirm the FPU-like nature of the energy returns demonstrated by the weakly nonlinear gravitational dynamics in global AdS 4 , which turn out to display a level of return accuracy much more striking than what has been seen in the past literature. In addition, our analytic work explains (at least in the regime of two-mode initial data dominated by one of the modes) the ubiquitous appearance of direct-reverse energy cascade sequences in the weakly nonlinear dynamics of AdS-like systems. Specific integer multiples of these directreverse cascade periods emerge from our considerations, for which the final reverse cascade results in a particularly close return of the energy distribution to the initial configuration. We conclude with a discussion of possible implications in section 5.

Preliminaries
We start with a very brief review of the basic setup of the spherically symmetric AdS-scalar field system and the resonant approximation to its weakly nonlinear dynamics. More details can be found in [9,10,16].
One considers Einstein's gravity with a negative cosmological constant in d spatial dimensions, coupled to a free massless scalar field. The equations of motion are One can consistently truncate to spherically symmetric configurations, corresponding to the metric ansatz where A(x, t), δ(x, t) and φ(x, t) only depend on the time coordinate t and the radial coordinate x, which is defined on the interval [0, π/2). We shall set 8πG = d − 1.
Following [1], we introduce Φ ≡ φ and Π ≡ A −1 e δφ (where dots and primes denote the t-and x-derivatives, respectively), and also the following two predefined functions The equations of motion are then written aṡ Static solutions of these equations are the AdS-Schwarzschild black holes A(x, t) = 1−M ν(x), δ(x, t) = 0 and φ(x, t) = 0, while unperturbed AdS d+1 corresponds to A = 1, δ = φ = 0. Note that because of the spherical symmetry imposed by our ansatz (4), the metric has no propagating degrees of freedom. On each given time slice, the constraint equations (6b) can be integrated to express the metric components in terms of the matter distribution given by φ(x, t) at the same moment of time. Effectively, A and δ can be completely integrated out leaving a nonlinear equation for φ of the form where 2 AdS is the Laplacian in (a non-dynamical) AdS d+1 and S symbolically denotes all the nonlinear terms (which are local in time but nonlocal in space) arising from integrating out A and δ. One can effectively construct S as an expansion in powers of φ, and only the cubic part is significant for our present discussion, since that is what affects the specific weakly nonlinear regime we focus on, see [10]. (One might find it strange that we talk about gravitational dynamics in AdS, while under the assumption of spherical symmetry the metric is non-dynamical and one ends up with a nonlinear wave equation for a scalar field in a fixed AdS background. Nonetheless, there are closely related constructions, which are purely gravitational, without any matter, and utilize a 'squashed' generalization of our ansatz, see [17]. Such extensions effectively result in nonlinear wave equations very similar to the one we have for the scalar field, but now satisfied by the metric components in the absence of any matter.) The problem of weakly nonlinear gravitational dynamics of the AdS-scalar field system under the assumption of spherical symmetry is thus reduced to a complicated nonlinear wave equation in a fixed AdS background. This highlights the relation between gravitational stability of AdS and simpler nonlinear wave equations in AdS that have been considered in the literature as toy models, for example, the λφ 4 theory in nondynamical AdS, see [13,21,23]. Our aim shall be to develop an effective treatment of this dynamics for small fields φ of order ε on long time scales of order 1/ε 2 . This is the regime in which interesting phenomenology, including black hole formation ('turbulent instability'), has been observed in numerical experiments starting with [1].
Before proceeding with weakly nonlinear analysis, one must thoroughly understand the linearized problem, i.e. the equation which under the assumption of spherical symmetry takes the form The eigenvalues for the operatorL are ω 2 n , with and the eigenfunctions are Here, P (a,b) n (x) are Jacobi polynomials of degree n. A remarkable feature, intimately linked to the isometries of AdS d+1 forming the conformal group SO(d, 2), is that all solutions oscillate with integer frequencies ω n . This, in fact, extends outside spherical symmetry, and the most general solution of (8) is time-periodic with period 2π. (The entire tower of eigenmodes of (8) fills an infinite-dimensional representation of SO(d, 2) -and it is also related by a simple transformation to wavefunctions of a superintegrable quantum-mechanical system known as the Higgs oscillator [38,39].) The fully resonant spectrum (10) is highly atypical, and it gives arbitrarily small nonlinearities an opportunity to have large effects, since the impact of resonant interactions between the modes tends to accumulate over time. In naive perturbation theory in terms of power series in ε, this feature is reflected in a breakdown of perturbative expansions on time scales of order 1/ε 2 , as noted already in [1]. More specifically, if one tries to develop solutions of (7) as solutions of (8) of order ε plus corrections of higher orders in ε, already at order ε 3 one ends up with 'secular' terms growing as ε 3 t, which overpower the leading (linear) term on time scales greater than 1/ε 2 . The origin of such terms can be directly traced back to the presence of resonances in the linearized spectrum.
The discussion above shows that naive perturbative expansions are not a valid way to approximate the weakly nonlinear dynamics of (7) on the relevant timescales. A particularly viable alternative approach is time-averaging, introduced for studies of AdS systems in [8] and developed analytically in [9,10] (various names were used for it). One starts by expanding solutions of (7) in terms of linearized eigenmodes (11) as φ(x, t) = −iε ∞ n=0 ω n α n (t)e −iωnt −ᾱ n (t)e iωnt e n (x), and equivalently rewriting (7) as a system of equations for α n . This system of equations for α n is of course completely equivalent to (7), but it matches what is known as the 'periodic standard form' in mathematical literature, and the time-averaging procedure may be applied to simplify this equation in a way that does not affect its accuracy on timescales of order 1/ε 2 . More specifically, substituting (12) in (7) and projecting on e n , one expresses iα n as a cubic combination of α k andᾱ k . Most of the terms in this cubic combination come with oscillatory factors inherited from the explicit time dependences in (12). Time-averaging (backed by precise mathematical theorems) discards all oscillatory terms, retaining only the resonant terms in which the oscillatory factors cancel. This is guaranteed not to affect accuracy on time scales of order 1/ε 2 and results in a resonant system of the form where, from now on, dot will stand for a derivative with respect to the slow time ε 2 t, which we shall call simply t from now on. A few remarks are in order: • The dependence on the small nonlinearity parameter has been completely scaled out of (13) with the introduction of the slow time, and there are no small parameters left. In fact, (13) enjoys a scaling symmetry: if α n (t) is a solution, so is λα n (λ 2 t), for any λ.
• The nonlinear physics of the problem is completely encoded in the interaction coefficients C nmkl which are a set of numbers expressed through the mode functions (11) and their spatial derivatives, and depending on the specific form of nonlinearity in (7). For the gravitational system we are considering here, the explicit expressions for C are very complicated and can be found in [10] (where these coefficients are denoted as T , R and S, depending on the number of coincident index values). In the Appendix, we extend them to the case of a massive scalar field, which will be discussed in section 5.
• The frequencies ω n can be absorbed in C by redefining α n , but it is often convenient to keep them explicit. Note that ω n + ω m = ω k + ω l is equivalent to simply n + m = k + l.
• A number of simpler related problems (nonlinear wave equations in non-dynamical AdS, variations of the Gross-Pitaevskii equation for Bose-Einstein condensates in harmonic traps) lead to resonant systems that only differ from (13) by specific values of the numerical coefficients. Some of these systems demonstrate perfect energy returns in the evolution of two-mode initial data [21,23,28,30].
• Purely from resonance analysis, (13) could have extra terms (for example, with three α's and noᾱ's), but such terms can be shown to vanish specifically in AdS due to special selection rules [9][10][11][12]. Correspondingly, there is an emergent U (1) symmetry rotating the phases of α n by a common shift, which is manifest in (13) but absent in (7).
• In immediate relation to the extra U (1) symmetry resulting from selection rules, there is an extra conservation law [10,13,14] for the associated 'particle number' in addition to the total 'linearized energy' generically conserved by resonant systems. They are associated with the symmetry transformations α n → e iϕ α n , α n → e iωnθ α n , respectively.
Our main objective is to examine solutions to (13) starting with two-mode initial data In particular, the energy distribution between the modes may be quantified by and we may scan the evolution of two-mode initial data for return moments when the deviation from the initial energy distribution is small. Due to the scaling symmetry of (13) we can transform any initial data (17) to satisfy Thereafter, it is sufficient to study the quantity If ∆(t) is 0, joint conservation of (14) and (15) guarantees that E 0 (t) = E 0 (0) and E 1 (t) = E 1 (0), i.e., we have found a perfect return, while 0 < ∆(t) 1 signifies an accurate but imperfect return. As we have already remarked, for a number of related resonant systems solved in the literature, ∆(t) periodically drops to 0 for any two-mode initial data of the form (17), signalling exact returns [21,23,28,30]. In our present work, we shall look closely into similar return phenomena for the gravitational resonant system (13), which is only possible numerically.

Numerics
Before describing the results of our own numerical experiments, we briefly summarize what has been done in this area before. Already in [8], one of the themes was the observation of accurate but imperfect returns of energy to the initial configuration for two-mode initial data in the resonant system for gravitational AdS 4 perturbations (there referred to as the Two-Time-Framework approximation). These numerical studies were performed at the resolution available at that time, which amounted to truncating the resonant system to the lowest 30 or 47 modes, and returns with precision of a few percent were observed. The evolution was seen to proceed in a sequence of direct and reverse energy cascades, and the third reverse cascade led to a more accurate return to the initial energy distribution than the first two. Subsequently, a few developments occured, including the derivation of analytic expressions for the interaction coefficients of the AdS resonant systems in terms of the AdS mode functions in [9,10] and an algorithm to convert these expressions into explicit (but very complicated) functions of the mode numbers in [37]. While this has allowed simulations of the resonant dynamics with much higher precision, the question of energy returns has not been properly readdressed in the literature following these developments.
We reiterate our reasons to revisit the question of energy returns in AdS 4 . Our current perspective is quite different from the predominant views at the time [8] was written. Indeed, in the years that have passed, a number of examples have emerged, where resonant systems closely related to the one studied in [8] display exact, rather than approximate, energy returns. These include nonlinear probe fields in AdS [21,23] as well as nonlinear Schrödinger equations in harmonic potentials that arise from AdS systems in a nonrelativistic limit [28,30,32]. This puts the results of [8] in a different light and calls for their re-examination. Note that the limited precision of [8] (mode number cut-off at a few dozens modes) makes it impossible to distinguish exact and approximate returns. In particular, in application to the systems of [21,23,28,30,32], where the returns are analytically known to be exact, such simulations would have indicated approximate returns (the imperfection in this case being a pure artifact of truncating the infinite-dimensional system).
To elucidate the issue of energy returns in global AdS 4 , we have performed simulations with up to 500 modes studying the dependence of the return accuracy on the mode number cut-off as well as on the arithmetic precision. For some initial energy distributions between the two lowest modes, the returns we observe are strikingly accurate. Fig. 1a provides an illustration in this regard: we study the evolution of two-mode initial data which are fairly generic (not particularly close to one-mode initial data) over a sequence of nine directreverse cascades. After each three cascades, we have a return to the initial configuration that is visually indistinguishable from exact, and furthermore the whole pattern periodically repeats after each three oscillations in a way that is visually indistinguishable from exact periodicity. All of this happens in (our high-precision numerical truncation of) the infinitedimensional nonlinear resonant system (13) that has no small parameters! Note that the plot given in fig. 1a is visually extremely similar to an analogous plot (see fig. 1b) in a closely related resonant system called the LLL equation, where the energy returns have been analytically proven to be exact [28]. Fig. 2 provides four illustrations where we can observe the behavior of the energy flow for different initial energy distributions between the first two modes: fig. 2a is particularly close to the second mode and the evolution appears periodic  after each direct-reverse cascade, while in fig. 2b and fig. 2c the apparent periodicity only emerges after three direct-reverse cascades. Finally fig. 2d describes initial data for which the evolution visibly deviates from perfect periodicity, although it is remarkable that the deviations are still pretty small.
And yet, spectacular as some of these returns are, they are not exact. This by itself is not conclusive, since returns in numerical simulations of truncated systems would not have been exact even if returns in the infinite-dimensional dynamical system were exact. We thus have to quantify the dependence of the accuracy of returns on the initial data, as well as on the mode number cut-off, arithmetic precision and other numerical imperfections. To set the stage, we present in fig. 3 plots showing the dependence of two essential quantities on the initial energy of mode 0 for two-mode initial data. The first quantity (fig 3a) is the minimal value of ∆ defined by (20) over the first three direct-reverse cascades, which quantifies the energy return precision (evidently, in computing the minimum we exclude the initial part of the first direct cascade where ∆ is small simply because the system has not had enough time yet to deviate from the initial conditions). We see that the worst returns are around E 0 ≈ 0.6 at a level slightly better than 10%, while for most other initial data away from that value the returns have precision better than 1%. The second quantity we plot (in fig 3b) measures the strength of the turbulent cascade starting with the given twomode initial data. The energy spectrum of configurations that undergo regular evolution is suppressed exponentially for large mode numbers, E n (t) ∼ n γ(t) exp(−ρ(t)n) with ρ(t) > 0, and turbulent singularity formation, in particular, would correspond to ρ(t) hitting zero at a certain time [15,33,34]. We quantify the strength of a turbulent cascade by determining how small ρ(t) becomes during the time interval of interest. Therefore, we first fit the logarithm of the spectrum log E n (t) to −ρ(t)n+γ(t) log n and then plot the minimal value of ρ attained over the first three oscillations. We see that the energy return imperfections are especially strong in the region where the turbulent cascade is strong as well. A strong turbulent cascade means, in particular, a stronger sensitivity to the mode number cut-off (as one would need Figure 2: This sequence of time evolutions of four two-mode initial data (17,19) shows that (a) the energy transfer between modes appears to evolve in a periodic way when the initial condition is close to mode 1 (a), that for larger E 0 the evolution is apparently periodic only after three oscillations (b,c), and that for equal-energy initial conditions, the returns are visibly not exact, but still more accurate after three oscillations (d).
to appreciably excite modes above the cut-off to keep the evolution exact, but those modes are excluded from simulations). This makes the problem of interpreting the deviations from exact returns rather subtle. Since we are entering the realm of precision questions, we must categorize and quantify the uncertainties incurred by our numerics. There are of course generic small errors arising from numerical integration of the equations of motion, but we feel that the following three aspects are particularly important, since they are specific to the sequences of direct and reverse cascades of energy we study: 1) Mode number truncation: Restricting to a finite number of modes is expected to have little effect if the energy mostly remains locked within a set of low-lying modes. In our case, an energy cascade develops, transferring the energy to higher modes. While the cascade is self-limiting (unlike in the fully turbulent case of [15], where ρ(t) develops a zero), it can be quite strong, as evident from fig. 3b. What we observe in numerical experiments is that, as the cascade hits the mode number cut-off and then starts receding, a comb-like pattern of spikes in the spectrum forms. This pattern undergoes an evolution of its own and remains visible in future direct-reverse cascade oscillations. As we can observe in fig. 4a, these comb patterns are pure numerical artifacts, and do not reflect the true dynamics plotted against the initial energy of mode 0 for two-mode initial data (17,19). Small values of the curve in (b) correspond to strong turbulent cascades of energy.
of the infinite-dimensional model, since the parts of the spectrum overrun by the comb artifacts disagree between simulations with different mode number cut-offs. Such artifacts are generically visible for numerical simulation of resonant systems, including those for which exact analytic solutions starting with two-mode initial data are known (and do not display such comb-like features). While the total energy in the comb-like artifacts is very small at the moment of their formation, it is difficult to predict how they affect future directreverse cascade sequences, since the system is nonlinear and prone to chaotic divergence of trajectories.
2) Arithmetic precision: There is a very peculiar way in which direct-reverse cascade sequences amplify numerical errors, including the most basic rounding errors caused by the finiteness of arithmetic precision. Indeed, assume that in the exact solution a particular mode number k oscillates between energies E max k and E min k . In our scenario of near-perfect returns E min k E max k , while for perfect returns E min k = 0. Since we are integrating the equations numerically and with finite precision, E max k contains a relative error given by at least the arithmetic precision, for example 10 −15 . As the cascade recedes and the energy drops, the absolute error in E k cannot decrease. Therefore, when one arrives at the bottom of the reverse cascade, the energy is E min k , but the absolute error is still 10 −15 E max k (in particular, values of E min k below 10 −15 E max k can never be reached, even if the exact solution corresponds to E min k = 0). The relative error is now 10 −15 E max k /E min k , much greater than the rounding errors themselves. This relative error will be transported upstream in the next direct cascade and will result in a large absolute error 10 −15 (E max k ) 2 /E min k at the peak of the cascade. Thus we can see that repeated transport of absolute errors downstream and relative errors upstream in direct-reverse cascade sequences leads to strong amplification of numerical imprecision (and would of course compromise exact returns even if they were present in the underlying dynamical system). In practice, we can see the formation of flat 'shelf' artifacts in the high mode number part of our spectrum as the first cascade recedes, as in fig. 4b (these artifacts precisely reflect the inability of the spectum to go below 10 −15 E max k in numerical simulations). Such artifacts continue to evolve in future spectrum  : Energy spectrum of the two mode initial data (17,19) at the bottom of the first reverse cascade for different mode number cut-off. In (a) the direct cascade is strong, in which case the most relevant effects are produced by the cut-off. In (b) the direct cascade is weak, and 'shelf' artifacts (discussed in the main text and originating from the finite arithmetic precision) dominate the imperfections of the spectrum.
oscillations. As the near-perfect returns we observe occur after three direct-reverse cascade cycles, assessing the ultimate impact of these imperfections on the return precision is subtle. 3) Interaction coefficients: Small errors in the interaction coefficients C can produce important errors in the evolution governed by (13). This is particularly important when high modes are involved, as computing the interaction coefficients typically requires numerical integrations of oscillatory functions with frequencies that grow with the mode numbers. These integrals have to be evaluated for each resonant quartet of modes satisfying n+m = j+ k, and the number of such quartets grows like O(N 3 max ), where N max is the mode number cutoff. In addition, in order to obtain accurate values of the coefficients C when N max is large, one needs to evaluate their expressions on a very dense grid and with increased arithmetic precision, which further increases the burden. To avoid these problems we computed the fully analytic expression for C in AdS 4 described in [37], which enables us to evaluate C safely for high modes. The remaining burden is then the computational cost of the simulations, which makes it hard to go beyond N max = 500.
We remark that improving on either 1) or 2) beyond the level of our current simulations would be very demanding in terms of computational costs. While a break-through in our numerical precision does not appear viable at this moment, we have performed some basic comparison of the return precision computed at different values of numerical approximation parameters (as well as a study of artifacts that we have briefly summarized in the passages above). In fig. 5, we show the dependence of return precision on the number of modes included in our simulations. As one can see, the return precision substantially increases as the cut-off is raised from values around 50 (typical of the simulations of [8]), but after 200 stabilizes at a small nonzero value. At a naive level, this can seen as an indication that the returns are not exact in the underlying dynamical system (13). One has to keep in mind that other factors hindering perfect returns (finite arithmetic precision, possible chaotic enhancement of artifacts over the course of three direct-reverse cascade oscillations) have not been taken (and are very difficult to take) into account. We conclude that our numerical analysis indicates strikingly accurate, but likely imperfect energy returns for two-mode initial data, though it is not possible as of now to control all possible sources of imperfection. There is, however, one particular aspect that can be further elucidated. Our motivation to take a closer look at energy returns in global AdS 4 has largely come from the similarity of the observed dynamics to the known related analytically tractable cases [21,23,28,30,31] where the returns are exact. This similarity extends in a few other ways, in particular the ultraviolet parts of the spectra remain approximately exponential for all times (only the slope of their logarithmic plot changes). However, for the systems of [21,23,28,30,31], energy returns for any two-mode initial data (17) are exact. This type of dynamics can be excluded for the resonant system (13) by inspecting solutions close to mode 0 or mode 1, which are analytically tractable. This will show that the resonant dynamics of AdS 4 cannot be strictly in the same class as those of the similar resonant systems of [21,23,28,30,31], which do exhibit perfect energy returns.

Analytics
While direct analytic investigations of the resonant system (13) are beyond practically imaginable limits, there are special regimes in which this system is analytically tractable. First of all, there are exact single-mode solutions, in which all amplitudes are zero except for one chosen mode [1]. The vicinities of such single-mode solutions form stability islands [1,3,4,8,37]. While it is common to linearize in the vicinity of such single-mode solutions and their generalizations [22,37], which results in linear systems that are potentially tractable, our approach here will also gain leverage by relying on single-mode solutions, but in a way different from linearization, and much more useful for our purposes.
Instead of linearizing around single-mode solutions, we shall assume that all other modes are exponentially suppressed in proportion to their distance from the dominant mode. We shall see that, in the limit when this exponential suppression becomes strong, the equations dramatically simplify. Such techniques have been applied in [36] to the analysis of stationary solutions of (13) and related resonant systems. While the equations are still nonlinear (unlike in the linearization approach mentioned above), they can be solved iteratively starting with the dominant mode. This structure is very convenient for our purposes, since we are trying to understand whether the energy returns may be exact. If the returns are exact, they are exact for all modes, while a failure of exact returns shall be seen in the approach we are adopting in a finite number of steps, since the discovery of any given non-returning mode guarantees that the returns in the whole system cannot be exact. This procedure has to be repeated two times, for solutions dominated by mode 0 and mode 1, respectively.

Solutions dominated by mode 0
We shall assume that the spectrum is exponentially suppressed as one moves away from the dominant mode 0: with δ 1. Then, after the redefinition C new At leading order in δ, only the terms with m = 0 survive: While this equation is still nonlinear, it offers a tremendous simplification over (13), since it can be solved recursively: once solutions for the modes up to q n have been constructed, finding q n+1 amounts to solving a single linear ODE. We display this structure by writing out the equations for the first few modes: By applying the scaling symmetry of (13), one can always set |q 0 | 2 to 1. In order to elucidate whether energy returns of the two-mode initial data of the form (17) are exact, we solve (23) with the initial conditions q 0 (0) = 1, q 1 (0) = 1, q n≥2 = 0. (The phases of q 0 (0) and q 1 (0) can be adjusted as necessary by the symmetry transformations (16), while the magnitude of q 1 (0) can be set to 1 by defining δ = |α 1 (0)|.) The solutions for q 0 and q 1 are Note that the energies in these first two leading modes are time-independent, which implies that close to mode 0 the cascade to higher energies is necessarily very weak (the situation will be more interesting for solutions dominated by mode 1). For higher q n , one then obtains The structure of solutions is easily understood recursively: the right-hand side consists of terms oscillating with frequencies given by linear combinations (with integer coefficients) of C 0000 , C 1010 , ..., C n−1,0,n−1,0 . Hence, q n will consist of terms that oscillate 1 with frequencies given by linear combinations (with integer coefficients) of C 0000 , C 1010 , ..., C n0n0 . It is only if all q n constructed in this manner (for two-mode initial conditions) oscillate with a common period that one gets exact energy returns. For the resonant system (13), one can straightforwardly proceed with this algorithm, after having computed a few low-lying interaction coefficients from the complicated formulas given in [10]. The result for q 2 (t), which is shows that |q 2 (t)| 2 is periodic with period T 2 = 10 123 π 2 . Thus, energy returns are perfect at this order, but the period T 2 is broken at the next order given by where In our particular case their specific values are: λ 1 = 3057 70π , λ 2 = 267 14π meaning that at |q 3 (t)| 2 the energy flow has a periodicity of T 3 = 140π 2 3 , so that T 3 is exactly 574 times greater than T 2 . Inspecting the next step of our iterative solution, we conclude that |q 4 | 2 is also periodic with T 4 = 420π 2 .
The pattern of common multiple periods will generalize to the higher modes q n because of the rational relations between the interaction coefficients C nmkl in global AdS 4 . Namely, one can apply the strategy of [37] to derive explicit, complicated formulas for C nmkl as functions of n, m, k, l. This leads, in particular, to the conclusion that 2πC 0n0n = 108 + 90n − 33n 2 − 58n 3 − 21n 4 − 2n 5 n(n + 1)(n + 2)(n + 3) Hence, all C 0n0n are rationals divided by π, and therefore q n , which oscillate with frequencies given by rational combinations of C 0n0n , have common periods of the form π 2 times a rational number. (As a direct consequence, |q n | 2 , which are what is important for our topic of energy returns, also have common periods of the same form, typically somewhat shorter than the common periods of q n .) The common period of q 0 , q 1 , . . . q n grows rather rapidly with n (for instance, we get estimates of order 10 22 for n = 25).
Our investigations indicate that common periods survive even if subleading corrections in δ are included (in the above analysis, we have only considered the leading contributions in δ for each mode), though we shall not explore this question systematically. The rational relations between the interaction coefficients are such that for any set of modes q 0 , q 1 , . . . q n and any given order of the small δ expansion, a common period of q 0 , q 1 , . . . q n exists. This suggests that, in the full resonant system (13), energy returns of arbitrary precision for initial data dominated by mode 0 will be seen, provided that one waits for an appropriate large number of direct-reverse cascade oscillations.
To summarize, the expansion of equation (13) in powers of δ around mode 0 enabled us to verify that the energy transfer between the modes is not exactly periodic for initial conditions close to mode 0. Nevertheless, there are rational relations between the AdS 4 interaction coefficients such that arbitrarily precise returns to the initial energy distribution starting from two mode initial data dominated by mode 0 occur if one waits long enough. This is very much in the spirit of the original FPU paradox.

Solutions dominated by mode 1
We shall now assume a hierarchically organized spectrum dominated by mode 1: with δ 1. We then get from (13), after redefining C new nmkl = C old nmkl / √ ω n ω m ω k ω l , the following equation for mode 0: and for all higher modes, (38) Retaining only the leading terms, we obtain and iq n =q 1 n k=1 C n1k,n+1−k q k q n+1−k +q 0 To appreciate the structure, we write out the first few equations explicitly: One first solves (42), where one can set |q 1 | 2 = 1 as before using the scaling symmetry. After that, (41) and (43) form a system of two coupled linear equations for q 0 and q 2 . Once q 0 , q 1 and q 2 have been thus obtained, all the higher equations are solved recursively one-byone, in a manner completely analogous to what we have previously described for solutions dominated by mode 0. Implementing this solution in practice for q 0 (0) = 1, q 1 (0) = 1 and q n>1 (0) = 0 (as before, the complex phases can be eliminated using the symmetries (16), while the magnitude of q 0 (0) can be fixed by redefining δ), one gets In what follows, we assume λ ∈ R, which is the case for the system we are dealing with, as well as for all the systems studied in [21,23,28,30,31]. Note that, here, one gets nontrivial flows of energy among the first subleading modes, unlike the case dominated by mode 0. These flows of energy are always periodic with period π/λ, as far as the three lowest modes are concerned. The first order at which violations of periodicity may enter is in the mode q 3 , which is described by the following solution: with For the resonant system (13), direct computation yields λ/γ = 7 √ 3665/149 ≈ 2.84. (Note that irrational numbers emerge here already in the periods of low-lying modes, unlike the case of initial data dominated by mode 0 we have considered previously.) Suppose for a moment that one had λ/γ = 3. Then E 3 ∼ |q 3 | 2 computed from (49) would have been proportional to a + be iλt + ce 2iλt + de 4iλt/3 2 , which oscillates with period of 3π/λ, which is thrice the period of |q 0 | 2 , |q 1 | 2 , |q 2 | 2 . One would thus have found exact returns after three direct-reverse cascades as far as the first four modes are concerned. Now, for the actual resonant system we study, λ/γ is of course not exactly 3. This shows that the third energy return cannot be exact, at least for initial data close to mode 1, no matter how close we get to exact returns in our numerical simulations. At the same time, the fact that 2.84... is close to 3 explains why we are seeing very accurate returns after three direct-reverse cascades (and also why the returns after the first two direct-reverse cascades are less exact), as in figs. 1a and 2d.
Note that our analysis gives a nice perspective on why direct-reverse cascade oscillations are ubiquitously seen in numerical simulations of various resonant systems. For solutions dominated by mode 1, for example, the first three modes perform an infinite sequence of direct-reverse cascades, while the higher modes hold only a small amount of energy and provide cosmetic modifications to the cascades. (Furthermore, even these higher modes oscillate with frequencies comparable to the lowest modes.)

Discussion
We have revisited the issue of energy returns to the initial configuration for two-mode initial data in the resonant approximation to weakly nonlinear gravitational dynamics of the AdS 4scalar-field system. Having performed numerical simulations with much higher precision than what has been previously seen in the literature, we have observed returns of striking accuracies, exemplified by fig. 1. The numerics also provided indications, however, that the small imperfections we see cannot be purely due to numerical artifacts. To elucidate the situation, we have performed an analytic study of solutions dominated by one of the two lowest modes and proved that the accurate returns we have observed numerically are inexact in this limit dominated by one of the two modes. This should be contrasted with a scenario observed in the recent literature on related systems [21,23,28,30,31], where perfect returns occur for all two-mode initial data of the form (17).
As usual in FPU-like situations, it is natural to expect that the near-perfect returns arise due to proximity to another dynamical system for which the returns are exact. One could first try to improve the quality of returns by variation of parameters of the model, such as the dimension of AdS or the mass of the scalar field. We have explored this scenario in the context of our analytic treatment of solutions dominated by mode 1 in section 4.2. We have seen that by adjusting the mass of the scalar field, it is possible to make the return of mode 3 exact ( fig. 6). This evidently qualitatively improves the precision of returns for initial data near mode 1, as the discrepancy is now in the strongly suppressed modes starting from mode 4, but the scenario still falls short of providing exact returns. We also note that a nonrelativistic version of the AdS dynamics (which technically corresponds to the limit of infinite scalar field mass) can be analytically proved to display perfect returns [32] in AdS 5 , rather than in AdS 4 (returns in AdS 5 at finite masses, on the other hand, are not close to being perfect). More broadly, one could look for arbitrary small modifications of the interaction coefficients C in (13), irrespectively of whether they originate from standard physically motivated PDEs, and perfect returns in any such system would be sufficient to clarify the origin of near-perfect returns in the (physically motivated) resonant system (13). In [31], a very large class of resonant systems displaying perfect returns for two-mode initial data has been constructed. This class is characterized, in particular, by the following identities satisfied by the interaction coefficients: n+m k=0 k 2 f k f n+m−k f n f m C nmk,n+m−k = c 2 (n 2 + m 2 ) + c 1 nm + c 0 (n + m), where f n can be either 1/ √ n! or G(G + 1) · · · (G + n − 1)/n!, and G, c 0 , c 1 , c 2 are arbitrary constants. Our numerical comparisons indicate that the interaction coefficients of the resonant system (13) do not appear close to satisfying such relations. One could look for modifications of (52) approximately satisfied by the AdS 4 interaction coefficients. Our preliminary study shows that some polynomial-type summation identities are satisfied with a good precision, but we feel that it is premature to judge what exact dynamical implications follow from such approximate identities.
The AdS/CFT paradigm provides an intriguing potential link between our results and dynamics of conformal field theories (CFTs) on spatial spheres. In particular, one would expect that very close returns to the initial state must occur at low energies in CFTs that can be accurately approximated in an appropriate 'holographic' limit by a single scalar field in the AdS bulk coupled to gravity. It would be very interesting to delineate this class of systems more precisely, and build explicit connections to the sort of dynamics we have described in this article. There is recent literature [40,41] providing detailed links between evolution of CFT states and weakly nonlinear gravitational dynamics in the bulk, which is likely to be useful in this regard.
In analogy with the explanation for the massless scalar field given in section 2, the current model has the resonant condition ω n + ω m = ω k + ω l , which through (65) is equivalent to n + m = k + l. However there could be two more resonant channels, ω n = ω m + ω k + ω l and ω n + ω m + ω k = ω l , equivalent to n = m + k + l + ∆ and n = l − m − k − ∆, respectively. If ∆ is integer these last two conditions can be satisfied and a priori two new terms must be included in the system of equations (13): In our situation we did not perform an analytic study of Q ijkl and U ijkl as in [10], where it was proven that for a massless scalar field (∆ = d) these coefficients vanish, but numerical calculations suggest that they also vanish for nonzero masses. On the other hand, when ∆ is not integer, the conditions n = m + k + l + ∆ and n = l − m − k − ∆ are not satisfied for any combination of the indices. Therefore, these interaction channels disappear upon time-averaging. We thus see that, for any ∆, the relevant dynamics is governed by equation (13) through T l , R il and S ijkl . We note that the given expressions for T l , R il and S ijkl are exactly the same as in [10], where all derivations are specialized to m 2 = 0. However, if we were to transform these integrals using integration by parts to the form of [9], the resulting expressions would have differed from those of [9] by additional terms with explicit dependence on m 2 .