Threshold for blowup for the supercritical cubic wave equation

TL;DR

This paper numerically investigates the threshold for blowup in the supercritical cubic wave equation, analyzing stability and spectral properties of self-similar solutions in dimensions five and seven.

Contribution

It provides numerical evidence for the blowup threshold and explores the spectral stability of self-similar solutions in supercritical dimensions.

Findings

Self-similar blowup solution is at the threshold between blowup and dispersion in dimensions 5 and 7.

Spectral analysis supports the stability properties of the blowup solution.

Numerical computations of the spectrum in various supercritical dimensions were performed.

Abstract

We consider the focusing cubic wave equation in the energy supercritical case, i.e., in dimensions $d \geq 5$. For this model an explicit nontrivial self-similar blowup solution was recently found by the first and third author in \cite{GlogicSchoerkhuber}. 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.

Figures

Tables

Table 1. Table 1. All eigenvalues of Eq. ( 2.19 ) appear to be real. In addition to two non-negative eigenvalues, there seem to be infinitely many negative ones, the largest three of them being shown.

$d$	${\lambda }_1$	${\lambda }_0$	${\lambda }_{-1}$	${\lambda }_{-2}$	${\lambda }_{-3}$
5	4.37213	1.00000	-0.53721	-1.88858	-3.17611
6	3.39524	1.00000	-0.54896	-1.96235	-3.32602
7	3.00000	1.00000	-0.55242	-2.00000	-3.41077
8	2.78200	1.00000	-0.55388	-2.02356	-3.46699
9	2.64296	1.00000	-0.55462	-2.03986	-3.50744
⋮	⋮	⋮	⋮	⋮	⋮
$\infty$	2.00000	1.00000	-0.55593	-2.13344	-3.76974

Table 2. Table 2. Results of the fit of theoretical prediction to the sub- and supercritical evolutions ( d = 5 𝑑 5 d=5 and d = 7 𝑑 7 d=7 cases), see ( 3.6 ). The value of λ 1 subscript 𝜆 1 \lambda_{1} is fixed according to the results of Sec. ( 2.1 ). Listed values are results of the fit.

$d=5$	subcritical	supercritical
$c$	5.653595752	5.653595793
$a_1$	$-5.423965718\cdot {10}^{-128}$	$5.002789061\cdot {10}^{-126}$
$a_0$	$2.184193927\cdot {10}^{-32}$	$2.144231964\cdot {10}^{-32}$
$a_{-1}$	0.7585737069	0.7585801575
${\lambda }_{-1}$	$-0.5387977976$	$-0.5387995060$
$d=7$	subcritical	supercritical
$c$	3.998463960	3.998463963
$a_1$	$-2.707332535\cdot {10}^{-98}$	$1.138265563\cdot {10}^{-95}$
$a_0$	$1.107517572\cdot {10}^{-34}$	$1.106153203\cdot {10}^{-34}$
$a_{-1}$	0.3001361956	0.3001362363
${\lambda }_{-1}$	$-0.5535087559$	$-0.5535088248$

Equations

$$({\partial }_t^2-\Delta )u(t,x)=u(t,x{)}^3$$

$$u_{\lambda }(t,x)={\lambda }^{-1}u(t/\lambda ,x/\lambda ),\lambda >0.$$

$$u_T(t,x)=\frac{\sqrt{2}}{T-t},T>0,$$

$$u_T^{\ast }(t,x)=\frac{1}{T-t}{U}^{\ast }\left(\frac{|x|}{T-t}\right),U^{\ast }(\rho )=\frac{2\sqrt{2(d-1)(d-4)}}{d-4+3{\rho }^2}$$

$${\partial }_t^{2}u(t,r)-{\partial }_r^{2}u(t,r)-\frac{d-1}{r}{\partial }_{r}u(t,r)=u(t,r{)}^3$$

$$u(0,\cdot )=u_1^{\ast }(0,\cdot )+f,{\partial }_{t}u(0,\cdot )={\partial }_t{u}_1^{\ast }(0,\cdot )+g.$$

$${\mathcal{C}}_T:=\bigcup _{t\in [0,T)}\{t\}\times \overline{{\mathbb{B}}_{T-t}^d}.$$

$$u(t,r)=u_T^{\ast }(t,r)+(T-t{)}^{-1}\phi (-\log (T-t),\frac{r}{T-t})$$

$$f_1(r)=(1+r^2{)}^{-2},g_1(r)=4(1+r^2{)}^{-3}.$$

$$\Vert (f,g){\Vert }_{H^4\times H^3({\mathbb{B}}_2^7)}\le \frac{\delta }{c}$$

$$u(0,\cdot )=u_1^{\ast }(0,\cdot )+f+\alpha f_1,{\partial }_{t}u(0,\cdot )={\partial }_t{u}_1^{\ast }(0,\cdot )+g+\alpha g_2$$

$$\begin{array}{r}(T-t{)}^{k-s_c}\Vert u(t,\cdot )-u_T^{\ast }(t,\cdot ){\Vert }_{{\dot{H}}^k({\mathbb{B}}_{T-t}^7)}\lesssim (T-t{)}^{\omega }\\ (T-t{)}^{k-s_c}\Vert {\partial }_{t}u(t,\cdot )-{\partial }_t{u}_T^{\ast }(t,\cdot ){\Vert }_{{\dot{H}}^{k-1}({\mathbb{B}}_{T-t}^7)}\lesssim (T-t{)}^{\omega }\end{array}$$

$$\tau =-\log (T-t),\rho =\frac{r}{T-t}.$$

$$\psi (\tau ,\rho ):=e^{-\tau }u(T-e^{-\tau },e^{-\tau }\rho )=U^{\ast }(\rho )+\phi (\tau ,\rho ),$$

$$\left({\partial }_{\tau }^2+3{\partial }_{\tau }+2\rho {\partial }_{\rho }{\partial }_{\tau }-{\Delta }_{\rho }+{\rho }^2{\partial }_{\rho }^2+4\rho {\partial }_{\rho }+2-V(\rho )\right)\phi (\tau ,\rho )=N(\phi (\tau ,\rho ))$$

$$\phi (\tau ,\rho )=e^{\lambda \tau }f(\rho )$$

$$\left({\partial }_{\tau }^2+3{\partial }_{\tau }+2\rho {\partial }_{\rho }{\partial }_{\tau }-{\Delta }_{\rho }+{\rho }^2{\partial }_{\rho }^2+4\rho {\partial }_{\rho }+2-V(\rho )\right)\phi (\tau ,\rho )=0.$$

$${\varphi }_0(t,r):=\frac{d}{d\epsilon }{u}_{T+\epsilon }^{\ast }(t,r){|}_{\epsilon =0}=(T-t{)}^{-2}{f}_0(\frac{r}{T-t})$$

$$f_0(\rho )=\frac{d-4-3{\rho }^2}{{\left(d-4+3{\rho }^2\right)}^2},$$

$$\left({\partial }_t^2-{\partial }_r^2-\frac{d-1}{r}{\partial }_r-3(T-t{)}^{-2}{U}^{\ast }(\frac{r}{T-t}{)}^2\right)\varphi (t,r)=0.$$

$$\begin{gathered} (1-{\rho }^2)f^{\prime \prime }(\rho )+\left(\frac{d-1}{\rho }-2(\lambda +2)\rho \right)f^{\prime }(\rho ) \\ -\left((\lambda +1)(\lambda +2)-\frac{24(d-1)(d-4)}{(d-4+3{\rho }^2{)}^2}\right)f(\rho )=0. \end{gathered}$$

$$x={\rho }^2.$$

$$f(\rho )=\frac{y(x)}{(d-4+3x{)}^2}.$$

$$\begin{gathered} y^{\prime \prime }(x)+\left(\frac{d}{2x}+\frac{2\lambda +5-d}{2(x-1)}-\frac{12}{3x+d-4}\right)y^{\prime }(x) \\ +\frac{3(\lambda -3)(\lambda -2)x+\lambda (\lambda +3)(d-4)-10d+16)}{4x(x-1)(3x+d-4)}y(x)=0. \end{gathered}$$

$$y_0(x)=1+\sum _{n=1}^{\infty }{a}_n{x}^n\text{and}{y}_1(x)=1+\sum _{n=1}^{\infty }{b}_n(1-x{)}^n$$

$$W[y_0,y_1](x):=y_0^{\prime }(x)y_1(x)-y_0(x)y_1^{\prime }(x)$$

$$f_0(\rho )=\frac{1-{\rho }^2}{(1+{\rho }^2{)}^2}\text{and}{f}_1(\rho )=\frac{1}{(1+{\rho }^2{)}^2}$$

$$x=dz$$

$$\begin{array}{r}z{y}^{\prime \prime }(z)+\left(-\frac{1}{2z}+\frac{4}{3z+1}+\lambda -\frac{3}{2}\right)y^{\prime }(z)+\frac{\lambda -2}{4}\frac{3(\lambda -3)z+\lambda +5}{z(3z+1)}y(z)=0.\end{array}$$

$$y(z)=\sum _{n=0}^{\infty }{a}_n(\lambda )z^n.$$

Full text

Irfan Glogić is supported by the Austrian Science Fund FWF, Project P 30076. Birgit Schörkhuber gratefully acknowledges support from the Klaus-Tschira Stiftung. Computations have been performed on Minerva cluster of the Max-Planck Institute for Gravitational Physics. Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID 258734477 - SFB 1173

1. Introduction

We consider the focusing cubic wave equation

[TABLE]

$(t,x)\in I\times\mathbb{R}^{d}$, where $I\subset\mathbb{R}$ is an interval containing zero. The model is invariant under rescaling $u\mapsto u_{\lambda}$

[TABLE]

The corresponding scale invariant Sobolev space for $(u(t,\cdot),\partial_{t}u(t,\cdot))$ is $\dot{H}^{s_{c}}\times\dot{H}^{s_{c}-1}(\mathbb{R}^{d})$, $s_{c}=\frac{d}{2}-1$. Note that $s_{c}=1$ for $d=4$ and in this case, Eq. (1.1) is referred to as energy critical. In this paper, we restrict ourselves to the supercritical case $d\geq 5$. Eq. (1.1) admits a trivial self-similar blowup solution,

[TABLE]

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

[TABLE]

for $d\geq 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 Sec. 2, we address the stability of $u^{*}_{T}$ and numerically investigate the underlying spectral problem for general dimensions. In Sec. 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.

1.1. 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 $\mathbb{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 investigation 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.

1.2. 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 Eq. (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 Sec. 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\geq 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\in\mathbb{N}_{0}\}$ in dimensions $5\leq d<13$, with $U_{0}=\sqrt{2}$. In $d\geq 13$, a result by Collot [13] proves the existence of non-self-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 Eq. (1.1) in $d=4$ have been obtained for example in [14], [15], [16] and [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.

2. Co-dimension one stable self-similar blowup

To investigate the stability of $u_{T}^{*}$ in the spherically symmetric setting, we fix $d\geq 5$ and consider the radial cubic wave equation

[TABLE]

for small radial perturbations $(f,g)$ of the blowup initial data

[TABLE]

Since we are interested in the blowup behavior near the origin, we consider the evolution in backward light cones

[TABLE]

We look for solutions that can be written as

[TABLE]

for some $T>0$, such that the perturbation $\varphi$ vanishes in a suitable sense as $t\to T^{-}$. The result of [1] proves this in $d=7$ under a co-dimension one condition on the initial data.

Theorem 2.1** ([1], radial version).**

Let $d=7$ and define

[TABLE]

There are constants $\omega,\delta,c>0$ such that for all smooth, radial $(f,g)$ with

[TABLE]

the following holds: There are $\alpha\in[-\delta,\delta]$ and $T\in[1-\delta,1+\delta]$ depending Lipschitz continuously on $(f,g)$ such that for initial data

[TABLE]

there is a unique solution $u$ in the backward light cone $\mathcal{C}_{T}$ blowing up at $t=T$ and converging to $u^{*}_{T}$ according to

[TABLE]

for $k=1,2,3$.

We note that the right-hand side of Eq. (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 $\varphi$ in self-similar coordinates

[TABLE]

Note that the backward light cone is mapped to $\bigcup_{\tau\in[0,\infty)}\{\tau\}\times[0,1]$. By defining

[TABLE]

the evolution for $\varphi$ is given by

[TABLE]

with $V(\rho)=3U^{*}(\rho)^{2}$ and $N(\varphi)=(U^{*}+\varphi)^{3}-3{U^{*}}^{2}$. 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) space-dimensions is the solution of the spectral problem for the corresponding linearization. As a matter of fact, it suffices to study mode solutions

[TABLE]

with $\lambda\in\mathbb{C}$, $\mathrm{Re}\lambda\geq 0$ and smooth profiles $f\in C^{\infty}[0,1]$, that satisfy the linearized equation

[TABLE]

In a rigorous formulation, the values $\lambda$ correspond to eigenvalues of a suitably defined differential operator $\mathbf{L}_{0}$, see [1].

2.1. The spectral problem

A fundamental observation is that the time translation invariance of the blowup solution induces a mode solution for $\lambda_{0}=1$. More precisely,

[TABLE]

with

[TABLE]

satisfies the linearized equation in physical variables $(t,r)$,

[TABLE]

Hence, $\varphi_{0}(\tau,\rho)=e^{-\tau}\phi_{0}(T-e^{-\tau},e^{-\tau}\rho)=e^{\tau}f_{0}(\rho)$ 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 Eq. (2.12) we obtain the following ordinary differential equation

[TABLE]

We are interested in values of $\lambda$ that yield solutions that are smooth on $[0,1]$. Note that (2.16) is a Fuchsian equation with six singular points. In particular, $\rho=0$ and $\rho=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 Eq. (2.16). This is done by the change of the independent variable

[TABLE]

In this way, the set of (regular) singular points becomes $\{0,1,(4-d)/3,\infty\}$. 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

[TABLE]

In this way we are lead to the following equation

[TABLE]

Since the set of Frobenius indices of Eq. (2.16) at $\rho=0$ is $\{0,2-d\}$ analytic solutions at $\rho=0$ are even and therefore (2.17) preserves analyticity at the origin. What is more, Eq. (2.16) and Eq. (2.19) have the same set of eigenvalues.

The sets of Frobenius indices for Eq. (2.19) at $x=0$ and $x=1$ are $\{0,1-d/2\}$ and $\{0,(d-3)/2-\lambda\}$ respectively. Therefore, if $(d-3)/2-\lambda$ is not a positive integer then the (normalized) analytic solutions at $x=0$ and $x=1$ are

[TABLE]

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

[TABLE]

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-\lambda$ 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 $\lambda$ 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 Tab. 1.

Interestingly, all eigenvalues appear to be real, even though there is no a priori reason for that. Actually, following the reasoning in [8], Sec. 3, one can prove that the eigenvalues for which $\text{Re}\>\lambda>(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 $\lambda_{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 $\lambda_{1}=3$ with an explicit corresponding eigenfunction, and we were in fact able to prove the following result.

Proposition 2.2** ([1], case $\ell=0$).**

Let $d=7$. For $\mathrm{Re\,}\lambda\geq 0$, the only solutions $f\in C^{\infty}[0,1]$ to Eq. (2.16) are (constant multiples of)

[TABLE]

for which $\lambda=1$ and $\lambda=3$ respectively.

On the other hand, it seems that there are infinitely many negative eigenvalues; in Tab. 1 we listed the largest three of them. Here we should point out that (for $d=7$) only the eigenvalues for which $\text{Re\,}\lambda>-1/2$ correspond to isolated eigenvalues of the operator $\mathbf{L}_{0}$, see Lemma 5.2 in [1]. This is related to the fact that we study Eq. (2.10) in $H^{3}(\mathbb{B}^{7})\times H^{2}(\mathbb{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 Tab. 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

[TABLE]

in Eq. (2.19) and by letting $d$ go to infinity we obtain the following equation

[TABLE]

In this process the interval $[0,1]$ contracts into a single point $z=0$. We therefore look for solutions of Eq. (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 $\lambda$. 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 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 Eq. (2.24) centered at zero

[TABLE]

Here the coefficients are given by the recurrence relation

[TABLE]

where

[TABLE]

and

[TABLE]

and the initial condition is $a_{0}(\lambda)=1$, $a_{1}(\lambda)=A_{-1}(\lambda)=(\lambda^{2}+3\lambda-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 Eq. (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 $\lambda\in\mathbb{C}$ there are two linearly independent solutions to Eq. (2.26) with the following behavior

[TABLE]

as $n\rightarrow\infty$. This, together with the fact that $\lim_{n\rightarrow\infty}A_{n}(\lambda)=\lim_{n\rightarrow}B_{n}(\lambda)=\infty$, implies that for every $\lambda\in\mathbb{C}$ there are $c_{1}(\lambda)$ and $c_{2}(\lambda)$ such that

[TABLE]

for large $n$. Also, note that since $a_{n}^{(1)}$ is the dominant solution the choice of the coefficient $c_{1}(\lambda)$ is unique, i.e. it does not depend on the particular choice of the solution $a^{(1)}_{n}$.

Now from Eqs. (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}(\lambda)$ is a minimal solution to Eq. (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

[TABLE]

written in the standard compact form. Pincherle’s theorem, see [24], p. 402, says that (2.29) is convergent if Eq. (2.26) has no minimal solution with initial condition $a_{0}=0$, $a_{1}=1$. In this case, by Pincherle’s theorem again, a solution to Eq. (2.26) is minimal if and only if the ratio $a_{1}/a_{0}$ equals (2.29). Therefore, the zeros of the following transcendental function

[TABLE]

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 Tab. 1.

The preceding conclusion was made under the assumption that $\lambda$ is not such that the solution to Eq. (2.26) with $a_{0}=0$, $a_{1}=1$ is minimal. As a matter of fact, such values of $\lambda$ 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, Tab. 1 suggests that one eigenvalue of Eq. (2.24) is $\lambda_{1}=2$. In fact, this is obvious since any constant function is a solution in that case. Furthermore, $\lambda_{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\rightarrow\infty$, having in mind (2.18) and (2.23). Then, since both (numerically observed) non-negative eigenvalues of Eq. (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 Eq. (2.19) is close to the one of Eq. (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 Eq. (2.19). We checked that it yields the same result as displayed in Tab. 1.

Finally, we remark that in the case $d=7$ the genuine instability $\lambda_{1}=3$ can be related to the conformal invariance of Eq. (2.15) which is due to the self-similar character of the potential. In fact, it is easy to check that the transformation $\phi\mapsto\tilde{\phi}$,

[TABLE]

leaves Eq. (2.15) invariant. In particular,

[TABLE]

gives rise to a solution of Eq. (2.12) given by

[TABLE]

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}(\rho)=\mathcal{O}(1)$ around $\rho=1$ in general. However, in $d=7$, $f_{0}(1)=0$ and thus

[TABLE]

corresponds to the mode solution for $\lambda_{1}=3$. This effect can be observed also in other supercritical wave equations, see e.g. [8] for the wave maps problem.

3. 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 self-similar 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).

3.1. Numerical approach

Following [8] we introduce new (computational) coordinates $(s,y)$ through

[TABLE]

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 $(\tau,\rho)$ defined in (2.8) when $h(s)\equiv 1$. We also define new dependent variables

[TABLE]

The time evolution of the new variables (3.2) follows from (3.1) and the equation (2.1); explicitly we have

[TABLE]

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\in\mathbb{R}$. Most importantly, $P(s,0)\rightarrow 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,\partial_{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 $\partial_{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 Fig. 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 self-adapting 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.

3.2. Results

For the family of initial data

[TABLE]

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/\sqrt{2}$ as $s\rightarrow\infty$, see Fig. 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_{*}\approx 1.710572581$ in $d=5$ and $a_{*}\approx 2.335609125$ 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

[TABLE]

To see the intermediate attractor at the threshold we have to find the approximate blowup time of the critical evolution corresponding to $a\approx a_{*}$. This procedure is given in [8]. By this, we observe that subcritical data approaches $u_{T}^{*}$ for intermediate times before it disperses, see Fig. 2. Furthermore, the behavior at the origin for different values of $a<a_{*}$ is displayed in Fig. 3.

As an additional evidence for the threshold nature of $u_{T}^{*}$, we pass to self-similar coordinates $(\tau,\rho)$ using the computed critical blowup time, see Fig.  4 and Fig.  5.

Finally, we perform a qualitative comparison of the numerical data with the analytical predictions of Sec. 2.1. We expect a solution in self-similar coordinates, recall Eq. (2.9), to behave like

[TABLE]

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}\sim a-a_{*}\equiv\varepsilon$ is small111Using 128 digits of precision for $d=5$ we were able to fine-tune $a$ up to $\varepsilon\approx 10^{-129}$ whereas for $d=7$ and higher with 96 digits we obtain $\varepsilon\approx 10^{-96}$. so that for long enough time we see the convergence of $\psi(\tau,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 Fig. 6, while the results of the fits are collected in Tab. 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, $\lambda_{-1}$ is in accordance with the numerically computed value in Tab. 1.

Acknowledgement. We would like to thank Piotr Bizoń and Roland Donninger for helpful discussions. Furthermore, we thank the anonymous referee for his/her useful comments and suggestions.