Conformal flow on $S^3$ and weak field integrability in AdS$_4$

We consider the conformally invariant cubic wave equation on the Einstein cylinder $\mathbb{R} \times \mathbb{S}^3$ for small rotationally symmetric initial data. This simple equation captures many key challenges of nonlinear wave dynamics in confining geometries, while a conformal transformation relates it to a self-interacting conformally coupled scalar in four-dimensional anti-de Sitter spacetime (AdS$_4$) and connects it to various questions of AdS stability. We construct an effective infinite-dimensional time-averaged dynamical system accurately approximating the original equation in the weak field regime. It turns out that this effective system, which we call the conformal flow, exhibits some remarkable features, such as low-dimensional invariant subspaces, a wealth of stationary states (for which energy does not flow between the modes), as well as solutions with nontrivial exactly periodic energy flows. Based on these observations and close parallels to the cubic Szego equation, which was shown by Gerard and Grellier to be Lax-integrable, it is tempting to conjecture that the conformal flow and the corresponding weak field dynamics in AdS$_4$ are integrable as well.


Introduction
Propagation of nonlinear waves in confining geometries presents significant challenges because the key mechanism stabilizing the evolution of waves on unbounded domains, which is the dispersion of energy by radiation, is missing in confined settings. Consequently, an arbitrarily small perturbation of a ground state can lead to complicated long-time behavior. The central physical problem in this context is that of energy transfer, namely, understanding how the energy injected into the system gets distributed over the degrees of freedom in the course of its evolution. This problem has been studied in the physics literature from a statistical viewpoint under the name of wave (weak) turbulence [1,2], but only recently emerged as an active research topic in dispersive PDE theory. In this deterministic mathematical approach, the energy cascades from low to high modes are quantified by the growth of higher Sobolev norms of solutions and the main question is whether these norms can become unbounded in finite or infinite time. The past few years have witnessed significant progress in the understanding of this issue in the context of nonlinear Schödinger equations (and their variations) on tori [3][4][5]. However, almost nothing is known, as far as we can tell, about manifestations of these phenomena in other interesting evolution equations of mathematical physics.
Our attention in this article will be focused on the weak field dynamics for the conformally invariant cubic wave equation on the three-sphere. This very natural geometrical model describes relativistic scalar waves propagating on a compact manifold and interacting with the background geometry and among themselves. The spherical geometry prevents the waves from dispersing to infinity and thus ensures that the nonlinear self-interactions remain important for all times, inducing complicated energy transfer patterns. Another important circumstance is that all normal mode frequencies of the corresponding linearized theory are integer in appropriate units. Such a fully resonant spectrum ensures that nonlinearities can produce significant effects over long times for arbitrarily small perturbation amplitudes, generating highly nontrivial weak field dynamics.
The weak field dynamics of the system we are considering is closely related to investigations of nonlinear stability of anti-de Sitter (AdS) spacetime initiated in [6] (see [7,8] for brief reviews and further references). A conformal transformation relates our considerations to the dynamics of the cubic Klein-Gordon equation with mass m 2 = −2 (conformally coupled self-interacting probe scalar) on the four-dimensional AdS space. In contrast to the bulk of AdS stability research, we do not consider the backreaction of perturbations on the AdS geometry. However, a number of features of the weak field dynamics remain unaffected by our simplification. (The weak field dynamics of a massless self-interacting probe scalar in AdS was introduced as a toy model for the AdS stability problem in [9].) Our main tool for analyzing the weak field dynamics, given a fully resonant spectrum of the linearized system, is the time averaging method that goes back to Bogoliubov and Krylov and is described in many books on perturbation theory (e.g., [10]). For a contemporary treatment adapted to studies of nonlinear PDEs, see [11]. Applied to the PDE we are considering, the time averaging method produces a simplified infinite-dimensional dynamical system, which we shall refer to as the conformal flow. By standard theorems underlying time averaging, this effective system (in PDE theory sometimes called the resonant system) accurately approximates the dynamics of the original PDE for small amplitude fields of order O(ε) on long time-scales of order O(1/ε 2 ), and in particular it accurately captures the energy transfer patterns on such time scales. (Time averaging and related techniques were introduced to the AdS stability problem closely related to our current studies in [12][13][14].) The conformal flow displays a number of highly special features not apparently present in the original PDE and suggestive of integrability. This includes additional conserved quantities, low-dimensional invariant subspaces, a variety of stationary states without any energy transfer, as well as solutions with exactly periodic energy transfer patterns.
Our original PDE can be viewed as an infinite system of oscillators with integer frequencies and a quartic potential. We do not see any indications that it might be integrable; only the time-averaged system describing its weak field dynamics displays signs of integrability. While extensive searches for quartic integrable systems (i.e., integrable mechanical analogs of our original PDE) have been conducted in the past (see [15] for a review and references to earlier work and [16] for some later contributions), we are not aware of any broad scans for systems displaying weak field integrability of the type suggested by our results.
The conformal flow introduced in this paper bears a strong resemblance to the cubic Szegő equation that has been designed and studied in a series of papers by Gérard and Grellier, of which [17][18][19][20] are particularly relevant for our purposes. This equation also emerges as the resonant approximation of a nonlocal transport equation, called the cubic half-wave equation [18], which is a special case of one-dimensional models of wave turbulence introduced in [21]. In [17][18][19][20] and related publications, a number of remarkable properties were established for the cubic Szegő equation, in particular, a Lax pair structure, existence of finite-dimensional invariant subspaces, as well as some weak turbulence phenomena. The bulk of our effort will be directed at recovering analogs of a subset of these results for the conformal flow, a system with a natural geometric origin. We furthermore believe that the analogy goes deeper than what we are able to explicitly demonstrate in this article.
There are also close similarities between the conformal flow and resonant approximations for the cubic Schrödinger equation, either with harmonic trapping [22] or in the infinite volume limit on the 2-dimensional torus [23,24]. Another system which shares some properties with the conformal flow and the cubic Szegő equation is the Lowest Landau Level evolution equation which appears in the studies of rapidly rotating Bose-Einstein condensates (see [25] and references therein). We furthermore mention the curious case of the Charney-Hasegawa-Mima (CHM) equation explored in [11]. Also in that case, effective integrable dynamics emerges in the time-averaged approximation (which the authors call the 'effective equation'). However, a substantial difference is that the linearized spectrum of the CHM equation is not fully resonant and splits into disconnected resonant clusters. In such a situation, the resulting effective time-averaged system literally separates into independent finite-dimensional integrable subsystems. This simplification does not happen in our case because our linearized spectrum is fully resonant.
We comment on how our work is related to extensive investigations of the last decade in the area of 'AdS/CFT Integrability' (see [26] for a review). One of the central ingredients of that line of research is integrability of sigma models (describing string worldsheets) on target spaces involving AdS factors. While we are also presently talking about integrability emerging due to special features of the AdS geometry, there are also essential differences. We are considering a field system evolving in AdS, rather than a sigma model (whose fields are maps from a two-dimensional surface into AdS). We are furthermore not talking about exact integrability of our AdS field system, but rather about the integrability of the effec-tive weak field theory (conformal flow) emerging from it within the standard time-averaged approximation.
Our exposition is organized as follows. In section 2, we formulate our scalar field dynamics on a three-sphere, describe its connection to AdS stability problems and apply time averaging to derive the corresponding effective system, which we name the conformal flow. In section 3, we give a very elementary pragmatic summary of some properties of the cubic Szegő equation relevant in our context, aimed mostly at the physics audience. In section 4, we construct a three-dimensional invariant subspace of the conformal flow and discuss some complex plane properties of the generating functions that encode the conformal flow amplitudes. In section 5, we analyze stationary states of the conformal flow, in which the initial conditions are adjusted to ensure that no energy transfer occurs. In section 6, we exhibit explicit solutions with periodic energy flows on the three-dimensional invariant subspace of the conformal flow and discuss their properties. Finally, in section 7, we give a summary of what we have practically demonstrated and an outlook in more ambitious directions, keeping the likelihood of Lax-integrability in mind.
2 Conformally coupled scalar field on the Einstein cylinder and in AdS 4 and its weak field dynamics As a simple model of confining geometry we consider the Einstein cylinder, which is the globally hyperbolic spacetime (M, g) with topology R × S 3 and metric where x ∈ [0, π], dω 2 is the round metric on the unit 2-sphere and r is the radius of S 3 . This spacetime has a constant scalar curvature R(g) = 6/r 2 . As a model of nonlinear dynamics on M, we consider the semilinear wave equation for the real scalar field φ : where g := g µν ∇ µ ∇ ν is the wave operator associated with g and λ is a constant. For concreteness, we assume that λ > 0 (which corresponds to a defocusing nonlinearity); however, in the small data regime that we focus upon here, all the results below hold true in the focusing case λ < 0 as well. Due to the identity equation (2) is conformally invariant. After rescaling t → t/r and φ → r √ λ φ, it takes the dimensionless form which can be interpreted as the cubic Klein-Gordon equation on the unit S 3 with a unit mass. For simplicity, we shall restrict our analysis to rotationally symmetric fields, which depend only on (t, x). Substituting φ(t, x) = v(t, x)/ sin(x) into (4) we obtain a nonlinear string equation with Dirichlet boundary conditions v(t, 0) = v(t, π) = 0, which are enforced by regularity of φ(t, x) on M.
Remark 1. We could have arrived at the same equation considering the conformally coupled self-interacting scalar (2) in the 4-dimensional anti-de Sitter (AdS 4 ) spacetime with the metricg = 1 where x varies between 0 and π/2. Indeed, g = cos 2 xg, hence it follows from (3) that the field redefinitionφ(t, x) = φ(t, x) cos x converts equation (2) from the Einstein cylinder to AdS. In order to define the evolution, one has to impose a boundary condition at x = π/2. Here, we shall impose the Dirichlet condition φ(π/2) = 0 on the equator, corresponding to the reflecting boundary conditions most commonly used in the AdS research context. This model is just a subsector of the model on the full Einstein cylinder, which can be implemented by imposing the reflection symmetry φ(t, x) = −φ(t, π − x). We emphasize this connection because problems involving nonlinear dynamics of small AdS perturbations have received a significant amount of attention in recent years, and our objectives here have much in common with that body of work.
Our goal is to understand the evolution of small smooth initial data v(0, x), ∂ t v(0, x). Decomposing the solution into a Fourier series, we get from (5) an infinite system of coupled oscillators with the interaction coefficients To factor out fast linear oscillations in (8), we change the variables using variation of constants This transforms the system (8) into where each c j in the sum is given by (10). Thus, each term in the sum has a factor e −iΩt , where Ω = (n+1)±(j +1)±(k +1)±(l +1), with all the three plus-minus signs independent. The terms with Ω = 0 correspond to resonant interactions, while those with Ω = 0 are nonresonant. Passing to slow time τ = ε 2 t and rescaling β n (t) = εα n (τ ), we see that for ε going to zero the non-resonant terms ∝ e −iΩτ /ε 2 are highly oscillatory, and therefore expected to be negligible (in a sense we are about to specify). Keeping only the resonant terms in (12) (which is equivalent to time-averaging), we obtain an infinite autonomous dynamical system where the summation runs over the set of indices {j, k, l} for which Ω = 0. Evaluating the integrals for the S-coefficients, one finds that this set reduces to {jkl | j + k − l = n} and for such resonant combinations of indices, Note that terms with n = j + k + l + 2 could have been present in principle, but the corresponding S-coefficients vanish, which can be verified by direct calculation. This is directly parallel to the selection rules that have been extensively discussed in the AdS stability literature [13,14,27,28]. It is part of the standard lore in nonlinear perturbation theory that solutions to (12) starting from small initial data of size O(ε) are well approximated by solutions of (13) on timescales of order O(ε −2 ). More precisely, if β n (t) and α n (τ ) are solutions of (12) and (13), respectively, and β n (0) = εα n (0) for each n, then |β n (t) − εα n (τ (t))| O(ε 2 ) for t O(ε −2 ). In other words, on this timescale the dynamics of solutions of equation (5) is dominated by resonant interactions. Straightforward proofs (normally phrased for systems with a finite number of degrees of freedom) can be found in textbooks on nonlinear perturbation theory, e.g., [10].
For convenience, we shall henceforth rescale τ to remove the numerical factors 2 and −3 from (13). Renaming so rescaled τ back to t and using · = d/dt, we finally arrive at the following system, which we call the conformal flow Remark 2. Note that the AdS formulation described in Remark 1 simply corresponds to setting all even-numbered modes to 0, while only keeping odd-numbered modes α 2m+1 with m running from 0 to infinity. It is easy to see that if one implements this constraint in (15) and expresses everything through the new mode counting index m, one gets back equation (15). Thus, the conformal flow is equally applicable to the dynamics in R × S 3 and in AdS 4 .
The conformal flow (15) is Hamiltonian with and symplectic form n 2i(n + 1) dᾱ n ∧ dα n : It enjoys the following one-parameter groups of symmetries (where λ, θ are real parameters): The latter two symmetries respect the Hamiltonian structure and give rise to two conserved quantities (in addition to the Hamiltonian itself): 'Linear energy:' We emphasize that these two conservation laws have no obvious counterparts in the original wave equation (4). Their prototypes were presented in the AdS stability context in [9,14].

An elementary introduction to the cubic Szegő equation
The main content of our current treatment will be in presenting a range of remarkably simple dynamical behaviors of the conformal flow system (15). In anticipation of these results, it is useful to contemplate for a moment the apparent special features of equation (15) itself. One can trivially rewrite (15) as One notable property of this representation is that many coefficients of the trilinear form on the right hand side are simply 1. Indeed, this happens for the contribution of any triplet of modes with all frequencies higher than the recepient mode n. It turns out useful to consider the following simpler equation, in which all numerical coefficients are set to 1, We shall see that, while the systems (23) and (24) are distinct, their dynamics is qualitatively similar, with a number of parallels between the emerging algebraic structures. Actually, the simplified system (24) is the Fourier representation of a simple paradifferential equation that has been studied before. To show this, let u(t, e iθ ) be a complex function on the circle whose nonnegative Fourier coefficients are given by α n (t) and negative Fourier coefficients vanish, i.e.
and assume that ∞ n=0 |α n | 2 < ∞. Mathematically, this means that u belongs to the Hardy space on the circle L 2 + (S 1 ) ⊂ L 2 (S 1 ) [29]. Furthermore, we define an orthogonal projector Π : L 2 (S 1 ) → L 2 + (S 1 ) (called the Szegő projector) that filters out negative frequencies Then, it is easy to see that the system (24) is equivalent to This equation was introduced by Gérard and Grellier in [17] under the name of the cubic Szegő equation. They showed that this equation has a remarkably deep structure, including a Lax pair, finite-dimensional invariant subspaces and weakly turbulent solutions. We shall demonstrate below that the conformal flow exhibits at least some of this structure (in a modified form), and conjecture that the similarity goes even further. In view of these parallels, we shall present below a very elementary summary of the features of the cubic Szegő equation that are important to us. Readers interested in more details are referred to the original publications [17][18][19][20].
The cubic Szegő equation is Hamiltonian with and symplectic form n 2i dᾱ n ∧dα n , and has the same symmetries (18)(19)(20) as (15). Hence, it possesses two conserved quantities analogous to (21) and (22), in addition to the Hamiltonian itself: Borrowing terminology from studies of the nonlinear Schrödinger equation, we shall refer to these charges as 'mass' and 'momentum,' respectively. Note that M + P = u 2 H 1/2 , where H s ∈ L 2 + is the Sobolev space on the circle equipped with the norm This fact can be used to control the growth in time of higher Sobolev norms and thereby prove that the Cauchy problem for the Szegő equation is globally well posed for smooth initial data in H s with s > 1/2 [17]. Any function u(e iθ ) of the form (25) in L 2 + (S 1 ) can be identified with a holomorphic function u(z) inside the unit disk |z| < 1 in the complex plane [29], hence we can write In terms of u(t, z), the cubic Szegő equation (27) takes the form General properties of the Cauchy-type integral appearing in (33) suggest the existence of meromorphic solutions for u with a finite number of isolated poles moving in the complex plane outside the unit disk. Indeed, it was shown in [17] that the cubic Szegő equation admits finite-dimensional invariant subspaces given by rational functions of arbitrarily high degree. Moreover, due to the Lax-integrability (and, consequently, infinitely many conserved quantities), the Szegő flow restricted to such subspaces is completely integrable in the sense of Liouville. The simplest invariant manifold with nontrivial dynamics is parametrized by the following single pole ansatz We shall now demonstrate explicitly by elementary means that the subspace given by such functions is indeed dynamically invariant under the flow, and analyze its evolution. The Fourier coefficients corresponding to (35) are α 0 = b and α n = (a + bp)p n−1 for n ≥ 1.
Expressed through the parameters of our ansatz, the mass and momentum take the form Substitution of (35) into (33) produces exactly three distinct dependences on z on both sides, namely, (1 − pz) −2 , (1 − pz) −1 and a z-independent term, resulting in the following three equations for a(t), b(t), and p(t) where we used the conservations laws (37) to simplify some terms. These equations describe a three-dimensional Liouville-integrable system which can be easily solved (in [19] an explicit formula for the general solution was derived using harmonic analysis tools). One first integrates (38) to find a, whereupon the remaining equations become linear. The solution corresponding to initial conditions a(0) = a 0 , b(0) = b 0 , p(0) = p 0 (for simplicity assumed to be real) is where ω = 1 2 (M + P ) 2 − 4P a 2 0 and the expressions for M and P in terms of (a 0 , b 0 , p 0 ) are given in (37). Thus, the solution is quasiperiodic for all initial conditions, while the corresponding 'mass' spectrum, is expressible through |p| 2 and the conserved quantities, and hence exactly periodic. In the special case M = P = a 2 0 , we have ω = 0 and then The evolution of two-mode initial data a 0 = 1, b 0 = 2ε, p 0 = 0 was used in [17] to illustrate a very interesting instability phenomenon. For these data, (43) gives hence, in the regime of ε → 0, |p(t n )| ∼ 1 − ε 2 /2 for a sequence of times t n = nπ 2ω , meaning that the energy goes periodically to arbitrarily high modes. Put differently, although all solutions are quasiperiodic, their radius of analyticity is not uniformly bounded from below. Note that (46) implies instability of the one-mode stationary state u(t, z) = a 0 e −iM t z.
As mentioned above, one can consider more general meromorphic solutions with an arbitrary number of time-dependent simple poles outside the unit disk One can convince oneself that the ansatz is consistent in the sense that substituting it in (33) results only in z-dependences of the form (1 − p k z) −1 and (1 − p k z) −2 on both sides, producing 2N ordinary differential equations for the 2N functions b k and p k . It is possible to consistently impose p N = 0, producing an ansatz with one parameter less. The single pole ansatz (35) is precisely of such a form. The equations of motion within the subspace defined by (47) were written down explicitly in [17]. As in the simplest case of the three-dimensional invariant subspace discussed above, the dynamics on all finite-dimensional invariant subspaces is Liouville-integrable and bounded, and hence quasiperiodic.
Among the quasiperiodic solutions there exist special solutions with time-independent amplitudes |α n |. Such solutions, which we call stationary states, were classified in [17] in the case of finite dimensional invariant manifolds. They are given by either finite Blaschke products or where N ≥ 1 and ≤ N − 1 are nonegative integers.
It was recently shown in [20] that outside the finite-dimensional invariant subspaces there exists a dense set of smooth solutions whose radius of analyticity tends to zero for a sequence of times t n → ∞. Consequently, the Sobolev norms with s > 1/2 are unbounded. Such a weakly turbulent behavior is somewhat surprising in a completely integrable model. The coexistence of integrability and turbulence is possible because the infinitely many conserved quantities are too weak to control higher regularity properties of solutions.

Three-dimensional invariant subspaces of the conformal flow
We now return to the conformal flow (15) and attempt to treat it along the lines applied above to the cubic Szegő equation. While the treatment is closely parallel, we shall be able to immediately recover only a part of the results available for the cubic Szegő equation. Using the generating function u(t, z) and its conjugateũ(t, z) defined as in (32) and (34), we find that the system (15) is equivalent to the integro-differential equation The following summation formula has been used while obtaining this complex plane representation: Derivations of (50) and (51) and further comments are given in appendix A. Although equation (50) looks like a somewhat more elaborate version of (33), we do not see an immediate way to characterize its meromorphic solutions and defer it to future work. Nonetheless, the lowest-dimensional nontrivial invariant subspace is easy to construct, as we shall explicitly demonstrate now.
The ansatz relevant for the three-dimensional invariant subspace, analogous to (36), is where b, a, p are complex-valued functions of time. The corresponding generating function has the following combination of poles outside the unit disk We note that the two-mode initial data (α n = 0 for n ≥ 2) are accommodated within this ansatz as a special limiting case p → 0 with b and (a−b)p finite. All of our general statements about the dynamics on the three-dimensional invariant subspace apply to solutions starting with such initial configurations. While we could have used the complex plane representation (50) to establish the validity of our ansatz, it is instructive to apply brute force summations in this case. Substitution of (52) into (15) yields Note that the p n factor has consistently cancelled on the two sides. It remains to show that different n-dependences on both sides can be matched and produce a sufficiently small number of equations. All summations in (54) can be performed by applying the relation [min(n, j, k, n + j − k) which follows from (51). While we give explicit expressions for the sums involved in (54) in appendix B, the only thing one needs to know about these sums to establish the validity of our ansatz is that they are all polynomials of degree L in n, times (n + 1), the latter factor coming from Faulhaber's sums 1 + 2 l + · · · + n l , with l ≤ L, originating from 1 + θ + · · · + θ n in (55). Therefore, explicit counting tells us that, upon substituting the summation formulas of appendix B in (54), the factor (n + 1) will cancel on the two sides, leaving behind a statement that two quadratic polynomials in n equal each other. Matching the coefficients of these polynomials produces three ordinary differential equations for three functions b, a, p, confirming the validity of our ansatz. The explicit equations for b, a and p are given by iȧ (1 + y) 2 = a 6 5|b| 2 + (18y 2 + 4y)|a| 2 + (6y − 1)ba + 10yāb , iḃ (1 + y) 2 = b |b| 2 + (6y 2 + 2y)|a| 2 + 2ybā + a 2y|b| 2 + (4y + 2)y 2 |a| 2 + y 2b a , where we have introduced the following notation, which will turn out to be useful later: The conservation laws (21) and (22) take the form Q =(1 + y) 2 |b| 2 + 4y Re(ba) + 2y(3y + 1)|a| 2 , (60) E =(1 + y) 2 (1 + 2y)|b| 2 + 4y(3y + 2) Re(ba) + 4y(6y 2 + 6y + 1)|a| 2 .
We remark that if u(t, z) is a solution of the conformal flow, so is z N u(t, z N +1 ) for any nonegative integer N . This automatically generates an infinite number of other threedimensional invariant subspaces of the conformal flow, in which only subsets of modes are activated.
We shall return to equations (56-58) in section 6, and analyze the dynamics within the three-dimensional dynamically invariant subspace (52) explicitly. Before proceeding in that direction, we shall discuss special solutions for which |α n | are time-independent.

Stationary states
The conformal flow (15) admits solutions of the form where the frequencies λ n and complex amplitudes A n are time-independent, and λ n are linear in n, that is λ n = λ−nω for some real λ and ω. For such solutions there is no energy transfer between the modes, hence we call them stationary states. The cubic Szegő equation admits a variety of stationary states [17]; such solutions have also been considered in the context of AdS stability [12,30,31], where they were referred to as 'quasiperiodic' solutions. Substituting (62) into (15) we get a nonlinear 'eigenvalue' problem [min(n, j, k, n + j − k) The simplest solutions of this algebraic system, easily seen by inspection, A n = c δ N n , λ = |c| 2 , ω = 0 for c ∈ C, give the one-mode stationary states for any non-negative integer N α n = c δ N n e −i|c| 2 t .
We will see shortly that these trivial solutions are the endpoints of two-parameter families of stationary solutions. We note in passing that they are seeds for time-periodic solutions of the original equation (5) whose construction will be described elsewhere along the lines of [32]. Within the ansatz (52), stationary states take the form Plugging this into the system (56-58), we obtain a system of algebraic equations that can be solved explicitly. The case a(0) = 0 yields For this solution Q = |c| 2 (1−|p| 2 ) 2 , hence λ = Q. The corresponding generating function reads For nonzero a we get the following two-parameter family of stationary states with ω = 0 For this solution hence λ = 1 2 Q. The corresponding generating function reads In addition, for |p| ≤ p * := 2 − √ 3 ≈ 0.268, there is a pair of two-parameter families of stationary states with nonzero ω. The range of |p| is restricted by the condition that κ := |p| 4 − 14|p| 2 + 1 should be real, so we have a tangential bifurcation at p * . These two solutions (labelled by ±) are given by and their conserved quantities are The corresponding generating function reads where p(t) = p(0)e iω ± t . Outside the three-dimensional invariant subspace given by (52) there exist many other stationary states. For instance, we verified that, in close parallel to the cubic Szegő equation, any finite Blaschke product yields a stationary state. There are also stationary states where only every N th mode is activated, for instance Remark 3. For all the above stationary states with ω = 0 we have λ = Q N +1 , where N is the number of zeroes of the generating function (counted with multiplicity). It would be interesting to find a reason that underlies this 'quantization' structure.
It is well known that stationary states of Hamiltonian systems admit a variational characterization. In our case, it follows from the Hamilton equations (17) that stationary states (62) are the critical points of the functional hence stationary states with ω = 0 are extrema of H for fixed Q, while stationary states with nonzero ω are extrema of H for fixed Q and E. This fact is very helpful in determining stability properties of stationary states, as will be described elsewhere.

Dynamics on the three-dimensional invariant subspace
We now return to equations (56-58) describing the dynamics on the three-dimensional invariant subspace of the conformal flow and demonstrate how to solve them explicitly. Once again, the situation is closely parallel to the three-dimensional invariant subspace of the cubic Szegő equation with its underlying one-dimensional periodic motion. From (56), we obtain for y defined by (59) This equation and (57) imply that is conserved. This quantity is related to the conformal flow Hamiltonian (16) by Equations (60), (61) and (80) can be resolved to express |b| 2 , |a| 2 and Re(ba) through Q, E, S and y as follows From these relations, together with (79), one obtainṡ Note that this equation has the algebraic form of energy conservation for an ordinary onedimensional harmonic oscillator. The solution of (83) reads where B and A are constants depending on E, Q, and S; ψ is a phase that must be determined from the entire set of initial conditions. By equation (82), this exactly periodic motion is transferred to |b| 2 , |a| 2 and Re(ba), and hence to the mode energy spectrum |α n | 2 = |b + na| 2 |p| 2n = |b| 2 + 2n Re(ba) + n 2 |a| 2 |p| 2n .
One will thus observe exact returns of the energy spectrum to the initial configurations for all solutions within the three-dimensional invariant subspace of the conformal flow. The turning points of the periodic motion of y described by (84), y ± = B ± A, are given by the roots of the quadratic polynomial appearing on the right-hand side of (83) (the special case A = 0 corresponds to stationary states). In terms of the energy spectrum (85), these turning points provide lower and upper bounds for the inverse and direct cascades of energy, respectively. One of the key questions in this context is how large can y grow starting from some small y. To answer this question, we consider the ratio An advantage of the last representation is that it only contains combinations of y + and y − directly expressible through the coefficients of the quadratic polynomial on the right-hand side of (83). One thus gets We furthermore notice that Therefore, dividing the numerator and denominator of (87) by Q 2 and replacing the numerator by its maximum and denominator by its minimum, we obtain a simple uniform upper bound (which very likely can be tightened with extra work) This proves that the transfer of energy to high frequencies (or, equivalently, the growth of the higher Sobolev norms (31)) is uniformly bounded. We recall from section 3 that there is no such bound for the Szegő flow on the three-dimensional invariant subspace (35). In this sense, the conformal flow is much less 'turbulent' than the Szegő flow. Heuristically, this is not very surprising in view of the fact that the interaction coefficients in (23) min(n, j, k, n + j − k) + 1 n + 1 decay when the recipient mode number n is large and at least one of the source mode numbers is small, undermining the efficiency of energy transfer from low to high frequencies. This is in contrast to the Szegő flow where all coupling coefficients are equal to 1.

Summary and open questions
Starting with a naturally defined geometric PDE (4) on a three-sphere describing a selfinteracting conformally coupled scalar field, we have considered its effective time-averaged weak field dynamics, arriving at the conformal flow (15). The conformal flow is both structurally similar and displays a number of dynamical parallels to the previously known cubic Szegő equation (24). Building on these analogies, we have revealed the dynamics on the threedimensional invariant subspace of the conformal flow described by the ansatz (52). Within this subspace, the dynamics is Liouville-integrable and bounded, and hence quasiperiodic, with exactly periodic energy flows (85) displaying an alternating sequence of direct and inverse cascades. Unlike the cubic Szegő case, the conformal flow dynamics within the threedimensional invariant subspace does not display turbulent behaviors (not even the 'weak weak turbulence' of [3]). Other similar three-dimensional subspaces can be immediately constructed, as follows from simple properties of the conformal flow equation. There are furthermore direct parallels between the structure of stationary states of the Szegő and conformal flows (these states are special solutions for which no energy transfer between the modes occurs). Of these, the Blaschke product states (76) are particularly intriguing, as they clearly demonstrate that the parallels between the Szegő and conformal flows extend beyond the three-dimensional invariant subspaces that have been our main focus here.
We conclude with a list of open questions that we leave for future investigations: • Is the conformal flow (15) Lax-integrable?
• Is there a way to re-express the conformal flow through a projector similar to the Szegő projector? Apart from the conceptual importance, such a projector should make calculations easier. The complex plane representation of the conformal flow (50) can be seen as a first step in this direction.
• Are there higher dimensional invariant subspaces? The three-dimensional subspace we have considered, with its double-pole generating function (53), and parallels to the cubic Szegő case strongly suggest that more elaborate meromorphic ansätze should work; however, we have not been able to find them.
• Are there weakly turbulent solutions? We have explicitly demonstrated that no turbulent behavior occurs within the three-dimensional invariant subspace (52). This does not in principle exclude turbulence for other initial conditions; however, we find it unlikely because the conformal flow generally appears to display less energy spread than the Szegő flow. • To what extent can the highly structured dynamics of the conformal flow be transferred to the original conformally invariant cubic wave equation (5)
Counting the number of relevant pairings of powers of θ in the two polynomials of the product reproduces the coefficient on the left-hand side of (93), thus verifying the summation formula. Using (93), one can rewrite (92) as (One must have |s| < |v| and |s| < |w| in order for the sum to converge.) Multiplying by z n (with |z| < |v| and |z| < |w|) and summing over n, one obtains after elementary simplification i∂ t ∂ z (zu) = 1 2πi Since by construction, the v-integration contour must not enclose any singularities of u, the integral inside the square brackets is simply given by a sum of two residues, resulting in the complex plane representation we have quoted in the main text:

B A few summation formulas
We assemble below the explicit summation formulas necessary for deriving (56-58) from (54). A crucial feature of these expressions in our context is that all of them are proportional to (n + 1).