Ultraviolet asymptotics for quasiperiodic AdS_4 perturbations

Spherically symmetric perturbations in AdS-scalar field systems of small amplitude epsilon approximately periodic on time scales of order 1/epsilon^2 (in the sense that no significant transfer of energy between the AdS normal modes occurs) have played an important role in considerations of AdS stability. They are seen as anchors of stability islands where collapse of small perturbations to black holes does not occur. (This collapse, if it happens, typically develops on time scales of the order 1/epsilon^2.) We construct an analytic treatment of the frequency spectra of such quasiperiodic perturbations, paying special attention to the large frequency asymptotics. For the case of a self-interacting phi^4 scalar field in a non-dynamical AdS background, we arrive at a fairly complete analytic picture involving quasiperiodic spectra with an exponential suppression modulated by a power law at large mode numbers. For the case of dynamical gravity, the structure of the large frequency asymptotics is more complicated. We give analytic explanations for the general qualitative features of quasiperiodic solutions localized around a single mode, in close parallel to our discussion of the probe scalar field, and find numerical evidence for logarithmic modulations in the gravitational quasiperiodic spectra existing on top of the formulas previously reported in the literature.

Ultraviolet dynamics is of crucial importance to the problem of evolution of AdS perturbations. Instability and black hole formation occur through non-linear energy transfer to very short wave-length modes, thus the question of ultraviolet asymptotics in weakly non-linear perturbation theory comes to the fore; see, for instance, the recent paper [10].
An important aspect of the non-linear perturbative AdS dynamics is the existence of periodic and approximately periodic solutions. For example, an algorithm to construct exactly periodic solutions order-by-order in the field amplitude has been presented in [4] (and extended to massive fields in [15]). More broadly, one can consider quasiperiodic solutions, in which the AdS normal modes oscillate with constant amplitudes at leading order in the field amplitude on time scales inversely proportional to the square of the field amplitude [6]. Since these time scales are precisely the ones on which collapse occurs in numerical experiments, as well as the shortest time scales on which substantial energy drift may occur in the weakly non-linear regime in AdS settings [6,7], the existence of such solutions presents an obstruction to efficient collapse. Islands of stability in the phase space of AdS perturbations are thought to be anchored by such quasiperiodic solutions [9,16].
At leading order, effectively resummed perturbation theory [6][7][8] takes the form of equations describing slow drifts (on time scales inversely proportional to the square of the field amplitude) of the amplitudes and phases of the AdS normal modes. This leading order approximation can be derived by many techniques (multi-scale or renormalization group analysis, averaging) and goes by a few different names in the literature ('two-time framework,' 'renormalization flow equations,' 'time-averaging,' 'resonant approximation'). Quasiperiodic perturbations appear in this context [6] as solutions for which the normal mode amplitudes remain constant in time, despite the presence of weak non-linear interactions. This condition can be reduced to an infinite system of cubic algebraic equations, solutions of which were analyzed numerically in [6], revealing exponential damping for high frequency modes. It was pointed out in [7] that the abundance of quasiperiodic solutions is intimately related to selection rules in perturbation theory, which, in turn, give rise to conserved quantities [8]. The solutions were explicitly parametrized by conserved quantities in [9] and numerically shown to be stable in [16].
Here, we commence our studies of the ultraviolet behavior of the AdS perturbation theory by analyzing the ultraviolet asymptotics of the spectra of such quasiperiodic solutions. We do it both for self-interacting probe scalar field dynamics [14] and for gravitational interactions [6]. For the self-interacting probe scalar field, we develop a fairly complete analytical treatment that recovers the essential features of the asymptotic spectra in a way that agrees with direct numerical evaluations. Specifically, one gets amplitude spectra approaching an exponential decay multiplied by the mode number for high-frequency modes. For the gravitational case, both the numerical behavior and the structure of the equations are more involved. We display the asymptotic properties of the coefficients in the equation determining the quasiperiodic spectrum, and explain the general qualitative features of the solutions, such as the existence of two spectral branches previously pointed out numerically [6]. We are unable, however, to draw clear analytic conclusions on the spectral asymptotics from these equations. At the same time, our numerical study of the gravitational quasiperiodic spectra gives evidence for a hierarchy of slow logarithmic modulations superposed on the simple fitting formulas mentioned in [6].
We present our results below by first reviewing (in section 2) the improved perturbation theory in AdS at first non-trivial order, leading to an accurate description of the dynamics on time scales inversely proportional to the square of the field amplitude, as well as the construction of quasiperiodic solutions in this context. We then describe our analysis of the quasiperiodic spectra asymptotics, first for a self-interacting scalar field (section 3), and then for the gravitational case (section 4). Other aspects of ultraviolet asymptotics of AdS dynamics are discussed in our parallel work [17], where, in particular, we analyze the ultraviolet behavior of the interaction coefficients of the time-averaged equations in arbitrary dimension.

Time-averaged effective theory and quasiperiodic solutions
We shall briefly review the construction of effective flow equations describing the slow energy transfer between the AdS scalar field modes, starting with the simpler set-up of a self-interacting scalar field [13,14] in AdS d+1 . A considerably more systematic discussion of these matters (phrased in a context including full gravitational non-linearities) can be found in [7,8].
We start with the global AdS d+1 metric written as where the radial x coordinate runs between 0 and π/2 and dΩ 2 d−1 is a (d − 1)-sphere line element. The equations of motion for a spherically symmetric probe scalar field φ(x, t) with a φ 4 /4! non-linearity take the form Expanding this in normal modes φ(x, t) = n c n (t)e n (x) yields where P n are Jacobi polynomials, ω n = d + 2n is the peculiar integer (fully resonant) AdS d+1 frequency spectrum, and are the 'interaction coefficients' that will play a key role in our discussions below.
One can now zoom in on the small amplitude regime, assuming amplitudes of order , and switch to the 'interaction picture,' introducing slowly varying complex amplitudes α k (which would be exactly constant if the non-linearities were neglected): c n = ε α n e iωnt +ᾱ n e −iωnt ,ċ n = iεω n α n e iωnt −ᾱ n e −iωnt , and hence (7) By the standard lore of averaging [8,13], when ε is small, only non-oscillating terms on the right-hand side can contribute significantly. This implies that, out of the sum over jkl, only terms satisfying the resonance condition ω n = ±ω j ± ω k ± ω l survive (with all the three plus-minus signs a priori independent). Furthermore, the C coefficients turn out to vanish (see footnote 3 of [8] for a simple proof) for all such resonance conditions other than ω n + ω j = ω k + ω l , i.e., n + j = k + l (and similar relations obtained by interchanging j with k or l). One thus arrives at the following effective time-averaged theory accurately approximating the original equations on time scales of order 1/ε 2 : where 3/3! = 1/2 accounts for two other contributions interchanging j with one of the two other summation indices in (7). One may look, following [6,14], for solutions of the time-averaged equation preserving the normal mode amplitudes |α k |, i.e., solutions in which there is no energy flow between the modes, despite the presence of non-linearities. 1 A large class of such solutions is given by imposing α n = a n e −iε 2 βnt , where β n = β 0 + n(β 1 − β 0 ), and a n , β 0 , β 1 are time-independent. This ansatz converts the differential equations (8) to a system of algebraic equations for a n , β 0 , β 1 : 2ω n a n (β 0 + n(β 1 − β 0 )) = 1 2 jk C njk,n+j−k a j a k a n+j−k .
It is our main objective in this paper to examine what we can say analytically about solutions to the quasiperiodicity conditions of this sort for the time-averaged system. For the self-interacting scalar field, we shall present a rather complete picture of the ultraviolet properties of solutions with spectra localized around a single mode. We shall also present some considerations for the case with full gravitational interactions and a free massless scalar field [1]. The ideology of our treatment is very similar though both technical aspects and phenomenology are more complicated and the amount of progress we can make is considerably more modest. In this case, one considers a spherically symmetric scalar field coupled to a dynamical spherically symmetric metric of the form with A and δ being functions of x and t. Due to spherical symmetry the metric is determined on each constant time slice by the constraint equations in terms of the matter distribution (i.e., the scalar field profile). One hence essentially obtains an effective theory for φ, albeit with non-linearity much more complicated than in (2). Improved perturbative expansions can be constructed for this system [6][7][8], though doing this analytically [7,8] requires considerable effort. See, in particular, [8] for a formulation in the language of timeaveraging. Once the effective averaged description of the slow energy transfer has been constructed, one can look for quasiperiodic solutions along the lines we described above for the probe scalar field.
3 Non-linear probe scalar field in AdS

Quasiperiodic solutions dominated by a single mode
The conditions (9) defining the spectra of quasiperiodic solutions form an infinite system of cubic algebraic equations, clearly beyond direct reach of any analytic methods.
A particular approach to constructing solutions of (9) has been proposed (in the gravitational context) in [6]. Namely one looks for solutions with one mode being much larger than any other mode. (There are obvious exact solutions with only one mode excited, and single-mode-dominated solutions can be seen as perturbations of such obvious single-mode solutions.) In [6], such solutions were constructed numerically. Here, we shall focus on solutions that are strongly localized in mode number, for which, as it turns out, a fairly complete analytic treatment can be given.
Consider first for simplicity a solution dominated by mode 0, namely a 0 = 1 (there is a scaling symmetry in the equation that always allows to set the amplitude of any given mode to 1), a 1 = δ 1. One can then see that there is a consistent ansatz a n = q n δ n , which results in the leading powers on the right-hand-side and left-hand-side matching in all of the equations (9). The procedure of constructing the solutions then goes as follows. In each equation (9) one only keeps the leading powers of δ, matching its coefficient on the two sides of the equation. For equation number 0, this defines β 0 , for equation number 1, this defines β 1 , each subsequent equation determines q 2 , q 3 , etc. Importantly, restricting to the leading power of δ (strong localization limit) in each equation truncates the infinite sums in (9) to a finite number of terms. Specifically, (One of the summation indices in (9) had to be set to 0 to ensure a leading, rather than subleading, power of δ.) Equation (11) presents a tremendous simplification over (9) and can be analyzed quite thoroughly. We shall do that after examining the large index properties of the C coefficients. For solutions dominated by mode number J > 0, a slighly more complicated, but still completely transparent picture emerges. Namely, one assumes a J = 1, a J+1 = δ 1. The self-consistent ansatz to substitute into (9) is then taken to be a k = q k δ |k−J| . With this ansatz, equation number J determines β 0 + J(β 1 − β 0 ). Equations number J + 1 and J − 1 thereafter determine β 0 , β 1 and q J−1 . These equation are effectively quadratic and have two solutions, which explains the two branches seen in [6,14]. Proceeding further iteratively, each pair of equations number J + n and J − n with n ≤ J forms a 2 × 2 linear system determining q J+n and q J−n . Finally, each equation with number greater than 2J determines the q with the corresponding number, much like what (11) does for the J = 0 case. We shall see below how this procedure works explicitly for J = 1.
We briefly comment on the general nature of the small δ expansion, of which we are mainly interested in the leading terms. (Systematic investigations of this issue are outside the scope of our present study.) The pattern we have displayed above suggests an asymptotic expansion around δ = 0 of the form Indeed, after the leading order contribution has been determined in the way we have described, subleading terms can be fixed order-by-order by substituting (12) into (9) and solving the resulting linear equations. Importantly, the prefactor δ |n−J| completely drops out when substituted in the equations (the leading order is δ-independent, as seen, for example, in (11) for J = 0). For that reason, the fact that the asymptotic expansions for different a n start with different powers of δ by no means upsets the consistency of determining subleading corrections simultaneously for all a n . Note that, for corrections of each given order in δ, the infinite sums in (9) will truncate to a finite number of terms. Indeed, in (9) written for mode number n, the right-hand side contains an infinite sum over j and k, with terms bearing factors of δ |j−J|+|k−J|+|n+j−k−J| . Since all terms in the power of δ are positive and the first two terms grow as j or k are increased for j, k > J, it is evident that, for each given n and J, only a finite number of choices of j and k exist that keep the power of δ below any given value. This leaves a finite number of terms out of the sum in (9) at each n and at each order in δ.
Asymptotic expansions of the form (12) strongly suggest that the leading order results we deal with here become as accurate as one needs for sufficiently small δ. The corrections in (12) are uniquely determined by solving linear equations and simply slightly shift the existing leading order spectral branches we have described above.

Asymptotics of the interaction coefficients
We shall now analyze the asymptotic behavior of the C coefficients appearing in (11) for AdS 4 . The mode functions for global AdS d+1 are given by (4). For large mode numbers, we can use the approximation (see, e.g., [18]): to write the mode functions as For generic dimension d, the approximate expression for the mode functions will diverge near x = 0, and we will not be able to employ the straightforward approach for extracting the asymptotics of C coefficients that we utilize below for d = 3. At d = 3, does not diverge. We pause to note that this formula is closely related to, and could also be easily derived from, the exact simplified expressions for the AdS 4 mode functions that appeared in [16] while we were preparing the current article for publication. Those expressions were used in [16] to obtain closed-form (albeit very long) expressions for the interaction coeffients in the gravitational case, and this could also be done for the case of the probe scalar. While the availability of such closed-form expressions for the interaction coefficients in this particular dimension is certainly a welcome development, the simpler analysis of the asymptotics provided here will be sufficient for our purposes in the present paper. In passing, we also note that an efficient recursive way to compute interaction coefficients in arbitrary dimension is provided in [17].
We shall start by displaying the asymptotic formula for C with four large indices, which is an instructive aside: In a general number of dimensions, we expect the first powers of the combinations of frequencies in this expression to generalize to power d − 2, in analogy to the considerations of [17]. The actual coefficients appearing in (11) have at least one small index, while the other indices may become large asymptotically. For three large indices, we have 2 and C nmk0 32 if max{n, k, l} > (n + k + l + 2)/2 π/16 if max{n, k, l} = (n + k + l + 2)/2 3π/16 if max{n, k, l} = (n + k + l + 1)/2 π/4 otherwise (18) In particular, We proceed with two small and two large indices, using if |m − n| = 3 7π/64 if |m − n| = 2 5π/32 if |m − n| = 1 or |m − n| = 0 0 otherwise (21) 2 Here and below, we are using the explicit formulas P  So in particular C 0,1,n,n+1 → 20 So in particular C 0,0,n,n → 12 π ≈ 3.8197. These observations for large mode number limits taken in different directions in the mode number space reveal a very simple plateau-like topography, see Fig. 1. Almost everywhere in the large index value region, with one index kept equal to 0 according to (11), the plot of C is completely flat at a fixed height of 8 √ 2/π. If one goes to the border regions of this plateau, where one of the three large indices becomes small, the plot rises somewhat, but neither is the rise numerically significant (only a few percent), nor is it statistically significant when appearing inside the sums, as in (11), since this rise is only attained for special, non-generic values of the indices. In the end, the large mode-number asymptotics of the solutions to (11) can be characterized almost entirely in terms of the single 'plateau' value of C, i.e., 8 √ 2/π. We shall see below how this works explicitly.

Solutions dominated by mode 0: asymptotic estimates
Armed with asymptotic expressions for the C coefficients, we can start our analysis of the large mode-number limit of the solutions to (11). As we already explained, one assumes in (9) a 0 = 1, a 1 = δ 1 and a n = q n δ n . From the equations number 0 and 1 in (9), one obtains .
For d = 3, this becomes At higher orders in δ, one has to solve (11). Note that the C coefficients inside the sum in (11) are nearly constant at large k, except for the endpoints of the sum, where they vary slightly. Since there are considerably more terms in the interior of the sum than near the endpoints, it is natural to assume that the equation can be approximated well by simply taking C = 8 √ 2/π (the plateau value we have identified above), solving the equation and then checking that our assumption has been justified. It is likewise natural to assume that the term proportional to k 2 in the coefficient of q k dominates the lower powers of k. All these assumptions can, of course, be explicitly checked by an a posteriori comparison to numerical evaluations of the solutions to (11).
Applying the above simplifications to (11), one obtains: One can assume the simple asymptotic behavior q n = Qn γ e κn , and check for what values of the parameters κ, γ, Q it is compatible with (25). We find that κ completely drops out from (25). It therefore cannot be determined by simple asymptotic matching and requires solving (11) iteratively mode-by-mode. After κ has dropped out, one is left with For large k, we can approximate the sum by an integral over ζ = j/k: Matching the leading powers on both sides of the equation fixes γ = 1. Thereafter, (26) fixes Q = 3 √ 2/10. We thus expect the asymptotic spectrum to behave as a n 3 √ 2 10 n (δe κ ) n .
Note that e κ cannot be absorbed in δ since the normalization of δ is fixed at the infrared end of the spectrum by imposing a 1 = δ. We have checked numerically that (28) accurately approximates exact solutions to (11). In particular, we have verified that no logarithmic modulations are present in addition to the power-law modulation in (28); we will explain the method we used in section 3.5. The situation will not be similarly simple when we turn to the gravitationally interacting case.

Solutions dominated by mode 0: systematic analysis
We have explained above that natural simplifying assumptions (which can be verified a posteriori by direct comparison with numerics, even in the absence of rigorous error bounds) reduce the quasiperiodicity conditions (11) for solutions dominated by mode 0 to (25). We have then verified that the asymptotics (28) is analytically consistent with the simplified equation (25), and also provides a valid approximation to the full numerical solutions of (11). At a practical level, one could rest content with (28). However, it would be better to have a more explicit analytic control over how this asymptotics emerges. (Note that being asymptotically consistent with an equation does not, strictly speaking, imply that solutions to this equation, with initial data specified in the infrared, approach asymptotics of this form.) Equation (25) provides a relatively simple iterative definition of q n via a quadratic recursive relation reminiscent of the one for the Catalan numbers. It can be analyzed quite exhaustively, with a curious mathematical structure emerging, which is what we shall explain next.
To gain further analytic insights in (25), we introduce a generating function (25), together with the initial conditions q 0 = q 1 = 1, is then equivalent to This equation has an obvious mechanical interpretation: a particle is moving in the q-space and λ-time with an external force e −λ and a potential −10 √ 2(q − 1) 3 /3. One is looking for a special solution of this equation that stops at λ = +∞ at the inflection point of the potential, i.e., q = 1. One-dimensional mechanics in a cubic potential is well-known to display a run-away behavior in finite time (all solutions run down the potential slope to q = +∞ at a certain finite moment of time). In the vicinity of that run-away moment, which we shall call κ, the e −λ term in (30) is negligible and one can integrate the equation to obtain Of course, the particular value of the run-away moment κ is sensitive to the entire evolution and requires solving (30) exactly, which cannot be done analytically.
What have we learned? q defined by (29) is an analytic function of z = e −λ with a double pole singularity at z = e −κ . q n are the coefficients of the Taylor series expansion of this function around z = 0. It is generally understood that high order behavior of the Taylor series is governed by the singularities of the expanded function in the complex plane. Darboux' theorem (see, e.g., [19]) gives a precise statement of this correspondence. In particular, a double pole at z = e −κ precisely corresponds to q n ∼ ne κn . Furthermore, the coefficient of (31) matches the coefficient of (28). Overall, our systematic analysis of (25) based on generating functions completely corroborates the simple asymptotic estimates of the previous section.

Solutions dominated by mode 1
As a proof-of-concept demonstration, we shall now describe how to analyze solutions dominated by mode 1. Compared to the solutions dominated by mode 0 considered above, this case presents the feature of having two spectral branches, originally seen for gravity in the numerical simulations of [6]. Extension to solutions dominated by still higher modes, on the other hand, is a pure technicality.
We assume a spectrum a 1 = 1, a 2 = δ 1 and a n = q n δ |n−1| , substitute it in (9) and retain only the leading powers of δ in each equation. (9) with n = 0, 1, 2 gives For d = 3, β 1 = 19/18π. Solving the two remaining equations simultaneously leads to with two solutions, corresponding to two spectral branches of the type seen numerically in [6]: .
The remaining equations reduce to In the limit of large k, we found that, for d = 3, the interaction coefficients are approximately constant: By evoking similar logic to what we have already used in our mode 0 derivations (without recapitulating all the subtleties), we assume an asymptotic ansatz q n = Qn γ e κn and substitute it in (38), obtaining Matching the leading power of k on both sides of this equation results in γ = 1. We also obtain the following relation between Q and κ, with the specific numerical values for β 0 , β 1 , q 0 , C 0 , C 1 given above. We then solve (38) numerically and fit the results to q n = Qn γ e κn (1 + c n + c 2 n 2 ). This is done by fitting 5 adjacent points to the above expression, and then moving these 5 points to higher and higher values of n in hope to see convergence. Quick convergence does, in fact, follow. For the solution (36), we get γ = 1, κ ≈ −0.06, Q ≈ 0.0846. For the solution (37), we get γ = 1, κ ≈ −0.648, Q ≈ 1.75. Both sets of values of κ and Q are in agreement with (41). We therefore conclude that our analytic treatment of the large mode-number limit provides a sound perspective. Attempts to augment the asymptotics with logarithmic modulations, for example, a factor (ln n) γ 1 , and obtaining estimates for γ 1 , result in a quick convergence of γ 1 to 0 as the mode number increases. This indicates the absence of logarithmic modulations, in accordance with our analytic results. We shall now turn to the case of gravitational interactions, where logarithmic modulations are, in fact, present.

Gravitational interactions 4.1 Gravitational quasiperiodic solutions
In the introduction, we have briefly explained that, compared to the self-interactive probe scalar field we have considered above, the effective time-averaged description for the full gravitational non-linearities is conceptually quite similar, but technically much more complicated. The averaged equations, as derived [8], are of the form (we are now using real amplitudes A n and phases B n , related to the complex amplitudes by A n e iBn ≡ α n ): The coefficients of the terms appearing in these time-averaged equations are expressed in terms of certain integrals of products of the AdS mode functions (4): The integrals in the expressions above are defined by x 0 dy e k (y)e l (y)µ(y), x 0 dy e k (y)e l (y)µ(y).
The construction of quasiperiodic solutions [6] of (42-43) is implemented analogously to the procedure we have described for (9). One first assumes B n = (β 0 + n(β 1 − β 0 ))ε 2 t, whereupon (42) ensures that A n are constant (this is how quasiperiodic solutions of the sort treated here are defined). Substituting the ansatz for B n into (43), on the other hand, leads to an infinite system of algebraic equations for β 0 , β 1 , A n , analogous to (9).
One can furthermore concentrate on solutions strongly localized around one chosen mode, as we did in section 3.1, with the localization controlled by an adjustible parameter δ. The resulting structure is straightforwardly analogous to (11), with the small δ limit truncating the infinite sums in quasiperiodicity condition to finite sums, and solution for β 0 , β 1 and (iteratively) A n becoming possible. In particular, the same explanation as in section 3.1 for the existence of two spectral branches holds. (The branches, present for solutions dominated by modes other than mode 0, were originally observed numerically in [6].)

Asymptotics of the interaction coefficients
As for the probe scalar case, the asymptotic analysis of the quasiperiodicity condition for gravity involves understanding the properties of the interaction coefficients appearing in (42-43). Of particular importance are the S coefficients with two or three large indices (since these are the coefficients controlling the complicated sums coupling a large number of different triplets of modes).
For d = 3 one can apply a similar sort of analysis to what we described in section 3.2 to extract the values of the X and Y integrals when three indices becomes large. For k > 0 and 0 j k, one gets With these values, one can compute the corresponding S coefficients. One discovers, however, that the dominant contribution, which would scale like the square of the index numbers according to (46), in fact exactly vanishes. This leaves S growing as the first power of the index number (under uniform scaling of the indices) when 3 indices become large, and when the resonance condition is satisfied. One can perform similar analysis for S with only two indices becoming large. In this case, cancellations do not occur and one discovers that S grows like the index number squared (under uniform scaling of the two large indices).
The above considerations leave one with a topography of the plot of the S coefficients (see Fig. 2) much more complicated than the simple plateau we have observed asymptotically for the case of the probe scalar field. For example, for S 0kj,k−j , if both indices k and j are scaled up uniformly in proportion to λ, the S-coefficient also scales in proportion to λ. But if j (or k − j) is allowed to stray close to small values, the plot rises dramatically, attaining values proportional to k 2 .

Solutions dominated by mode 0
We can now describe the analytic structure of gravitational quasiperiodic solutions strongly dominated by mode number 0. Just as we did for the probe scalar field, one assumes a spectrum of the form A 0 = 1, A 1 = δ 1 and A n = q n δ n and substitutes it in (43), together with B n = (β 0 + n(β 1 − β 0 ))ε 2 t, keeping only the lowest powers of δ in each equation. Equations number 0 and 1 then determine β 0 and β 1 . The remaining equations take the form One cannot, however, proceed with this equation in the same way we did for the probe scalar field case. The reason is that the dependence of S on its indices favors the summation index regions close to the end-points in (54), where q itself cannot be approximated by its asymptotics. At this stage, we cannot see any obvious ways to proceed with our asymptotic analysis and extract the large mode number behavior of q, even if one has detailed control over the asymptotics of S and R (for example, via the exact analytic expressions of [16]). One might be surprised that the exact asymptotics for the gravitational case proves so elusive, given that simple fitting formulas, q n ∼ e κn /n in our language, have been proposed for solutions dominated by mode 0 in [6]. (Note that, just like for the scalar case, e κn drops out if substituted in (54) and it simply corrects δ n to the actual asymptotic exponential fall-off, given that the normalization of δ is fixed in the infrared by imposing δ = A 1 .) If the asymptotic solution is as simple as q n ∼ e κn /n, it should seem strange that one cannot at least verify it using the equations and encounters all the analytic difficulties we have observed. Given this perception, we have decided to re-examine the issue numerically. Even in the absence of successful asymptotic analysis, (54) provides a straightforward iterative way to evaluate q n , which can be easily programmed in Mathematica. We have implemented this procedure up to mode number 70 in hope of getting a better sense of the asymptotic behavior of q n .
We first attempted to fit our numerical findings to q n = Qn γ e κn , which is what we have used for the scalar case. All of our fits are performed by choosing a group of adjacent points with the number of points equal to the number of unknowns in the putative asymptotic formula, fitting exactly to this group of points, and then moving the group of points to higher and higher index values in search of convergence. For q n = Qn γ e κn , we observed very slow convergence for γ. While not being very far from −1, γ fails to reach a plateau all the way until mode number 70.
We have then attempted the same procedure with a logarithmic modulation q n = Qn γ (ln n) γ 1 e κn and observed that the convergence of γ is considerably improved, while γ 1 does not appear to go anywhere close to 0. (Note that a similar procedure applied to the probe scalar field spectrum yields γ 1 rapidly converging to 0.) We finally attempted to test the presence of ln(ln n) corrections, though admittedly our mode range is not particularly large to allow tracing very slow modulations reliably. Nonetheless, applying the fitting procedure to q n = Qn γ (ln n) γ 1 (ln ln n) γ 2 e κn 1 + c 1 n + c 1 n ln n + c 2 n 2 + c 2 n 2 ln n + c 2 n 2 ln 2 n (55) results in good convergence within the examined mode range. γ rapidly reaches −1, within a percent, while γ 1 and γ 2 appear to be around 0.3 and 0.1, respectively. (These values change by only a few percent between estimating around mode 50 and around mode 70 with our procedure, indicating adequate convergence. Removing the corrections decreasing as 1/n 2 and faster from the above ansatz makes the plots less flat, though the actual estimated values for the powers remain close to the ones quoted above.) We thus find evidence that quasiperiodic spectra of AdS perturbations with gravitational non-linearities display slow logarithmic modulations in their ultraviolet asymptotics, making them considerably more challenging for both analytic and numerical treatments than the power laws previously mentioned in the literature. It would be interesting to analyze such logarithmic corrections further, especially with the explicit expressions for the interaction coefficients in AdS 4 that have just appeared in [16]. The numerical optimization allowed by these expressions makes it possible to go to much higher mode numbers than what one can get by brute force evaluation. (See also our parallel work [17] that presents effective iterative techniques for evaluating the interaction coefficients in a general number of dimensions.)