Threshold for blowup for the supercritical cubic wave equation

We consider the focusing cubic wave equation in the energy supercritical case, i.e. in dimensions . For this model an explicit nontrivial self-similar blowup solution was recently found by the first and third author in Glogić and Schörkhuber (2018 (arXiv:1810.07681)). Furthermore, the solution is proven to be co-dimension one stable in d  =  7. In this paper, we study the equation from a numerical point of view. For d  =  5 and d  =  7 in the radial case, we provide evidence that this solution is at the threshold between generic ODE blowup and dispersion. In addition, we investigate the spectral problem that underlies the stability analysis and compute the spectrum in general supercritical dimensions.


Introduction
We consider the focusing cubic wave equation The corresponding scale invariant Sobolev space for (u(t, ·), ∂ t u(t, ·)) is Ḣ sc ×Ḣ sc−1 (R d ), s c = d 2 − 1. Note that s c = 1 for d = 4 and in this case, equation (1.1) is referred to as energy critical. In this paper, we restrict ourselves to the supercritical case d 5. Equation (1.1) admits a trivial self-similar blowup solution, (1.2) which is known to be stable, see [2,3]. In the recent work [1] by the first and the third author a non-trivial self-similar solution was found for d 5. Furthermore, in d = 7 this solution is proven to be co-dimension one stable. More precisely, there exists a co-dimension one Lipschitz manifold of initial data in a small neighbourhood of u * T , whose solutions blow up in finite time and converge asymptotically to u * T (modulo space-time shifts and Lorentz boosts) in the backward light cone of the blowup point.
In section 2, we address the stability of u * T and numerically investigate the underlying spectral problem for general dimensions. In section 3, we provide evidence that this solution is an attractor within a co-dimension one manifold that is a threshold between ODE blowup and dispersion. To finish the introduction, we give a short overview on threshold phenomena in energy supercritical wave equations, which is our main motivation to study the model at hand.

Threshold for blowup in energy supercritical models
In the past decades, a vast body of literature has been concerned with the question of global existence versus finite-time blowup of solutions in nonlinear PDEs. Intuitively, one expects 'small' initial data to lead to global in time solutions, while singularities are supposed to form from sufficiently 'large' data. Consequently, one could expect the existence of a certain threshold between these two basins of attraction. From a mathematical point of view, this is highly non-trivial to be made precise; in particular in energy supercritical models, where little is known in general.
In the context of gravitational collapse, this question has been investigated numerically in the 90s by Choptuik [4] for a simple matter model. By considering one-parameter families of solutions interpolating between dispersion and black hole formation, he found a (discretely) self-similar solution as an intermediate attractor for solutions close to the threshold between these two end states. From a physical point of view, this solution violates conjectures about the structure of gravitational singularities. In the past years, numerous simulations have investigated so-called critical phenomena in more involved models, see e.g. [5] for a review. It was found that threshold solutions between different stable regimes are either static or self-similar (discretely or continuously). Remarkably, the existence of the Choptuik critical solution has recently been proven by Reiterer and Trubowitz in [6] by using computer-assisted methods.
In the recent years it has become evident that self-similar threshold solutions seem to be a feature of energy supercritical models (including Einstein's equation) rather than being specific to gravitational collapse: for the wave maps equation, a critical self-similar blowup solution was observed numerically in d = 3 by Bizoń, Chmaj and Tabor [7] and by Biernat, Bizoń and the second author in higher space dimensions [8]. Other examples are the Yang-Mills equation on R 5+1 [9], and the supercritical focusing wave equation in three dimensions [10]. From an analytic point of view, these problems are widely open. In fact, the model considered in this paper seems to be the first wave equation for which a candidate for a critical self-similar solution is known in closed form. Thus, it provides a good starting point for the analytic invest igation of threshold phenomena in supercritical wave equations.
We remark in passing that critical self-similar solutions have been observed numerically in other problems such as the three-dimensional parabolic-elliptic Keller-Segel model [11]. There, the situation seems to be even more complex due to the existence of different stable blowup regimes and at least two different critical self-similar solutions.

Blowup dynamics for the supercritical cubic wave equation
To complete our discussion, we mention some known results about singularity formation for the model under investigation. For a more thorough discussion on the focusing wave equation, we refer the reader e.g. to [3].
For equation (1.1), the generic blowup behavior is conjectured to be governed by the ODE blowup solution (1.2). In fact, this is also corroborated by our numerical simulations presented in section 3. From a rigorous point of view, the stability of this solution (in a local sense) has been proved by Donninger and the third author in [3] for all odd d 5 in spherical symmetry. The non-radial case was addressed by Donninger and Chatzikaleas [2] in d = 5.
Non-trivial self-similar solutions have been investigated numerically by Kycia [12]; from that one can expect the existence of infinitely many profiles {U n : n ∈ N 0 } in dimensions 5 d < 13, with U 0 = √ 2. In d 13, a result by Collot [13] proves the existence of nonself-similar blowup solutions with more than two unstable directions.
We conclude this section by briefly commenting on the energy critical case. There, the threshold is characterized in terms of the static ground state solution and this is fairly well understood from an analytic point of view. Results for equation (1.1) in d = 4 have been obtained for example in [14][15][16][17]. We also refer the reader to [18] for a characterization of the threshold for an energy critical wave equation in three dimensions. However, these techniques are specific to the critical case and cannot be transferred directly to supercritical problems.

Co-dimension one stable self-similar blowup
To investigate the stability of u * T in the spherically symmetric setting, we fix d 5 and consider the radial cubic wave equation for small radial perturbations ( f , g) of the blowup initial data Since we are interested in the blowup behavior near the origin, we consider the evolution in backward light cones We look for solutions that can be written as (2.4) for some T > 0, such that the perturbation ϕ vanishes in a suitable sense as t → T − . The result of [1] proves this in d = 7 under a co-dimension one condition on the initial data. (2.5) There are constants ω, δ, c > 0 such that for all smooth, radial ( f , g) with δ c the following holds: There are α ∈ [−δ, δ] and T ∈ [1 − δ, 1 + δ] depending Lipschitz continuously on ( f , g) such that for initial data there is a unique solution u in the backward light cone C T blowing up at t = T and converging to u * T according to for k = 1, 2, 3.
We note that the right-hand side of equation (2.7) is normalized to the behavior of u * T in the respective norm. Furthermore, by a solution, we mean a solution to the corresponding operator equation in adapted coordinates, see [1]. The co-dimension one condition is formulated in terms of explicit functions ( f 1 , g 1 ) which arise as eigenfunctions in the spectral problem for the linearization around the blowup solution. This will be explained in more detail below.
The proof of theorem 2.1 relies on the analysis of the time evolution for the perturbation ϕ in self-similar coordinates the evolution for ϕ is given by It can be shown that the dynamics are governed by the linearized problem and that the nonlinearity N on the right-hand side can be treated perturbatively. The key to a generalization of theorem 2.1 to other (odd) spacedimensions is the solution of the spectral problem for the corresponding linearization. As a matter of fact, it suffices to study mode solutions with λ ∈ C, Reλ 0 and smooth profiles f ∈ C ∞ [0, 1], that satisfy the linearized equation In a rigorous formulation, the values λ correspond to eigenvalues of a suitably defined differ ential operator L 0 , see [1].

The spectral problem
A fundamental observation is that the time translation invariance of the blowup solution induces a mode solution for λ 0 = 1. More precisely, satisfies the linearized equation in physical variables (t, r), is a solution of (2.12). In the nonlinear time-evolution this gauge mode can be controlled by suitably adjusting the blowup time T > 0 and is not an obstruction to stability. However, as will be discussed in the following, there is an additional genuine instability that yields the co-dimension one condition on the data. By inserting the mode ansatz (2.11) into equation (2.12) we obtain the following ordinary differential equation We are interested in values of λ that yield solutions that are smooth on [0, 1]. Note that (2.16) is a Fuchsian equation with six singular points. In particular, ρ = 0 and ρ = 1 are singularities and Frobenius theory implies that eigenfunctions are not just smooth but analytic on [0, 1]. This observation allows for an effective use of the shooting method to compute the eigenvalues. This standard technique relies on approximately computing analytic solutions emanating from two singular points and then adjusting the underlying parameter so that solutions smoothly match at some third, conveniently chosen point. Our particular approach follows the work of Biernat and Bizoń [19] on an analogous problem relative to supercritical wave maps.
To simplify the analysis we first reduce the number of singular points of equation (2.16). This is done by the change of the independent variable x = ρ 2 . (2.17) In this way, the set of (regular) singular points becomes {0, 1, (4 − d)/3, ∞}. We remark that Fuchsian equations with four singular points go under the name of Heun, and to bring the equation to its canonical Heun form, see [20], we also scale the dependent variable (2.18) In this way we are lead to the following equation Since the set of Frobenius indices of equation (2.16) respectively. Note that the scenario that yields an eigenvalue is precisely when y 0 and y 1 are constant multiples of each other, which is in turn equivalent to Wronskian being identically zero. Now we approximately compute the solutions y 0 and y 1 by truncating the series (2.20) for a large n and then we evaluate the Wronskian (2.21) at the midpoint x = 1/2. Of course, the choice of this point may depend on other singularities, as their position (in general) determines the radius of convergence of series (2.20). The (approximate) eigenvalues are then given by the zeroes of the Wronskian. We remark that for the Heun equation there is no closed form expression for the Wronskian, unlike for example for the hypergeometric equation where (2.21) is given in terms of Gamma functions. However, functions the y 0 and y 1 are built into Maple, and we can therefore numerically compute the zeros of W[y 0 , y 1 ](1/2) with relatively high precision. Still, recall that the shooting is done with the assumption that (d − 3)/2 − λ is not a positive integer. Therefore, this complementary case has to be treated separately. Namely, in this situation it can happen that both Frobenius solutions at x = 1 are analytic, in which case the underlying λ is obviously an eigenvalue. Otherwise only the subdominant Frobenius solution is analytic and by factoring out its asymptotic behavior at x = 1 we obtain a new equation which is amenable to the shooting method. Our findings are displayed in table 1.
Interestingly, all eigenvalues appear to be real, even though there is no a priori reason for that. Actually, following the reasoning in [8], section 3, one can prove that the eigenvalues for which Re λ > (d − 3)/2 are necessarily real. Note, however, that this is of little significance as already for d = 9 we observe no eigenvalues in this region. Note also that in each dimension, in addition to λ 0 = 1 there is exactly one more non-negative eigenvalue. However, it is very difficult to rigorously prove this observation. Nonetheless, for d = 7 this eigenvalue happens to be λ 1 = 3 with an explicit corresponding eigenfunction, and we were in fact able to prove the following result.
On the other hand, it seems that there are infinitely many negative eigenvalues; in table 1 we listed the largest three of them. Here we should point out that (for d = 7) only the eigenvalues for which Re λ > −1/2 correspond to isolated eigenvalues of the operator L 0 , see lemma 5.2 in [1]. This is related to the fact that we study equation (2.10) in H 3 (B 7 ) × H 2 (B 7 ) and the spectral cut-off −1/2 is dictated by the choice of (the regularity of) this space. What is more, by increasing the Sobolev exponents the spectral cut-off is pushed to the left and in this way more and more negative eigenvalues get uncovered.
A particularly interesting feature of our numerical results in table 1 is that eigenvalues seem to decrease as the dimension increases. It is therefore natural to ask as to whether they have limiting values as d goes to infinity. To further investigate this, we rescale the independent variable x = dz (2.23) in equation (2.19) and by letting d go to infinity we obtain the following equation (2.24) In this process the interval [0, 1] contracts into a single point z = 0. We therefore look for solutions of equation (2.24) that are analytic at z = 0 only. There is indeed a formal power series solution centered at zero, which is however generically not convergent as zero is an irregular singular point. Therefore, by requiring convergence we impose a quantization condition on λ. This is however a different problem from the one above (where the relevant solutions are a priori analytic) and we therefore need a different approach. For this, we use an adaptation of the so-called continued fraction method. This method relies on a remarkable observation Table 1. All eigenvalues of equation (2.19) appear to be real. In addition to two nonnegative eigenvalues, there seem to be infinitely many negative ones, the largest three of them being shown. which was (to the best of knowledge of the authors) first time made by Jaffé in [21]. He makes use of the connection between the three term recurrence relations and continued fractions to compute the bound states of the hydrogen atom. This technique was later popularized inside the general relativity community by Leaver [22], who used it to calculate the quasinormal modes of Kerr black holes. Interestingly, as pointed out by Bizoń [23], although relatively old, this method does not seem to be well-known within the mathematics community. We start with the normalized formal power series solution to equation (2.24) centered at zero y(z) = ∞ n=0 a n (λ)z n . (2.25) Here the coefficients are given by the recurrence relation a n+2 (λ) = A n (λ)a n+1 (λ) + B n (λ)a n (λ) and the initial condition is a 0 (λ) = 1, a 1 (λ) = A −1 (λ) = (λ 2 + 3λ − 10)/2. The radius of convergence of the series (2.25) is determined by the asymptotics of the coefficients a n , and since this sequence satisfies equation (2.26) we make use of the difference equation theory to determine all possible scenarios. For this we refer the reader to an excellent book by Elaydi [24]. First, it can be proved that for every λ ∈ C there are two linearly independent solutions to equation (2.26) with the following behavior a (1) n (λ) ∼ n! 2 n n λ− 1 2 and a (2) n (λ) ∼ (−3) n n −6 (2.27) as n → ∞. This, together with the fact that lim n→∞ A n (λ) = lim n→ B n (λ) = ∞, implies that for every λ ∈ C there are c 1 (λ) and c 2 (λ) such that a n (λ) = c 1 (λ)a (1) n (λ) + c 2 (λ)a (2) n (λ) (2.28) for large n. Also, note that since a (1) n is the dominant solution the choice of the coefficient c 1 (λ) is unique, i.e. it does not depend on the particular choice of the solution a (1) n . Now from equations (2.28) and (2.27) it is clear that the eigenvalues are precisely the zeros of the function c 1 . In the language of difference equations theory this is equivalent to saying that a n (λ) is a minimal solution to equation (2.26). Here minimal solution denotes the one that is asymptotically negligible relative to any other solution that is not its constant multiple (for large values of n). Now, consider the following continued fraction  29). Therefore, the zeros of the following transcendental function are eigenvalues. Furthermore, the roots of this function can be approximately computed by truncating the continued fraction for some large n and then solving the corresponding equation. We displayed our results in the last row of table 1. The preceding conclusion was made under the assumption that λ is not such that the solution to equation (2.26) with a 0 = 0, a 1 = 1 is minimal. As a matter of fact, such values of λ exist and they can be computed by truncating (2.29) for large n and then calculating the zeros of the reciprocal of the resulting function. What is more, it appears that there are infinitely many of these values and in general it is possible that one of them coincides with an eigenvalue causing its 'loss', in a sense that it then fails to be a zero of the function (2.30). However, this does not seem to happen in our case.
Interestingly, table 1 suggests that one eigenvalue of equation (2.24) is λ 1 = 2. In fact, this is obvious since any constant function is a solution in that case. Furthermore, λ 0 = 1 is an eigenvalue with eigenfunction y 0 (z) = 1 − 3z, as can be found by calculating the limit of (2.14) as d → ∞, having in mind (2.18) and (2.23). Then, since both (numerically observed) non-negative eigenvalues of equation (2.24) are known explicitly together with their corresponding eigenfunctions it is likely that an adjustment of spectral techniques developed in [25,26] and [1] would yield a rigorous proof of non-existence of other unstable eigenvalues. Subsequently, by using some kind of perturbative argument to prove that for large values of d the spectrum of equation (2.19) is close to the one of equation (2.24) one would obtain an analog of theorem 2.1 for all large dimensions d.
We note that the continued fraction method can also be used to compute the spectrum of equation (2.19). We checked that it yields the same result as displayed in table 1.
Finally, we remark that in the case d = 7 the genuine instability λ 1 = 3 can be related to the conformal invariance of equation (2.15) which is due to the self-similar character of the potential. In fact, it is easy to check that the transformation φ →φ, leaves equation (2.15) invariant. In particular, gives rise to a solution of equation (2.12) given bỹ For d = 5, this is just an identity. In higher space dimensions, this transformation induces a singularity at the boundary of the cylinder, since f 0 (ρ) = O(1) around ρ = 1 in general. However, in d = 7, f 0 (1) = 0 and thus corresponds to the mode solution for λ 1 = 3. This effect can be observed also in other supercritical wave equations, see e.g. [8] for the wave maps problem.

Threshold behavior
To understand the nature of the threshold for blowup we study the time-evolution from a numerical point of view. The simulations were performed by the second author. We use two different numerical approaches. The first one is very efficient for numerical studies of selfsimilar blowup and follows the scheme introduced in [8]. It uses a suitable redefinition of variables and provides very accurate results for solutions which exhibit self-similar blowup.
Although being very efficient for observing blowup, it is not suitable to simulate solutions that disperse. Therefore we solve the equation using the original (t, r) variables to provide evidence for dispersion of subcritical initial data.
In the following, we present results for d = 5 and d = 7 as an illustration. However, we expect analogous behavior for higher dimensions (we did simulations in d = 9 and d = 11 indicating this).

Numerical approach
Following [8] we introduce new (computational) coordinates (s, y) through t = e −s h(s)ds, r = e −s y. (3.1) The function h introduced in (3.1) is used to make the coordinate transformation adapt to the blowing up solution; note that (s, y) are the self-similar coordinates (τ , ρ) defined in (2.8) when h(s) ≡ 1. We also define new dependent variables V(s, y) = e −s u(t, r), P(s, y) = e −2s ∂ t u(t, r). 3) The main advantage of this rescaling and the redefinition of dependent variables is that V stays finite if we set h(s) = 1/P(s, 0). More precisely, this choice of h leads to the following long time asymptotics at y = 0: V(s, 0) = 1 + Ce −s , for some constant C ∈ R. Most importantly, P(s, 0) → 1/U(0) in the case of self-similar blowup, where U is a self-similar profile. Note that this behavior, finiteness of the dependent variables, is in stark contrast to the original variables (u, ∂ t u) which blow up in finite time. Some care has to be taken to assure that (3.1) defines a coordinate transformation, in particular for our choice of h it is evident that P(s, 0) = 0 for finite s is problematic. It particular, this happens when ∂ t u(t, 0) vanishes for some t. This is the reason why these coordinates are not suitable for the study of solutions that exhibit oscillating behavior. We overcome this by choosing initial conditions which avoid this behavior, as it is apparent from figure 1.
We also note that in order to speed up the numerical calculations we choose the initial data for which the decision on either blowup or dispersion can be made relatively early so that we optimise for computational resources. Other than that the initial data we use are generic. We note that the evolution is stopped when P(s, 0) reaches a small positive value. To understand the behavior beyond this point, we evolve the corresponding initial data in physical coordinates (t, r) and observe dispersion.
As pointed out in [8] the transformation (3.1) with a proper choice of h introduces selfadapting coordinates that accurately resolve both spatial and temporal scales of blowing up solution. In this way we avoid using an adaptive spatial mesh and a time rescaling techniques. In fact to solve (3.3) we use a standard method of lines with 6th order finite difference approximation in space and a 6th order Runge-Kutta method with fixed time step as a time stepping algorithm. Staggered spatial grid with fixed mesh size is used to deal with the y = 0 singularity in (3.3). In addition, a symmetry of the variables (3.2) at the origin is used to construct finite difference stencils close to the coordinate singularity. We add a standard dissipation term to suppress high frequency noise in the data introduced by spatial discretisation. Together with (3.3) we solve the differential of (3.1) with the initial condition t(0) = 0 to get the relation t(s) needed for the subsequent analysis. The code was written in Mathematica whose flexibility and functionality allowed us to use arbitrary precision arithmetics seamlessly. We remark that higher precision is crucial to get close enough to the critical solution and obtain a detailed description of near critical evolutions. To speed up computations we parallelized the bisection search (discussed in the following section) probing the search interval using multiple (typically 64) cores simultaneously.

Results
For the family of initial data we integrate (3.3) forward in s with a > 0 as the only free parameter. For large amplitudes a the solution approaches the ODE blowup (1.2). In particular, P(s, 0) goes to 1/U 0 (0) = 1/ √ 2 as s → ∞, see figure 1. We perform a bisection search in the amplitude a based on the criteria outlined in the previous section; in this way, we find the critical value of a, which we call a * , at threshold. For this value, P(s, 0) is approximately 1/U * (0), which indicates that U * is an attractor within the threshold. Explicitly, a * ≈ 1.710 572 581 in d = 5 and a * ≈ 2.335 609 125 in d = 7.
To analyze the subcritical evolution, we evolve initial data for a < a * in (t, r). We solve (2.1) using the same algorithms as used for solving the system (3.3). We note that at the initial time s = t = 0 the transformation in (3.1)-(3.2) implies that r = y , and V(0, y) = u(0, r), P(0, y) = ∂ t u(0, r). (3.5) To see the intermediate attractor at the threshold we have to find the approximate blowup time of the critical evolution corresponding to a ≈ a * . This procedure is given in [8]. By this, we observe that subcritical data approaches u * T for intermediate times before it disperses, see figure 2. Furthermore, the behavior at the origin for different values of a < a * is displayed in figure 3.  As an additional evidence for the threshold nature of u * T , we pass to self-similar coordinates (τ , ρ) using the computed critical blowup time, see figures 4 and 5.
Finally, we perform a qualitative comparison of the numerical data with the analytical predictions of section 2.1. We expect a solution in self-similar coordinates, recall equation (2.9), to behave like ψ(τ , 0) = c + a 1 e λ1τ + a 0 e τ + a −1 e λ−1τ + · · · , (3.6) close to the critical solution, where the dots stand for even faster decaying modes. The bisection procedure ensures that the coefficient of the unstable mode a 1 ∼ a − a * ≡ ε is small 5 so that for long enough time we see the convergence of ψ(τ , 0) to U * (0). We find that the amplitude of the gauge mode is typically larger than the amplitude of the genuine unstable mode. This is due to propagation of errors on different stages of data analysis. We stress that    the appearance of the gauge mode is an artefact of the data preprocessing and is due to the uncertainty of the blowup time T which is then used in the coordinate transformation.
The comparison of the theoretical prediction (3.6) with the numerical results is presented in figure 6, while the results of the fits are collected in table 2. Note that for c we obtain approximately U * (0). The magnitude of a 0 is small and a 1 is small and of opposite sign in super-and subcritical evolutions. Furthermore, λ −1 is in accordance with the numerically computed value in table 1.