Numerical simulations for the energy-supercritical nonlinear wave equation

TL;DR

This paper presents numerical simulations of the energy-supercritical nonlinear wave equation, showing that the critical Sobolev norm remains bounded over time, supporting recent theoretical scattering results.

Contribution

It provides the first numerical evidence that solutions to the energy-supercritical nonlinear wave equation have bounded Sobolev norms, supporting conjectures about their long-term behavior.

Findings

Critical Sobolev norm remains bounded in simulations

Supports conditional scattering results for nonlinear wave equations

Numerical evidence for long-term stability of solutions

Abstract

We carry out numerical simulations of the defocusing energy-supercritical nonlinear wave equation for a range of spherically-symmetric initial conditions. We demonstrate numerically that the critical Sobolev norm of solutions remains bounded in time. This lends support to conditional scattering results that have been recently established for nonlinear wave equations.

Tables

Table 1. Table 1 . Comparison of numerical solutions for various R max subscript 𝑅 R_{\max} where d = 3 𝑑 3 d=3 , p = 6 𝑝 6 p=6 , and u 0 = 4 ​ exp ⁡ ( − r 2 ) subscript 𝑢 0 4 superscript 𝑟 2 u_{0}=4\exp(-r^{2}) and u 1 = 0 subscript 𝑢 1 0 u_{1}=0 in ( 2.3 ).

Time	$R_{\max }$	$u\left(t,r=0\right)$	$u\left(t,r=19\right)$	$\underset{r\in [0,R_{\max }]}{\max }\left|u(t,r)\right|$
$t=5$	$R_{\max }=20$	0.00417100	0.00000000	0.30992600
	$R_{\max }=30$	0.00417100	0.00000000	0.30992600
	$R_{\max }=50$	0.00417100	0.00000000	0.30992600
$t=10$	$R_{\max }=20$	0.00008800	0.00000000	0.15864800
	$R_{\max }=30$	0.00008800	0.00000000	0.15864800
	$R_{\max }=50$	0.00008800	0.00000000	0.15864800
$t=15$	$R_{\max }=20$	0.00000900	0.00000000	0.10659100
	$R_{\max }=30$	0.00000900	0.00000000	0.10659100
	$R_{\max }=50$	0.00000900	0.00000000	0.10659100

Table 2. Table 2 . Numerical illustration of approximate energy conservation, where Δ max ​ E ​ ( t ) := max t ∈ [ 0 , 15 ] ⁡ | E ​ ( t ) − E ​ ( 0 ) | assign subscript Δ 𝐸 𝑡 subscript 𝑡 0 15 𝐸 𝑡 𝐸 0 \Delta_{\max}E(t):=\max_{t\in[0,15]}|E(t)-E(0)| .

Case	$d=3,p=6$		$d=5,p=2$
	$E\left(0\right)$	${\Delta }_{\max }E\left(t\right)/E\left(0\right)$	$E\left(0\right)$	${\Delta }_{\max }E\left(t\right)/E\left(0\right)$
1	164.184	1.6468e-5	6.02927	2.3386e-7
2	1980.14	2.0506e-5	86.1537	1.7213e-7
3	1996.30	1.9631e-5	128.380	1.7939e-7
4	424.502	3.0777e-5	11.8951	1.9420e-7
5	436.317	2.7890e-5	22.3556	2.9567e-7

Equations

$$i{\partial }_{t}u+\Delta u=\mu |u{|}^{p}u,u:{\mathbb{R}}_t\times {\mathbb{R}}_x^d\to \mathbb{C}$$

$$-{\partial }_t^{2}u+\Delta u=\mu |u{|}^{p}u,u:{\mathbb{R}}_t\times {\mathbb{R}}_x^d\to \mathbb{R}.$$

$$E(u)={\int }_{{\mathbb{R}}^d}\frac{1}{2}|\nabla u{|}^2+\frac{\mu }{p+2}|u{|}^{p+2}dx\text{for NLS}$$

$$E(u,{\partial }_{t}u)={\int }_{{\mathbb{R}}^d}\frac{1}{2}|{\partial }_{t}u{|}^2+\frac{1}{2}|\nabla u{|}^2+\frac{\mu }{p+2}|u{|}^{p+2}dx\text{for NLW.}$$

$$u(t,x)\mapsto \{\begin{array}{ll}{\lambda }^{\frac{2}{p}}u({\lambda }^{2}t,\lambda x)& \text{for NLS}\\ {\lambda }^{\frac{2}{p}}u(\lambda t,\lambda x)& \text{for NLW,}\end{array}$$

$$s_c=\frac{d}{2}-\frac{2}{p},$$

$$\underset{t\to \infty }{\lim }\Vert (u(t),{\partial }_{t}u(t))-(v(t),{\partial }_{t}v(t)){\Vert }_{{\dot{H}}^{s_c}\times {\dot{H}}^{s_c-1}}=0.$$

$$-{\partial }_t^{2}u+{\partial }_r^{2}u+\frac{d-1}{r}{\partial }_{r}u=\mu |u{|}^{p}u,u:{\mathbb{R}}_t\times (0,\infty )\to \mathbb{R},$$

$$(d,p,s_{c})=\big{(}3,6,\tfrac{7}{6}\big{)}\quad\text{and}\quad(d,p,s_{c})=\big{(}5,2,\tfrac{3}{2}\big{)}.$$

$$-{\partial }_t^{2}u+\Delta u-u=|u{|}^{p}u,$$

$$\|u\|_{L_{x}^{\rho}(\mathbb{R}^{d})}=\biggl{(}\int_{\mathbb{R}^{d}}|u(x)|^{\rho}\,dx\biggr{)}^{\frac{1}{\rho}}$$

$$\|u\|_{L_{t}^{q}L_{x}^{\rho}(I\times\mathbb{R}^{d})}=\bigl{\|}\,\|u(t)\|_{L_{x}^{\rho}(\mathbb{R}^{d})}\,\|_{L_{t}^{q}(I)}.$$

$$\hat{f}(\xi )=(2\pi {)}^{-\frac{d}{2}}{\int }_{{\mathbb{R}}^d}{e}^{-ix\xi }f(x)dx,\text{so that}f(x)=(2\pi {)}^{-\frac{d}{2}}{\int }_{{\mathbb{R}}^d}{e}^{ix\xi }\hat{f}(\xi )d\xi .$$

$$\Vert u{\Vert }_{{\dot{H}}_x^s({\mathbb{R}}^d)}=\Vert |\nabla {|}^{s}u{\Vert }_{L_x^2({\mathbb{R}}^d)}=\Vert |\xi {|}^s\hat{u}{\Vert }_{L_{\xi }^2({\mathbb{R}}^d)}.$$

$$P_{N}u={\mathcal{F}}^{-1}{\phi }_N\hat{u}.$$

$$\|u\|_{\dot{B}^{s}_{q,\alpha}(\mathbb{R}^{d})}=\bigl{\|}\|N^{s}P_{N}u\|_{L^{q}(\mathbb{R}^{d})}\bigr{\|}_{\ell_{N}^{\alpha}(2^{\mathbb{Z}})}.$$

$$\Vert (u,{\partial }_{t}u){\Vert }_{{\dot{H}}_x^{s_c}\times {\dot{H}}_x^{s_c-1}}^2:=\Vert u{\Vert }_{{\dot{H}}_x^{s_c}}^2+\Vert {\partial }_{t}u{\Vert }_{{\dot{H}}_x^{s_c-1}}^2,\text{where}{s}_c>1,$$

$$\Vert u{\Vert }_{{\dot{B}}_{2,\infty }^{s_c}}=\underset{N\in 2^{\mathbb{Z}}}{\sup }{N}^{s_c}\Vert P_{N}u{\Vert }_{L_x^2}\text{and}\Vert {\partial }_{t}u{\Vert }_{{\dot{B}}_{2,\infty }^{s_c-1}}=\underset{N\in 2^{\mathbb{Z}}}{\sup }{N}^{s_c-1}\Vert P_N{\partial }_{t}u{\Vert }_{L_x^2}.$$

$$N^s\Vert {\phi }_N\hat{u}{\Vert }_{L_{\xi }^2}\lesssim \Vert u{\Vert }_{{\dot{H}}^s}$$

$$\Vert u{\Vert }_{{\dot{B}}_{2,\infty }^{s_c}}\lesssim \Vert u{\Vert }_{{\dot{H}}^{s_c}}\text{and}\Vert {\partial }_{t}u{\Vert }_{{\dot{B}}_{2,\infty }^{s_c-1}}\lesssim \Vert {\partial }_{t}u{\Vert }_{{\dot{H}}^{s_c-1}}.$$

$$\hat{u}(\xi )$$

$$=|\xi {|}^{-\frac{d-2}{2}}{\int }_0^{\infty }{J}_{\frac{d-2}{2}}(r|\xi |)u(r)r^{\frac{d}{2}}dr,$$

$$\displaystyle\partial_{tt}u=\frac{1}{r^{d-1}}\frac{\partial}{\partial r}\Big{(}r^{d-1}\frac{\partial u}{\partial r}\Big{)}-\mu|u|^{p}u,$$

$${\partial }_{r}u(t,r){|}_{r=0}=0,u(t,R_{\max })=0,$$

$$u(0,r)=u_0(r),{\partial }_{t}u(t,r){\mid }_{t=0}=u_1(r)$$

$$r_j=j\Delta r,\text{mbox}forj=0,1,\dots ,N;t_n=n\Delta t,\text{mbox}forn=0,1,\dots .$$

$$\frac{U_j^{n+1}-2U_j^n+U_j^{n-1}}{\Delta t^2}$$

$$\displaystyle\quad\quad=\frac{(\eta_{j}^{+})^{d-1}\big{(}U_{j+1}^{n}-U_{j}^{n}\big{)}+(\eta_{j}^{-})^{d-1}\big{(}U_{j-1}^{n}-U_{j}^{n}\big{)}}{\Delta r^{2}}-\mu|U_{j}^{n}|^{p}U_{j}^{n},$$

$$\displaystyle\frac{U_{0}^{n+1}-2U_{0}^{n}+U_{0}^{n-1}}{\Delta t^{2}}=\frac{d\big{(}U_{1}^{n}-2U_{0}^{n}+U_{-1}^{n}\big{)}}{\Delta r^{2}}-\mu|U_{0}^{n}|^{p}U_{0}^{n},\qquad n\geq 1,$$

$$\frac{U_1^n-U_{-1}^n}{2\Delta r}=0,U_N^n=0,n\ge 0,$$

$$U_j^0=u_0(r_j),\frac{U_j^1-U_j^0}{\Delta t}=u_1(r_j),\text{mbox}for0\le j\le N-1.$$

$$\Delta t\le \Delta r\sqrt{\frac{2^{d-1}}{1+3^{d-1}}}.$$

Full text

1. Introduction

In recent years, there has been a great deal of progress in understanding the long-time behavior of solutions to nonlinear dispersive partial differential equations. One line of research has focused on the scattering problem for large solutions under optimal regularity assumptions on the initial conditions, particularly in the setting of defocusing nonlinear Schrödinger equations (NLS) and wave equations (NLW). Progress in this direction was precipitated especially by the development of new techniques (e.g. the concentration compactness approach to induction on energy) that were developed in order to establish global well-posedness and scattering (i.e. asymptotically linear behavior) for certain special cases, e.g. the mass- and energy-critical NLS. Outside of these special cases, however, current techniques are often limited to proving conditional results, in which one shows that scattering occurs under the assumption of a priori bounds for a critically-scaling Sobolev norm. In this paper, we will present numerical simulations for the energy-supercritical NLW that lend support to the veracity of these assumed critical bounds. A similar study was carried out in [10] in the setting of the energy-supercritical NLS. A related study also appeared in [63] in the setting of the nonlinear Klein–Gordon equation. See also [22] for a numerical study of the boundary of the forward scattering region for the focusing nonlinear Klein–Gordon equation.

To describe the problem and our results more precisely, we introduce the equations

[TABLE]

and

[TABLE]

In each case, the parameter $\mu$ yields either the defocusing ($\mu>0$) or focusing ($\mu<0$) case, and $p>0$ is the power of the nonlinearity. These are Hamiltonian PDE, with the conserved energy given by

[TABLE]

and

[TABLE]

Both equations also enjoy a scaling symmetry, namely

[TABLE]

which defines a notion of critical regularity for these equations. In particular, if we define

[TABLE]

then one finds that the $\dot{H}^{s_{c}}$-norm of $u|_{t=0}$ for NLS and the $\dot{H}^{s_{c}}\times\dot{H}^{s_{c}-1}$ norm of $(u,\partial_{t}u)|_{t=0}$ for NLW are invariant under the rescaling (1.2).111Here $\dot{H}^{s}$ denotes the homogeneous $L^{2}$-based Sobolev space; see Section 1.1. Generally speaking, these are the optimal spaces for initial data in terms of the well-posedness theory of (NLS) and (NLW); see e.g. [7, 47, 8].

The main topic of this paper is the question of scattering. We say that a forward-global solution $u$ to (NLW) scatters (in $\dot{H}^{s_{c}}\times\dot{H}^{s_{c}-1}$) if there exists a solution $v(t)$ to the linear wave equation such that

[TABLE]

An analogous definition holds for solutions to (NLS).

A special case of (NLS) and (NLW), called the energy-critical case, occurs when the scaling symmetry (1.2) leaves the energy of the solution invariant as well. This corresponds to choosing $p=\tfrac{4}{d-2}$ in dimensions $d\geq 3$, or equivalently $s_{c}=1$. For the case of NLS, there is also the mass-critical case corresponding to $p=\tfrac{4}{d}$, in which case the scaling symmetry preserves the mass (i.e. the $L^{2}$ norm), which is a conserved quantity for NLS (but not for NLW). For these special cases, conservation of energy/mass yields a priori control over the critical Sobolev norm (in the defocusing case, at least). Ultimately, this provides enough control over solutions to establish global well-posedness and scattering, although proving this is a very challenging problem that required the work of many mathematicians over many years to settle definitively (see [3, 9, 12, 13, 14, 15, 16, 29, 33, 38, 39, 43, 45, 57, 66, 67, 68, 69, 70, 1, 27, 28, 31, 32, 46, 54, 58, 64]):

Theorem 1.1** (Global well-posedness and scattering).**

For the defocusing case of the mass- and energy-critical NLS or energy-critical NLW, arbitrary initial data in the critical Sobolev space lead to global solutions that scatter. Similar results hold in the focusing case, provided one imposes suitable size restrictions on the mass/energy.

The resolution of Theorem 1.1 required the development of a powerful new set of techniques. The initial breakthrough was due to Bourgain, who introduced the method of ‘induction on energy’ [3]. This technique has been significantly developed and refined. Presently, the typical approach to problems as in Theorem 1.1 follows the so-called ‘Kenig–Merle roadmap’ developed in [33]. One proceeds by contradiction: Assuming the theorem to be false, one constructs a minimal energy counterexample, which (due to minimality) enjoys certain compactness properties. One then shows that such compactness properties are at odds with the dispersive/conservative nature of the equation and ultimately lead to a contradiction; this is often achieved through the use of conservation laws together with certain nonlinear estimates known as virial or Morawetz estimates. For an expository introduction to these techniques, we refer the reader to [44, 71].

Beginning with the work of Kenig and Merle [34], a great deal of recent research has focused on establishing analogous results beyond the mass- and energy-critical cases. In such cases, the ‘Kenig–Merle roadmap’ naturally leads to a proof of scattering under the assumption of a priori bounds in the critical Sobolev space, where the assumed bounds play the role of the ‘missing conservation law’ at critical regularity. Stated roughly, we have the following conjecture:

Conjecture 1.2**.**

For the defocusing NLS or NLW, any solution that remains bounded in the critical Sobolev space is global-in-time and scatters.

By now, the range of positive results of this type is extensive. For the case of NLS, see [34, 40, 36, 37, 51, 52, 53, 72, 20, 49, 48, 26, 73]; for the case of NLW, see [19, 4, 5, 6, 23, 24, 25, 21, 35, 41, 42, 56, 59, 60]. While some recent remarkable work of Dodson [17, 18] has actually established unconditional scattering results at critical regularity for the energy-subcritical NLW with radial initial data, the majority of the scattering results for large data at ‘non-conserved’ critical regularity are conditional in nature. We would also like to mention the recent work of D’Ancona [11], who has established some well-posedness results for the defocusing energy-supercritical NLW in the exterior of a ball.

In this paper, we carry out numerical simulations for the energy-supercritical NLW with radial (i.e. spherically symmetric) initial conditions222Note that radiality is preserved in time. This is a consequence of the fact that the Laplacian commutes with rotations, together with the uniqueness of solutions., where energy-supercritical refers to the condition $p>\tfrac{4}{d-2}$, or equivalently $s_{c}>1$. In the radial setting, (NLW) takes the form

[TABLE]

where we write $u=u(t,r)$, with $r>0$ and impose the Neumann boundary condition $\partial_{r}u|_{r=0}\equiv 0$. As in [10, 63], the restriction to radial solutions provides a significant simplification in the numerical analysis of (NLW). Our main result is to demonstrate (numerically) boundedness of the critical Sobolev norms for large time, thus lending support to the conditional scattering results discussed above. We expect that similar results will hold in the non-radial setting and plan to address this case in future work.

For the sake of concreteness, we focus on two representative cases, namely,

[TABLE]

These particular cases were considered in the works [35, 41, 42, 5, 6], which established scattering under the assumption of a priori bounds for $(u,\partial_{t}u)$ in $\dot{H}^{s_{c}}\times\dot{H}^{s_{c}-1}$. We study a range of choices for $u_{0}=u(0,r)$ and $u_{1}=\partial_{t}u(0,r)$ (see Section 3), and in all cases we observe (numerically) that the critical Sobolev norm converges after a short time and, in particular, remains bounded for large times. Additionally, we compute numerically the potential energy (i.e. the $L^{p+2}$ norm), the $L^{\infty}$ norm, and certain scale-invariant Besov norms. We observe that the Besov norms become relatively small (compared to the Sobolev norms), and that the higher Lebesgue norms decay at a rate that matches solutions to the linear wave equation. All of this behavior is consistent with scattering (see e.g. Section 4). We describe our results in detail in Section 3.

We would also like to compare our results with the work of Strauss and Vazquez [63] on the closely-related defocusing nonlinear Klein–Gordon equation (NLKG)

[TABLE]

which they studied in dimension $d=3$ with radial initial conditions. They considered several choices of nonlinearity, including $p\in\{2,4,6,8\}$ as well as the nonlinearity $\sinh(5u)-5u$. They also observed numerically that solutions should decay in $L^{\infty}$ at a rate matching solutions to the underlying linear problem. Our results may be viewed in part as an extension of those in [63] (albeit in the setting of NLW, rather than NLKG). Indeed, we also observe the sharp $L^{\infty}$ decay in the NLW setting; additionally, we have simulated the solutions over a longer time interval and have computed several other quantities that are related to the problem of scattering (e.g. the Sobolev and Besov norms). We have also treated the case $d=5$, in addition to the three-dimensional case.

At present, existing analytic techniques are generally insufficient to rigorously establish the boundedness in time of the critical Sobolev norm of solutions, unless such norms can be controlled by conserved quantities (although we should mention again the remarkable work of [17, 18] for the case of the radial energy-subcritical NLW). In particular, this problem seems to be especially challenging in the energy-supercritical regime, as there is no known coercive conserved quantity above the regularity of the energy. In general, such boundedness has generally been expected to hold true in the defocusing setting, as both the dispersion of the underlying linear equation and the defocusing nature of the nonlinearity tend to cause solutions to spread out and decay. Our results provide additional numerical evidence in support of the belief of boundedness and lend support to the wide range of conditional scattering results for NLW that have been established in recent years.

It remains an important open problem in the analysis of nonlinear dispersive PDE to determine whether energy-supercritical Sobolev norms do indeed remain bounded in time in general. In the negative direction, we would like to mention the recent preprint of Merle, Raphaël, Rodnianski, and Szeftel [50], who have constructed radial solutions to the defocusing energy-supercritical NLS with energy-supercritical Sobolev norms blowing up in finite time at a polynomial rate! Their result relies on the hydrodynamical formulation of NLS, making use of suitable underlying dynamics for the compressible Euler equations in order to produce a highly oscillatory blowup profile.

The rest of this paper is organized as follows: In Section 1.1, we collect some basic notation and preliminaries. In Section 2, we describe the numerical methods we use in this work. In Section 3, we describe the sets of initial conditions used in the numerical simulations and discuss our numerical findings. In Section 4, we prove a simple scattering result (namely, scattering holds if the critical Besov norm is sufficiently small compared to the critical Sobolev norm), which is relevant to the discussion in Section 3. Finally, in Appendix A we discuss the notion of incoming/outgoing waves, which play a role in our choice of initial conditions.

1.1. Notation and preliminaries

We use the standard notation for Lebesgue norms, e.g.

[TABLE]

for $1\leq\rho<\infty$. We denote space-time norms by $L_{t}^{q}L_{x}^{\rho}$, i.e.

[TABLE]

We define Sobolev and Besov norms by utilizing the Fourier transform, denoted by $\mathcal{F}f$ or $\hat{f}$. We define

[TABLE]

The homogeneous $L^{2}$-based Sobolev spaces are then defined by

[TABLE]

Besov spaces are defined using the standard Littlewood–Paley multipliers. In particular, for $N\in 2^{\mathbb{Z}}$ we let $\varphi_{N}$ denote a smooth bump function supported where $|\xi|\sim N$, with $\sum\varphi_{N}\equiv 1$. We then define the Littlewood–Paley projections $P_{N}u$ through the Fourier transform, i.e.

[TABLE]

The Besov norm $\dot{B}^{s}_{q,\alpha}$ is defined via

[TABLE]

We will frequently consider the Sobolev norm

[TABLE]

as well as the Besov norms

[TABLE]

Note that for any $s$, we have

[TABLE]

uniformly in $N$, which implies

[TABLE]

As mentioned above, the restriction to radial (i.e. spherically symmetric) solutions leads to some simplifications. We have already mentioned the simplification of the PDE (and hence the numerical analysis). Additionally, we may change to spherical coordinates and write

[TABLE]

where $J_{\nu}$ denotes the standard Bessel function (see e.g. [62]). This informs our numerical computation of the Fourier transform, and therefore the relevant Sobolev and Besov norms.

2. Numerical methods

We first truncate (1.4) to a computational domain $[0,R_{\rm max}]$ with $R_{\rm max}$ sufficiently large that the truncating effect can be neglected. We then reformulate (1.4) and solve the following problem:

[TABLE]

Here the homogeneous Dirichlet boundary conditions are considered at $r=R_{\max}$.

In [63], an implicit finite difference method is introduced to solve the NLKG (1.5). This method exactly conserves the discrete energy, but at each time step one has to solve a nonlinear system with iteration methods (e.g. Newton’s method). Hence, its computational costs are high. To avoid solving nonlinear systems at each time step, explicit methods are popular in practice. In [22], an explicit method with a second-order difference scheme for both temporal and spatial discretization is used to study the boundary of the forward scattering region for NLKG. In [10], the authors study the question of boundedness of critical Sobolev norms for NLS using a finite difference scheme in space and the explicit fourth-order Runge–Kutta method in time. Compared to implicit methods, these explicit methods are easier to implement and have less computational costs. They do not have the exact energy conservation, but could provide a good approximation to it if proper numerical parameters are used in simulations [22, 10].

Next, we introduce a second-order finite difference method to solve (2.1)–(2.3). Our method is different from that in [22] where the change of variable $v=ru$ is introduced to move the singularity term from linear to nonlinear part. Instead, we will directly discretize the equation (2.1). Denote the mesh size $\Delta r=R_{\max}/N$ with $N$ a positive integer, and let $\Delta t>0$ be the time step. Then we define the spatial grid points and time sequence as

[TABLE]

We denote the numerical approximation of $u(t_{n},r_{j})$ by $U_{j}^{n}$. Then the equation (2.1) can be approximated by the difference scheme:

[TABLE]

for $1\leq j\leq N-1$ and $n\geq 1$, where we denote $\eta_{j}^{\pm}=(r_{j}\pm{h}/{2})/r_{j}$. For $j=0$, the approximation is given by:

[TABLE]

where $U_{-1}$ denotes the solution at the ghost point $r=r_{-1}$.

The discretization of the boundary conditions (2.2) gives

[TABLE]

which implies that $U_{-1}^{n}\equiv U_{1}^{n}$ for $n\geq 0$. At $t=0$, the initial condition (2.3) is discretized as

[TABLE]

Combining (2.4)–(2.7) yields an explicit second-order finite difference scheme to the problem (2.1)–(2.2). In the simulations, one sufficient stability condition of this scheme is

[TABLE]

This suggests that the simulations of (1.4) for higher dimensions are generally more time-consuming, as a smaller time step is required to ensure the numerical stability.

The discrete energy is calculated by

[TABLE]

for $n\geq 1$, while the Besov norms are calculated following the discretization method in [10].

2.1. Verification of numerical method

We will test the performance of our method for different numerical parameters (i.e. $R_{\rm max}$, $\Delta r$, and $\Delta t$).

First, we study the truncating effect of $R_{\max}$ by fixing the time step $\Delta t$ and mesh size $\Delta r$. In general, a large computational domain $[0,R_{\max}]$ is required to avoid the artifacts from truncation. However, the larger the computational domain, the larger the number of unknowns, and thus the larger the computational costs.

In Table 1, we compare the numerical solutions $u$ that are computed with different choices of $R_{\max}$. For the different $R_{\max}$ chosen here, our method yields the same solutions of $u(t,r=0)$ and $\max_{r\in[0,R_{\max}]}|u(t,r)|$. Moreover, the solution $u(t,r=19)$ remains zero even through time $t=15$, suggesting that $R_{\max}=20$ is large enough for the computation time $t\in[0,15]$.

Next, we will fix $R_{\rm max}=20$ and study the numerical errors for different $\Delta t$ and $\Delta r$. Since the exact solution of (1.4) is unknown, we will use the numerical solution with a very fine mesh (here $\Delta r=2^{-12}$ and $\Delta t=2^{-14}$) as a proxy for the “exact” solution when computing the numerical errors.

Figure 1 shows the $\ell^{2}$-norm errors in solution $u(2,r)$, where an order line of $O(\Delta r^{2})$ is included for easy comparison. We choose the time step and mesh size to satisfy $\Delta t=\Delta r/4$. The results in Figure 1 confirm the second-order accuracy of our methods. It additionally shows that for fixed numerical parameters, the numerical errors for $d=5$ are smaller than those of $d=3$. For example, the numerical error of $\Delta r=0.0098$ is $2.54\times 10^{-4}$ for $d=3$ and $9.38\times 10^{-7}$ for $d=5$.

3. Numerical results

We now discuss our numerical results. To begin, let us discuss the specific choices of initial conditions used in this paper. The key assumption on the data is that of spherical symmetry, which reduces the equation to a one-dimensional problem and thereby greatly simplifies the numerical analysis. We begin by considering the cases of Gaussians that are large enough to be safely outside of the small data regime (for otherwise scattering is a known consequence of the small-data well-posedness theory). As in [10], we would also like to consider some initial data for which the underlying linear equation would experience some initial focusing toward the origin; for such data, we would then like to observe that the defocusing nonlinearity counters this effect. The authors of [10] achieved this by multiplying the Gaussian initial data by $e^{\alpha ir^{2}}$ for some $\alpha>0$. In the setting of the radial wave equation, it seems natural to consider ‘incoming’ initial data for this purpose (see e.g. [2]), which in our setting refers to the condition

[TABLE]

where $u_{0}=u|_{t=0}$ and $u_{1}=\partial_{t}u|_{t=0}.$ We discuss the origin of this condition in more detail in Appendix A.

Apart from choosing incoming/outgoing initial data (which requires defining $u_{1}$ precisely in terms of $u_{0}$), we found that varying the choice of initial velocity $u_{1}$ plays almost no role in terms of the long-time behavior of the solution. Thus, other than the cases for which we take ‘incoming’ data, we will be content to work with the simplest choice $u_{1}=0$. Similar to [10], we will then choose Gaussian data ($u_{0}=C\exp(-r^{2})$), ‘ring’ data ($u_{0}=Cr^{2}\exp(-r^{2})$), or an oscillating Gaussian ($u_{0}=C\exp(-r^{2})\sin(ar)$).

As far as the incoming condition, we note that (3.1) includes the singular term $1/r$ in the expression for $u_{1}$. While one can verify that $u_{1}$ still belongs to $L^{2}\cap\dot{H}^{s_{c}-1}$ (see Lemma A.1), we found that unless $u_{0}$ vanishes to high enough order at $r=0$, it is difficult to simulate this condition numerically. Thus we were led only to consider the incoming condition for the ‘ring’ initial data $u_{0}=Cr^{2}\exp(-r^{2})$ and for the oscillatory data $u_{0}=C\exp(-r^{2})\sin(ar)$, for which cases the presence of $1/r$ is harmless.

Altogether we considered five choices of initial conditions, which are detailed in the subsections below. As described in the introduction, for each case we consider the combinations $(d,p)=(3,6)$ and $(d,p)=(5,2).$ In the following, we summarize the quantities studied, as well as the corresponding findings. For each case, we study the time evolution of:

• (i)

The solution $u(t,r)$. We find that the solution decays over time and travels outward at a constant speed.

• (ii)

The critical Sobolev norms $\|u(t)\|_{\dot{H}^{s_{c}}}$ and $\|\partial_{t}u(t)\|_{\dot{H}^{s_{c}-1}}$, with $s_{c}$ as in (1.3). We find that the critical Sobolev norms may initially oscillate, but quickly settle down and converge to a constant. Moreover, our numerical results show that

[TABLE]

• (iii)

The potential energy $\|u(t)\|_{L^{p+2}}$ and the supremum norm $\|u(t)\|_{L^{\infty}}$. We find that both quantities tend to zero as $t\to\infty$. More precisely, we observe

[TABLE]

for sufficiently large $t$, matching the decay rates for the underlying linear wave equation.

• (iv)

The critical Besov norms $\|u(t)\|_{\dot{B}_{2,\infty}^{s_{c}}}$ and $\|\partial_{t}u(t)\|_{\dot{B}_{2,\infty}^{s_{c}-1}}$. We find that the Besov norms remain bounded and, in fact, become relatively small compared to the critical Sobolev norm as $t\to\infty$. As discussed in Section 4, if solutions are sufficiently dispersed in frequency relative to their Sobolev norm, then one can prove scattering by a small-data type argument.

In what follows, we present our numerical results for our five representative cases, where we will always choose $R_{\max}=20$, $\Delta r=4e$-4, and $\Delta t=1.25e$-4. As discussed previously, our explicit finite difference method does not exactly conserve the energy. However, Table 2 shows that it has a good approximate conservation of energy for a long time (e.g. through time $t=15$).

3.1. Case 1. Gaussian data.

In this case, we take the initial condition

[TABLE]

Figure 2 presents the time evolution of the solution, which shows that the solution decays over time and behaves essentially as a wave packet traveling outward with constant speed.

Next, we further study the decay of the solution $u$. We would like to show decay of both the potential energy (i.e. the $L^{p+2}$-norm), as well as pointwise decay (i.e. the $L^{\infty}$-norm), both of which are consistent with scattering. In fact, we can show that the decay rate for these quantities matches the decay rate for solutions to the linear equation. To see this, we present in Figure 3 the time evolution of the following quantities and observe boundedness (in fact, convergence) for large $t$:

[TABLE]

Figure 4 represents the main result of this paper, namely, numerical evidence for the boundedness of the critical Sobolev norms. In fact, after some initial oscillation, we see that both the $\dot{H}^{s_{c}}$-norm of $u$ and the $\dot{H}^{s_{c}-1}$ norm of $\partial_{t}u$ quickly settle down and converge. Moreover, our numerical results show that

[TABLE]

Finally, we plot the critical Besov norms (namely, $\dot{B}^{s_{c}}_{2,\infty}$ for $u$ and $\dot{B}^{s_{c}-1}_{2,\infty}$ for $\partial_{t}u$) in Figure 5. While these norms do not converge to zero, we can observe that they become relatively small compared to the critical Sobolev norms in Figure 4. As discussed in Section 4, it is possible to prove scattering in such a scenario.

3.2. Case 2. Ring data.

In Case 2 we take

[TABLE]

As we observed the same phenomena as in Case 1, we will be somewhat brief in our presentation.

Figure 6 presents the time evolution of the solution. It shows that the solution is more oscillating than the case with Gaussian initial conditions; however, similar to Case 1, the solution still decays over time and travels outward.

Figure 7 presents the critical Sobolev norms, Besov norms, and higher Lebesgue norms over time. It shows that the critical Sobolev norms $\|u\|_{\dot{H}^{s_{c}}}$ and $\|\partial_{t}u\|_{\dot{H}^{s_{c}-1}}$ quickly converge and stay at the same constant after a relatively small time, while the Besov norms are eventually bounded by a relatively small constant compared to the Sobolev norms. The time evolution of the higher Lebesgue norms shows that the decay of $\|u\|_{L^{p+2}}$ and $\|u\|_{L^{\infty}}$ are qualitatively the same as in Case 1.

3.3. Case 3. Incoming ring data.

In Case 3 we take

[TABLE]

As described above, such initial data will lead to some initial ‘focusing’ at the level of the linear wave equation. This effect is countered by the defocusing nonlinearity.

Figure 8 illustrates the time evolution of the solution $u$, where the plot for $t=0$ is the same as that in Figure 6. Due to the different initial velocity $u_{1}$, the evolution of solution in Cases 2 and 3 is initially quite different (cf. Figures 6 and 8 for $t=1$), but after a long time the influence of initial velocity greatly decreases.

The time evolution of the norms in Figure 9 is quite similar to the observations in Figure 7. This suggests that the initial velocity $u_{1}$ ultimately plays an insignificant role in the study of long-time behavior of solution.

3.4. Case 4. Oscillatory Gaussian data

Next, we consider the initial data as follows:

[TABLE]

The numerical results are consistent with all previous cases: the pointwise decay and outward travel of the solution, with higher norms decaying at rates consistent with the linear wave equation. Moreover, the Sobolev norms settle down, and the Besov norms become small compared to the Sobolev norms.

3.5. Case 5. Incoming oscillatory Gaussian data

Finally, we consider the same initial position as $u_{0}$, but with incoming data:

[TABLE]

Once again, we find that making the data incoming ultimately has little effect on the long-time behavior of solutions, and the numerical results are consistent with all that we have seen above.

4. A simple scattering result

In this section, we demonstrate that if a solution to (NLW) is sufficiently dispersed in frequency compared to its critical Sobolev norm, then it scatters. We present this result because the numerical results presented above demonstrate that this phenomenon occurs, at least for the sets of initial conditions that we considered. We note, however, that this result is still fundamentally a small-data result. It would require some significant new input to rigorously demonstrate that such dispersion in frequency occurs in general.

For the sake of concreteness, we focus on the particular case $(d,p)=(3,6)$ considered above, in which case $s_{c}=\frac{7}{6}$.

Proposition 4.1**.**

Suppose $\|(u_{0},u_{1})\|_{\dot{H}^{\frac{7}{6}}\times\dot{H}^{\frac{1}{6}}}=E$. Then there exists $\eta_{0}=\eta_{0}(E)>0$ such that if

[TABLE]

then the solution to $-\partial_{t}^{2}u+\Delta u=|u|^{6}u$ with data $(u_{0},u_{1})$ is global in time and obeys the space-time bounds

[TABLE]

In particular, $u$ scatters.

To prove this result, it is useful to complexify the equation and introduce the variable

[TABLE]

Then $(u,\partial_{t}u)\in\dot{H}^{s_{c}}\times\dot{H}^{s_{c}-1}$ if and only if $v\in\dot{H}^{s_{c}}$ (with comparable norms), with a similar statement concerning the Besov norms. The equation (NLW) is then equivalent to

[TABLE]

which in turn is equivalent to the Duhamel formulation

[TABLE]

The proof of Proposition 4.1 will rely primarily on Strichartz estimates for the operator $e^{it|\nabla|}$. To begin, we have the following estimates (see, for example, [71]): Let $(q,r)$ and $(\tilde{q},\tilde{r})$ be wave admissible in $3d$, i.e.

[TABLE]

Defining

[TABLE]

we have:

[TABLE]

The specific space-time norms appearing Proposition 4.1 are critical norms for the case $(d,p)=(3,6)$. That is, they are invariant under the rescaling (1.2).

We define

[TABLE]

By the Strichartz estimates above (noting that $\gamma(4,4)=\tfrac{1}{2}$), we have

[TABLE]

A less standard ingredient for the proof of Proposition 4.1 will be the following refined Strichartz estimate for this norm.

Lemma 4.2** (Refined Strichartz estimate).**

There exists $\theta\in(0,1)$ such that

[TABLE]

Proof.

As the Besov norm is controlled by the Sobolev norm, it suffices to prove an estimate of this form for each norm appearing in (4.5). Let us focus on the $L_{t,x}^{4}$ norm, as the proof for the $L_{t,x}^{12}$ norm follows along similar lines.

We employ the Littlewood–Paley frequency decomposition

[TABLE]

as described in Section 1.1. Let us denote

[TABLE]

By the Littlewood–Paley square function estimate, we may write

[TABLE]

Now let $\varepsilon>0$ be a small parameter and define two pairs of sharp wave-admissible exponents $(q_{l},r_{l})$ and $(q_{h},r_{h})$ by

[TABLE]

Note that

[TABLE]

In particular, by (4.3) and Bernstein inequalities333Bernstein inequalities refer to the following general estimates for frequency-localized functions:

$\displaystyle\|P_{N}f\|_{L^{r_{2}}(\mathbb{R}^{d})}$ $\displaystyle\lesssim N^{\frac{d}{r_{2}}-\frac{d}{r_{1}}}\|P_{N}f\|_{L^{r_{1}}(\mathbb{R}^{d})},\quad 1\leq r_{1}\leq r_{2}\leq\infty,$

$\displaystyle\||\nabla|^{s}P_{N}f\|_{L^{r}(\mathbb{R}^{d})}$ $\displaystyle\sim N^{s}\|P_{N}f\|_{L^{r}(\mathbb{R}^{d})},\quad 1<r<\infty,\quad s\in\mathbb{R}.$

, we have

[TABLE]

Thus, continuing from above and applying (4.3) and Cauchy–Schwarz, we have

[TABLE]

which yields the desired estimate.

In the case of the $L_{t,x}^{12}$ norm, the situation is actually simpler because we work away from the ‘strictly admissible’ line $\tfrac{1}{q}+\tfrac{1}{r}=\tfrac{1}{2}$. In this case we put one of the lowest frequency pieces $u_{N_{1}}$ in $L_{t}^{12}L_{x}^{12+}$ and one of the highest frequency pieces $u_{N_{6}}$ in $L_{t}^{12}L_{x}^{12-}$. All of the remaining pieces are taken out with a supremum in $N$ of the $L_{t,x}^{12}$ norm, which yields the Besov norm after an application of Strichartz (as above). Using Bernstein and Strichartz for the low and high frequency pieces, we get a gain of $(\tfrac{N_{1}}{N_{6}})^{\varepsilon}$, which can be used to defeat the logarithmic loss from the sums in $N_{2},\dots,N_{5}$ and to sum in $N_{1},N_{6}$ via Cauchy–Schwarz, just as above. ∎

We turn to the proof of Proposition 4.1.

Proof of Proposition 4.1.

The standard well-posedness theory for (NLW) yields a local-in-time solution to (4.1) satisfying the Duhamel formula (4.2). It will therefore suffice to prove the estimate

[TABLE]

on any interval $I\subset\mathbb{R}$ containing $t=0$, where $\theta\in(0,1)$ is as in Lemma 4.2. Indeed, choosing $\eta_{0}=\eta_{0}(E)$ sufficiently small (where $E,\eta_{0}$ are as in the statement of Proposition 4.1), we may make the first term on the right-hand side of (4.6) as small as we wish. Then, by a standard continuity argument, (4.6) implies

[TABLE]

giving the desired global space-time bounds for $v$.

It therefore remains to prove (4.6). Recalling the Duhamel formula (4.2), we apply Lemma 4.2, (4.4), and the fractional chain rule444This refers to the estimate

$\||\nabla|^{s}F(u)\|_{L^{r}}\lesssim\|F^{\prime}(u)\|_{L^{r_{1}}}\||\nabla|^{s}u\|_{L^{r_{2}}}$

for $0<s<1$ and $1<r,r_{1},r_{2}<\infty$ satisfying $\tfrac{1}{r}=\tfrac{1}{r_{1}}+\tfrac{1}{r_{2}}$. to estimate

[TABLE]

giving (4.6), as desired. ∎

Appendix A Incoming waves for the radial wave equation

In this section we describe the notion of incoming/outgoing waves for the radial wave equation. In particular, we wish to explain the origin of the ‘incoming’ condition

[TABLE]

Our discussion will be brief—for more details, see [2].

First, consider the radial wave equation in three space dimensions:

[TABLE]

If $u$ and $v$ are related through

[TABLE]

then one finds that $u$ solving (A.2) is equivalent to $v$ solving the $1d$ wave equation

[TABLE]

however, we should now view $v$ as a solution on the half-line $r\in(0,\infty)$ with the Neumann boundary condition $\partial_{r}v|_{r=0}=0$.

If instead $u$ solves the radial equation in five space dimensions, i.e.

[TABLE]

then (after a short computation) one finds that the same situation arises when $u$ and $v$ are related via

[TABLE]

Now, for the $1d$ wave equation (with Neumann boundary conditions) one can use the explicit formula for solutions to see that every solution may be decomposed into an incoming and outgoing piece, i.e.

[TABLE]

where $v_{in}$ moves toward $r=0$ and $v_{out}$ moves toward $r=\infty$ (both with speed one). As $t\to\infty$, the solution becomes increasingly outgoing, while as $t\to-\infty$ the solution becomes increasingly incoming. Writing $(v_{0},v_{1})$ for the initial conditions of $v(t)$, one finds that at time zero the incoming component is

[TABLE]

while the outgoing component is

[TABLE]

where $\partial_{r}^{-1}$ denotes the antiderivative. For more detail, see [2].

We are interested in prescribing ‘incoming’ initial data. In this case, at the linear level there is an initial ‘focusing’ of the solution toward the origin, which we can then observe is countered by the defocusing effect of the nonlinearity. In particular, an initial data pair for $v$ is incoming if $v_{+}=0$, that is, if

[TABLE]

Let us now derive equivalent conditions at the level of $(u_{0},u_{1})$. These are the conditions appearing in our choices of initial data.

For the three-dimensional case, we have the simple relation $v=\partial_{r}(ru)$, and the incoming condition becomes

[TABLE]

For the five-dimensional case, we instead get $\partial_{r}^{3}(r^{3}u)=2v+\partial_{r}v$, and the incoming condition becomes

[TABLE]

In general, one derives the condition (A.1).

We close this section with the following lemma, which shows that for sufficiently regular $u_{0}$, the initial velocity velocity $u_{1}$ still belongs to $L^{2}\cap\dot{H}^{s}$ for suitable $s$ (despite the presence of the singular term $1/r$).

Lemma A.1**.**

Let $u$ be a Schwartz function in $\mathbb{R}^{d}$, with $d\geq 3$. Then $\tfrac{1}{|x|}u$ belongs to $\dot{H}^{s}$ for any $0\leq s\leq\min\{\tfrac{d}{2}-1,1\}$. In fact,

[TABLE]

for any such $s$.

Proof.

It is enough to prove (A.6). This estimate is very similar to Hardy’s inequality, which states

[TABLE]

(see e.g. [65]). In particular, an application of Hardy’s inequality (and the boundedness of Riesz potentials) reduces (A.6) to the commutator estimate

[TABLE]

For this we will use the Fourier transform. The key fact that we need is

[TABLE]

which can be deduced using scaling and symmetry properties of the Fourier transform, or by a direct computation using the Gamma function (see [61, Lemma 1, p.117]). In particular, we are faced with estimating

[TABLE]

As $0\leq s\leq 1$, we can estimate this by

[TABLE]

Using the Lorentz-space555Here $L^{q,\alpha}$ denotes the Lorentz-space defined via the quasi-norm

\|f\|_{L^{q,\alpha}}=\bigl{\|}\lambda|\{|f|>\lambda\}|^{\frac{1}{q}}\bigr{\|}_{L^{\alpha}((0,\infty),\frac{d\lambda}{\lambda})}.

refinements of Young’s inequality and Hölder’s inequality (cf. [30, 55]), we may estimate

[TABLE]

which is acceptable. In the estimates above, we have used $s\leq\tfrac{d}{2}-1$ to guarantee that the exponents above fall into acceptable intervals. ∎

Acknowledgements

Y. Zhang was supported by the US National Science Foundation under grant number DMS-1620465. We are grateful to M. Beceanu for explanations related to the notion of incoming/outgoing waves for radial wave equations. We would also like to thank J. Colliander for introducing us to the related work of Strauss and Vazquez [63], and to an anonymous referee for pointing out the work of [22].