A hyperboloidal method for numerical simulations of multidimensional nonlinear wave equations

We start with the $(n+1)$-dimensional Minkowski metric in spherical polar coordinates,

$$\eta=-\rmd t^{2}+\rmd r^{2}+r^{2}\sigma^{(n-1)},\quad t\in\mathbb{R},\quad r\in(0,\infty).$$

(6)

Here $\sigma^{(n-1)}$ is the standard round metric on $S^{n-1}$, the $(n-1)$-dimensional unit sphere. Explicitly, we have

	$\displaystyle\sigma^{(1)}=\rmd\varphi^{2},$		(7)
	$\displaystyle\sigma^{(2)}=\rmd\theta^{2}+\sin^{2}\theta\,\rmd\varphi^{2},$		(8)

where $\theta\in[0,\pi]$ and $\varphi\in[0,2\pi)$. More generally for $n>3$,

$$\sigma^{(n-1)}=\rmd\theta_{1}^{2}+\sin^{2}\theta_{1}\left(\rmd\theta_{2}^{2}+\sin^{2}\theta_{2}(\ldots(\rmd\theta_{n-2}^{2}+\sin^{2}\theta_{n-2}\rmd\varphi^{2}))\right)$$

(9)

with $\theta_{1},\ldots,\theta_{n-2}\in[0,\pi]$ and $\varphi\in[0,2\pi)$.

We introduce a new time coordinate

$$\tilde{t}=t-\sqrt{a^{2}+r^{2}}$$

(10)

with a constant $a=n/C>0$. The hypersurfaces $\tilde{t}=\mathrm{const}$ are hyperboloids and $C$ is their constant mean curvature (up to a sign as discussed further below). In the new coordinates the Minkowski metric reads

$$\eta=-\rmd\tilde{t}^{2}-\frac{2r}{\sqrt{a^{2}+r^{2}}}\,\rmd\tilde{t}\,\rmd r+\frac{a^{2}}{a^{2}+r^{2}}\,\rmd r^{2}+r^{2}\sigma^{(n-1)}.$$

(11)

Let us also introduce a new radial coordinate $\tilde{r}$ depending on $r$ only. We can arrange for the spatial metric to be conformally flat if we demand

$$\frac{a^{2}}{a^{2}+r^{2}}\rmd r^{2}=\frac{r^{2}}{\tilde{r}^{2}}\rmd\tilde{r}^{2}\;\Leftrightarrow\;\frac{\rmd r}{\rmd\tilde{r}}=\frac{r\sqrt{a^{2}+r^{2}}}{\tilde{r}a}.$$

(12)

This ordinary differential equation can be solved explicitly. If we impose the boundary condition $\tilde{r}\to 1$ as $r\to\infty$, we obtain

$$r=\frac{2a\tilde{r}}{1-\tilde{r}^{2}}=\frac{2n\tilde{r}}{C(1-\tilde{r}^{2})}.$$

(13)

From (10) we observe that on a slice $\tilde{t}=\mathrm{const}$, as the conformal boundary $\tilde{r}\to 1$ is approached, we have

$$t,r\to\infty,\quad t-r\to\tilde{t}=\mathrm{const},$$

(14)

so the slice becomes asymptotically characteristic and $\tilde{r}=1$ corresponds to $\mathrsfs{I}^{+}\,$.

The metric now takes the form

	$\displaystyle\eta$	$\displaystyle=$	$\displaystyle-\rmd\tilde{t}^{2}-\frac{2r^{2}}{\tilde{r}a}\rmd\tilde{t}\,\rmd\tilde{r}+\frac{r^{2}}{\tilde{r}^{2}}\left(\rmd\tilde{r}^{2}+\tilde{r}^{2}\sigma^{(n-1)}\right)$		(15)
		$\displaystyle=$	$\displaystyle\Omega^{-2}\left[-\tilde{\alpha}^{2}\rmd\tilde{t}^{2}+(\rmd\tilde{r}+\tilde{\beta}^{\tilde{r}}\rmd\tilde{t})^{2}+\tilde{r}^{2}\sigma^{(n-1)}\right]=:\Omega^{-2}\tilde{\eta}.$

Here we have introduced the conformal factor

$$\Omega=\frac{\tilde{r}}{r}=\frac{C}{2n}(1-\tilde{r}^{2}),$$

(16)

the conformal shift vector $\tilde{\beta}$, which only has a radial component

$$\tilde{\beta}^{\tilde{r}}=-\frac{\tilde{r}}{a}=-\frac{C}{n}\tilde{r},$$

(17)

and the conformal lapse function

$$\tilde{\alpha}=\frac{C}{2n}(1+\tilde{r}^{2})$$

(18)

so that $\tilde{\alpha}^{2}-(\tilde{\beta}^{\tilde{r}})^{2}=\Omega^{2}$.

From now on we will work entirely with the conformal metric $\tilde{\eta}$. Its induced metric $\tilde{\gamma}$ on the $\tilde{t}=\mathrm{const}$ slices is flat,

$$\tilde{\gamma}=\rmd\tilde{r}^{2}+\tilde{r}^{2}\sigma^{(n-1)}.$$

(19)

The conformal extrinsic curvature of the $\tilde{t}=\mathrm{const}$ slices is¹¹1 We use the sign convention $\mathcal{L}_{\tilde{\nu}}\tilde{\gamma}_{ab}=-2\tilde{K}_{ab}$. For the $\tilde{t}=\mathrm{const}$ slices to intersect $\mathrsfs{I}^{+}\,$we need the constant $C>0$ so that the shift (17) is negative at $\tilde{r}=1$. Spacetime indices $a,b$ are raised and lowered using the conformal metric $\tilde{\eta}$, and repeated indices are summed over.

$$\tilde{K}_{ab}=\frac{\tilde{K}}{n}\tilde{\gamma}_{ab}$$

(20)

with the conformal mean curvature given by

$$\tilde{K}=-\tilde{\alpha}^{-1}C=-\frac{2n}{\tilde{r}^{2}+1}.$$

(21)

Its Lie derivative along the conformal future-directed timelike unit normal $\tilde{\nu}$ to the $\tilde{t}=\mathrm{const}$ slices is

$$\mathcal{L}_{\tilde{\nu}}\tilde{K}=-\tilde{\alpha}^{-1}\tilde{\beta}^{\tilde{r}}\partial_{\tilde{r}}\tilde{K}=\frac{8n\tilde{r}^{2}}{(\tilde{r}^{2}+1)^{3}}.$$

(22)

Combining the contracted Gauss and Ricci equations, we can relate the Ricci scalar $\tilde{R}$ of $\tilde{\eta}$ to the Ricci scalar $\tilde{R}_{\tilde{\gamma}}$ of $\tilde{\gamma}$:

$$\tilde{R}=\tilde{R}_{\tilde{\gamma}}-2\tilde{\alpha}^{-1}\tilde{\Delta}\tilde{\alpha}-2\mathcal{L}_{\tilde{\nu}}\tilde{K}+\tilde{K}^{2}+\tilde{K}_{ab}\tilde{K}^{ab}.$$

(23)

Here

$$\tilde{\Delta}=\partial_{\tilde{r}}^{2}+(n-1){\tilde{r}}^{-1}\partial_{\tilde{r}}+\tilde{r}^{-2}\mathring{\Delta}^{(n-1)}$$

(24)

is the standard flat Laplacian of the flat conformal metric $\gamma$ and $\mathring{\Delta}^{(n-1)}$ is the Laplace-Beltrami operator on the sphere $S^{n-1}$. Explicitly, we have

	$\displaystyle\mathring{\Delta}^{(1)}=\partial_{\varphi}^{2},$		(25)
	$\displaystyle\mathring{\Delta}^{(2)}=\partial_{\theta}^{2}+\cot\theta\,\partial_{\theta}+\frac{1}{\sin^{2}\theta}\partial_{\varphi}^{2}.$		(26)

In our numerical implementation for higher dimensions $n>3$, we shall impose an SO$(n-1)$ symmetry so that all but one of the angular coordinates are suppressed. The angular Laplacian then reduces to

$$\mathring{\Delta}^{(n-1)}=\partial_{\theta}^{2}+(n-2)\cot\theta\,\partial_{\theta}.$$

(27)

Using the fact that $\tilde{R}_{\tilde{\gamma}}=0$ and the form of the conformal extrinsic curvature (20), we obtain from (23)

$$\tilde{R}=-2\tilde{\alpha}^{-1}\tilde{\Delta}\tilde{\alpha}-2\mathcal{L}_{\tilde{\nu}}\tilde{K}+\frac{n+1}{n}\tilde{K}^{2}=\frac{4n[-\tilde{r}^{4}+(n-5)\tilde{r}^{2}+n]}{(\tilde{r}^{2}+1)^{3}}.$$

(28)

The scalar field $\Phi$ is supposed to obey the NLW

$$\Box\Phi=\mu|\Phi|^{p-1}\Phi,$$

(29)

where $\Box$ is the d’Alembertian of the Minkowski metric $\eta$.

Defining a conformally rescaled scalar field

$$\tilde{\Phi}=\Omega^{(1-n)/2}\Phi$$

(30)

and using the conformal transformation properties of the d’Alembertian and the Ricci scalar, we have

$$\tilde{\Box}\tilde{\Phi}-\frac{n-1}{4n}\tilde{R}\tilde{\Phi}=\Omega^{-(n+3)/2}\left(\Box\Phi-\frac{n-1}{4n}R\Phi\right)=\mu\Omega^{[p(n-1)-n-3]/2}|\tilde{\Phi}|^{p-1}\tilde{\Phi},$$

(31)

where on the right-hand side we have used (29) and the fact that $\eta$ is flat, $R=0$. This equation is regular at $\mathrsfs{I}^{+}\,$, where $\Omega=0$, iff $n\geqslant 2$ and

$$p(n-1)-n-3\geqslant 0\;\Leftrightarrow\;p\geqslant\frac{n+3}{n-1}=:p_{\mathrm{conf}}.$$

(32)

For any $n\geqslant 3$ we have $p_{\mathrm{conf}}<p_{\mathrm{crit}}$, where $p_{\mathrm{crit}}$ was defined in (4), cf. table 1. Thus the conformal approach is capable of treating both energy-subcritical and supercritical NLWs.

Performing an $(n+1)$-decomposition of the conformal d’Alembertian, (31) becomes

$$-\mathcal{L}_{\tilde{\nu}}^{2}\tilde{\Phi}+\tilde{\Delta}\tilde{\Phi}+\tilde{K}\mathcal{L}_{\tilde{\nu}}\tilde{\Phi}+\tilde{\alpha}^{-1}\tilde{\alpha}^{,i}\tilde{\Phi}_{,i}-\frac{n-1}{4n}\tilde{R}\tilde{\Phi}=\mu\Omega^{[p(n-1)-n-3]/2}|\tilde{\Phi}|^{p-1}\tilde{\Phi},$$

(33)

where spatial indices $i$ are raised and lowered using the flat conformal metric $\tilde{\gamma}$. We can write this equation in first-order in time form by introducing an auxiliary field $\tilde{\Pi}$,

	$\displaystyle\mathcal{L}_{\tilde{\nu}}\tilde{\Phi}$	$\displaystyle=:$	$\displaystyle\tilde{\Pi},$		(34)
	$\displaystyle\mathcal{L}_{\tilde{\nu}}\tilde{\Pi}$	$\displaystyle=$	$\displaystyle\tilde{\Delta}\tilde{\Phi}+\tilde{K}\tilde{\Pi}+\tilde{\alpha}^{-1}\tilde{\alpha}^{,i}\tilde{\Phi}_{,i}-\frac{n-1}{4n}\tilde{R}\tilde{\Phi}-\mu\Omega^{[p(n-1)-n-3]/2}|\tilde{\Phi}|^{p-1}\tilde{\Phi}.$		(35)

Writing out the Lie derivative $\mathcal{L}_{\tilde{\nu}}=\tilde{\alpha}^{-1}(\partial_{\tilde{t}}-\mathcal{L}_{\tilde{\beta}})$, we arrive at

	$\displaystyle\tilde{\Phi}_{,\tilde{t}}=\tilde{\beta}^{\tilde{r}}\tilde{\Phi}_{,\tilde{r}}+\tilde{\alpha}\tilde{\Pi},$		(36)
	$\displaystyle\tilde{\Pi}_{,\tilde{t}}=\tilde{r}^{1-n}\left[\tilde{r}^{n-1}\left(\tilde{\beta}^{\tilde{r}}\tilde{\Pi}+\tilde{\alpha}\tilde{\Phi}_{,\tilde{r}}\right)\right]_{,\tilde{r}}+\tilde{\alpha}\tilde{r}^{-2}\mathring{\Delta}^{(n-2)}\tilde{\Phi}-\frac{n-1}{4n}\tilde{\alpha}\tilde{R}\tilde{\Phi}$
	$\displaystyle\qquad-\mu\tilde{\alpha}\Omega^{[p(n-1)-n-3]/2}|\tilde{\Phi}|^{p-1}\tilde{\Phi}.$		(37)

We have double-checked using computer algebra that (36)–(37) are equivalent to (29) if the coordinate transformation given by (10) and (13) is applied directly.

Equations (36)–(37) are perfectly regular at future null infinity $\tilde{r}=1$, where $\Omega=0$, provided that $p\geqslant p_{\mathrm{conf}}$ (or $\mu=0$). No boundary conditions are needed there because all characteristics point towards the exterior of the domain.

It should be noted that we have written the radial principal part of (37) in conservative form. The finite differencing (section 3.1) will be applied in precisely this order. This was found to be essential for the numerical evolutions to be stable.

We use a hybrid method for the spatial discretisation, which combines a finite-difference method in the radial direction with a pseudo-spectral method in the angular directions. The latter has the advantage that it is easier to build in parity conditions by an appropriate choice of expansion functions, as will be explained below. Furthermore, if we introduce a logically rectangular grid with equidistant grid points in the standard spherical coordinates $\theta$ and $\varphi$ on the two-sphere, the distance between neighbouring $\varphi$-grid points becomes very small near the poles. This leads to a severe restriction on the time step due to the Courant-Friedrichs-Lewy (CFL) condition if a finite-difference method is used. With a pseudo-spectral method, it is straightforward to remove an increasing number of higher (more oscillatory) $\varphi$-modes as the axis is approached (section 3.2), thus allowing for larger time steps.

We take the spatial dimension to be $n=3$ without symmetries first. Higher dimensions with $\mathrm{SO}(n-1)$ symmetry will be treated as a special case.

We introduce equidistant grid points in each dimension,

	$\displaystyle\tilde{r}_{i}=(i+\textstyle\frac{1}{2})h_{\tilde{r}},\quad\;i=0,1,\ldots,N_{\tilde{r}}-1,\quad h_{\tilde{r}}=1/(N_{\tilde{r}}-\textstyle\frac{1}{2}),$		(38)
	$\displaystyle\theta_{j}=(j+\textstyle\frac{1}{2})h_{\theta},\quad j=0,1,\ldots,N_{\theta}-1,\quad h_{\theta}=\pi/N_{\theta},$		(39)
	$\displaystyle\varphi_{k}=kh_{\varphi},\qquad\;\;k=0,1,\ldots,N_{\varphi}-1,\quad h_{\varphi}=2\pi/N_{\varphi}.$		(40)

The radial grid is staggered at the origin, whereas the outermost grid point coincides with $\mathrsfs{I}^{+}\,$($\tilde{r}=1)$. The grid in the elevation $\theta$ is staggered about the axis. The reason for this choice is that there are terms in the wave equation that are formally singular at the origin ($\tilde{r}=0$) and on the axis ($\sin\theta=0$), which would be difficult to evaluate if there were grid points at the origin and on the axis.

Fourth-order finite differences are used in order to approximate the radial derivatives. In the interior a centred stencil is used. Setting $u_{i}:=u(\tilde{r}_{i})$ for a smooth function $u$ and leaving out its angular dependence for simplicity,

	$\displaystyle u^{\prime}_{i}\approx(u_{i-2}-8u_{i-1}+8u_{i+1}-u_{i+2})/(12h_{\tilde{r}}),$		(41)
	$\displaystyle u^{\prime\prime}_{i}\approx(-u_{i-2}+16u_{i-1}-30u_{i}+16u_{i+1}-u_{i+2})/(12h_{\tilde{r}}^{2})$		(42)

for $i=0,1,\ldots,N_{\tilde{r}}-3$, where $\approx$ stands for equality up to a truncation error of the order $O(h_{\tilde{r}}^{4})$. Evaluating the above stencils at $i=0$ and $i=1$ requires values of $u$ at the ghost points $\tilde{r}_{-1}=-h_{\tilde{r}}/2$ and $\tilde{r}_{-2}=-3h_{\tilde{r}}/2$, which we copy from interior points according to

$$u(-\tilde{r},\theta,\varphi)=u(\tilde{r},\pi-\theta,\pi+\varphi),$$

(43)

satisfied by any smooth function $u$. We require $N_{\varphi}$ to be even, which ensures that if there is a grid point at $\varphi$, there is also a grid point at $\pi+\varphi$.

Near the outer boundary at $\tilde{r}=1$, one-sided differences are used,

	$\displaystyle u^{\prime}_{N_{\tilde{r}}-2}\approx(-u_{N_{\tilde{r}}-5}+6u_{N_{\tilde{r}}-4}-18u_{N_{\tilde{r}}-3}+10u_{N_{\tilde{r}}-2}+3u_{N_{\tilde{r}}-1})/(12h_{\tilde{r}}),$		(44)
	$\displaystyle u^{\prime}_{N_{\tilde{r}}-1}\approx(3u_{N_{\tilde{r}}-5}-16u_{N_{\tilde{r}}-4}+36u_{N_{\tilde{r}}-3}-48u_{N_{\tilde{r}}-2}+25u_{N_{\tilde{r}}-1})/(12h_{\tilde{r}}),$		(45)
	$\displaystyle u^{\prime\prime}_{N_{\tilde{r}}-2}\approx(u_{N_{\tilde{r}}-6}-6u_{N_{\tilde{r}}-5}+14u_{N_{\tilde{r}}-4}-4u_{N_{\tilde{r}}-3}-15u_{N_{\tilde{r}}-2}+10u_{N_{\tilde{r}}-1})/(12h_{\tilde{r}}^{2}),$		(46)
	$\displaystyle u^{\prime\prime}_{N_{\tilde{r}}-1}\approx(-10u_{N_{\tilde{r}}-6}+61u_{N_{\tilde{r}}-5}-156u_{N_{\tilde{r}}-4}+214u_{N_{\tilde{r}}-3}-154u_{N_{\tilde{r}}-2}+45u_{N_{\tilde{r}}-1})$
	$\displaystyle\qquad\qquad/(12h_{\tilde{r}}^{2}),$		(47)

the truncation error again being $O(h_{\tilde{r}}^{4})$.

Our pseudo-spectral method is based on Fourier expansions in both $\theta$ and $\varphi$. This allows us to employ fast Fourier transform (FFT) techniques [14], which would not be available if we expanded the fields in spherical harmonics.

Leaving out its radial dependence for simplicity, a smooth function $u$ is expanded as

$$u(\theta,\varphi)\approx\sum_{l=0}^{N_{\theta}-1}\left(\cos(l\theta)\sum_{m=0\atop{m\,\mathrm{even}}}^{N_{\varphi}/2-1}a_{lm}\,\rme^{\rmi m\varphi}+\sin(l\theta)\sum_{m=1\atop{m\,\mathrm{odd}}}^{N_{\varphi}/2-1}a_{lm}\,\rme^{\rmi m\varphi}\right).$$

(48)

Any smooth function $u$ on the two-sphere must obey

$$u(-\theta,\varphi)=u(\theta,\pi+\varphi),\qquad u(\pi+\theta,\varphi)=u(\pi-\theta,\pi+\varphi),$$

(49)

and the way we expand the even and odd $m$-modes separately in (48) enforces this. For $u$ to be real, the coefficients $a_{lm}$ in (48) must be complex, and we may use

$$\mathrm{Re}(a_{lm}\rme^{\rmi m\varphi})=\mathrm{Re}\,a_{lm}\cos(m\varphi)-\mathrm{Im}\,a_{lm}\sin(m\varphi).$$

(50)

In order to compute the expansion coefficients $a_{lm}$ from the point values $u_{jk}:=u(\theta_{j},\varphi_{k})$, we first perform a real FFT w.r.t. $\varphi$ using scipy.fft.rfft. Next, we perform a discrete cosine transform w.r.t. $\theta$ for the even-$m$ modes (scipy.fft.dct) and a discrete sine transform for the odd-$m$ modes (scipy.fft.dst). The (default) Type II versions of the discrete cosine and sine transforms implemented in scipy are appropriate for our staggered $\theta$ grid.

Derivatives of the expansion can now be computed analytically by simply differentiating the basis functions $\cos(l\theta)$, $\sin(l\theta)$ and $\rme^{\rmi m\varphi}$. Finally in order to obtain the values of the derivatives at the grid points $(\theta_{j},\varphi_{k})$, we apply the inverse discrete cosine and sine transforms w.r.t. $\theta$, followed by an inverse real FFT w.r.t. $\varphi$.

In the code the main object is the array of values of the unknowns $\tilde{\Phi}$ and $\tilde{\Pi}$ at the grid points $(\tilde{r}_{i},\theta_{j},\varphi_{k})$. Terms in the evolution equations, especially nonlinear terms, are evaluated and added pointwise.

We remark that a purely spectral (rather than pseudo-spectral collocation) method such as the one used in [10] cannot be employed in our case because of the nonlinear term in the wave equation we consider.

We employ the method of lines, whereby the evolution equations (36)–(37) are first discretised in space as described in the previous subsection, leading to a system of ordinary differential equations (ODEs), one for each grid point. These ODEs are then integrated forward in time using the standard fourth-order Runge-Kutta method.

In order to obtain a stable finite-difference method, fifth-order Kreiss-Oliger dissipation [15] is applied in the radial direction, whereby the term

	$\displaystyle(Qu)_{i}=\frac{\varepsilon}{64h_{\tilde{r}}}(u_{i-3}-6u_{i-2}+15u_{i-1}-20u_{i}+15u_{i+1}-6u_{i+2}+u_{i+3})$
	$\displaystyle\qquad\;=\frac{\varepsilon}{64}h_{\tilde{r}}^{5}\left(\frac{\rmd^{6}u}{\rmd\tilde{r}^{6}}\right)_{i}+\Or(h_{\tilde{r}}^{7})$		(53)

for $i=0,1,\ldots,N_{\tilde{r}}-4$ is added to the right-hand side of the discretised evolution equations for both $\tilde{\Phi}$ and $\tilde{\Pi}$. Ghost points near the origin are filled according to (43). No dissipation is added at the outermost three radial grid points. We typically take the parameter $\varepsilon=0.2$. It should be noted that the artificial dissipation term (53) is of higher order in $h_{\tilde{r}}$ than the truncation error of the finite-difference method.

In a pseudo-spectral method, nonlinear terms such as the one in our wave equation can cause high frequencies to be represented incorrectly due to aliasing [16]. We address this by applying spectral filtering according to the 2/3 rule, whereby the top third of all $\theta$-modes are removed at each substep of the Runge-Kutta algorithm.

In order to combat the clustering of grid points near the poles of the two-sphere, we filter out an increasing number of the highest Fourier modes w.r.t. $\varphi$ as $\theta$ approaches $0$ or $\pi$. Specifically, we remove the proportion $1-\sin\theta$ of all $\varphi$-Fourier modes (beginning at the highest frequencies) [17].

The time step is restricted by the smallest distance between neighbouring grid points, and this occurs between neighbouring $\theta$-grid points at the smallest radius: $\Delta x_{\min}=\tilde{r}_{0}h_{\theta}$, where $\tilde{r}_{0}=h_{\tilde{r}}/2$. The time step is then taken to be $\Delta\tilde{t}=\lambda\Delta x_{\min}$, where $0<\lambda<1$ according to the CFL condition. We typically choose $\lambda=0.8$.

For expressions such as the energy and flux (section 4.2) we will need to compute integrals on a $\tilde{t}=\mathrm{const}$ slice and on its spherical boundary. The radial integration is performed using Simpson’s rule (taking $N_{\tilde{r}}$ to be even),

	$\displaystyle\int_{0}^{1}u(\tilde{r})\rmd\tilde{r}\approx\frac{h_{\tilde{r}}}{3}\left(\frac{27}{8}u_{0}+\frac{17}{8}u_{1}+4u_{2}+2u_{3}+4u_{4}+2u_{5}+\ldots\right.$
	$\displaystyle\left.+4u_{N_{\tilde{r}}-4}+2u_{N_{\tilde{r}}-3}+4u_{N_{\tilde{r}}-2}+u_{N_{\tilde{r}}-1}\right),$		(54)

where the error is of the same order as the truncation error of the finite-difference scheme, $\Or(h_{\tilde{r}}^{4})$. The modified coefficients of $u_{0}$ and $u_{1}$ are due to the fact that the grid is staggered at the origin, with a first grid point at $\tilde{r}_{0}=h_{\tilde{r}}/2$, and $u$ is assumed to be an even function of $\tilde{r}$ here.

The spectral expansion on the right-hand side of (48) can easily be integrated exactly:

$$\int_{S^{2}}u\,dS^{(2)}=\int_{0}^{2\pi}\int_{0}^{\pi}u(\theta,\varphi)\,\sin\theta\,\rmd\theta\,\rmd\varphi\approx 2\pi\sum_{l=0\atop{l\,\mathrm{even}}}^{N_{\theta}-1}\frac{2a_{l0}}{1-l^{2}}.$$

(55)

Given the complexity of the algebraic operations involved in deriving the evolution equations, it is very useful to have an exact solution at hand in order to check that the implementation is correct, and that the numerical method converges as expected.

In A we construct exact solutions to the linear wave equation. They are based on an expansion in spherical harmonics and are given in terms of a mode function, which we take to be

$$F(x)=A\,x\,\exp\left[-\frac{1}{2}\left(\frac{x}{\sigma}\right)^{2}\right].$$

(59)

We choose $A=1$ and $\sigma=1$.

For the convergence test in the case of $n=3$ without symmetries, we consider a superposition of two modes with spherical harmonic indices $(l,m)=(1,1)$ and $(2,-2)$. In the case of $n=5$ with SO$(4)$ symmetry, we take two modes with $l=1$ and $2$.

The initial data are set according to $\tilde{\Phi}_{0}:=\tilde{\Phi}(\tilde{t}=0)=\tilde{\Phi}_{\mathrm{exact}}(\tilde{t}=\tilde{t}_{0})$ and similarly for $\tilde{\Pi}$, where we choose $\tilde{t}_{0}=-15$ so that the solution is initially ingoing and centred roughly in the middle of the conformal radial interval $[0,1]$.

We now evolve these initial data and compare with the exact solution. In figure 1 we plot the $L^{2}$ norm of the error

$$\|\tilde{\Phi}-\tilde{\Phi}_{\mathrm{exact}}\|_{L^{2}}=\left[\int_{0}^{1}\tilde{r}^{n-1}\rmd\tilde{r}\int_{S^{n-1}}(\tilde{\Phi}-\tilde{\Phi}_{\mathrm{exact}})^{2}\,dS^{(n-1)}\right]^{1/2}$$

(60)

as a function of time for three different radial resolutions at fixed angular resolution. We observe that doubling the resolution decreases the error nearly by a factor of $2^{4}=16$, as expected for a fourth-order accurate finite-difference method.

We remark that at the chosen angular resolution, the spherical harmonics are represented exactly by our pseudo-spectral method, so there is no discretisation error in the angular dimensions for these linear evolutions.

In this section we derive and verify numerically the energy balance that the scalar field obeys on our hyperboloidal foliation. On standard slices of Minkowski spacetime approaching spacelike infinity, the energy (3) is conserved. In contrast, the energy computed on hyperboloidal slices decreases during the evolution. This is due to the fact that hyperboloidal slices intersect $\mathrsfs{I}^{+}\,$, so that outgoing waves pass through this conformal boundary, carrying away energy. We will compute this energy flux explicitly. A similar energy balance was derived in [10] in the context of linear scalar fields in Kerr spacetime, where an additional inner boundary arises at the black hole horizon.

The scalar field $\Phi$ is associated with an energy-momentum tensor

$$T_{ab}=\nabla_{a}\Phi\nabla_{b}\Phi-\textstyle\frac{1}{2}\eta_{ab}\nabla_{c}\Phi\nabla^{c}\Phi-\frac{\mu}{p+1}|\Phi|^{p+1}\eta_{ab},$$

(61)

where $\nabla$ is the covariant derivative of the Minkowski metric $\eta$. The vanishing of its divergence, $\nabla^{b}T_{ab}=0$, is equivalent with the NLW (29). The Killing vector $k:=\partial/\partial t=\partial/\partial\tilde{t}$ gives rise to the conserved current

$$E^{a}=T^{a}{}_{b}k^{b}=T^{a}{}_{t}$$

(62)

satisfying

$$\nabla_{a}E^{a}=0.$$

(63)

Let $\Sigma_{\tilde{t}}$ denote the slice of constant hyperboloidal time $\tilde{t}$ with future-directed unit normal $\nu$ (normalised w.r.t. $\eta$) and induced Riemannian metric $\gamma_{ab}=\eta_{ab}+\nu_{a}\nu_{b}$. Furthermore let $s$ denote the outward-directed normal to the timelike surface $r=\mathrm{const}$ (again normalised w.r.t. $\eta$) and $h_{ab}=\eta_{ab}-s_{a}s_{b}$ the induced pseudo-Riemannian metric on such a surface. Applying Gauss’ theorem to (63), the energy balance reads

	$\displaystyle E(\tilde{t}_{2})-E(\tilde{t}_{1}):=\int_{\Sigma_{\tilde{t}_{2}}}\nu_{a}E^{a}\sqrt{\det\gamma}\,\rmd r\,\rmd^{n-1}\theta-\int_{\Sigma_{\tilde{t}_{1}}}\nu_{a}E^{a}\sqrt{\det\gamma}\,\rmd r\,\rmd^{n-1}\theta$		(64)
	$\displaystyle\qquad=\int_{\tilde{t}_{1}}^{\tilde{t}_{2}}\int_{S^{n-1}}s_{a}E^{a}\sqrt{\det|h|}\,\rmd^{n-1}{\bf\theta}\,\rmd\tilde{t}\,\Big{|}_{\mathrsfs{I}^{+}}=:F(\tilde{t}_{1},\tilde{t}_{2}),$		(65)

where $\theta$ stands collectively for the $n-1$ angular coordinates on $S^{n-1}$.

Expressing $\Phi$ in terms of the conformally rescaled scalar field $\tilde{\Phi}$ defined in 30, we obtain for the energy

	$\displaystyle E(\tilde{t})=\frac{C}{4n}\int_{0}^{1}\tilde{r}^{n-1}\rmd\tilde{r}\int_{S^{n-1}}dS^{(n-1)}\left\{(1+\tilde{r}^{2})\left[\tilde{\Pi}^{2}+\tilde{\Phi}_{,\tilde{r}}^{2}+\frac{1}{\tilde{r}^{2}}\|\mathring{\nabla}^{(n-1)}\tilde{\Phi}\|^{2}\right.\right.$
	$\displaystyle\left.\qquad+2\mu\Omega^{[p(n-1)-n-3]/2}\frac{1}{p+1}|\tilde{\Phi}|^{p+1}\right]$
	$\displaystyle\qquad\left.+2(n-1)\frac{\tilde{r}^{2}-1}{\tilde{r}^{2}+1}\,\tilde{r}\tilde{\Phi}_{,\tilde{r}}\tilde{\Phi}-4r\tilde{\Phi}_{,\tilde{r}}\tilde{\Pi}+(n-1)^{2}\frac{\tilde{r}^{2}}{\tilde{r}^{2}+1}\tilde{\Phi}^{2}\right\}.$		(66)

Here $\mathring{\nabla}^{(n-1)}$ refers to the gradient on $S^{n-1}$ and $dS^{(n-1)}=\sqrt{\det\sigma^{(n-1)}}\rmd^{n-1}\theta$ to its area element, e.g. for $n=3$

$$\|\mathring{\nabla}^{(2)}\Phi\|^{2}=\Phi_{,\theta}^{2}+\frac{1}{\sin^{2}\theta}\,\Phi_{,\varphi}^{2},\qquad dS^{(2)}=\sin\theta\,\rmd\theta\,\rmd\varphi.$$

(67)

In the first line of (4.2), we recognize the usual energy expression for the linear scalar field. The second line of (4.2) contains the contribution from the nonlinear potential, and the terms in the third line arise from the conformal rescaling.

The flux integral takes the form

$$F(\tilde{t}_{1},\tilde{t}_{2})=-\frac{C^{2}}{n^{2}}\int_{\tilde{t}_{1}}^{\tilde{t}_{2}}\rmd\tilde{t}\int_{S^{(n-1)}}dS^{(n-1)}(\tilde{\Phi}_{,\tilde{r}}-\tilde{\Pi})^{2}\,\Big{|}_{\tilde{r}=1}.$$

(68)

It is is manifestly negative, reflecting the fact that the waves carry away energy as they pass through $\mathrsfs{I}^{+}\,$.

We now verify the energy balance numerically for two cases: $n=3,\,p=5$ (critical) and $n=5,\,p=3$ (supercritical). In both cases we impose an SO$(n-1)$ symmetry. The initial data are taken to be momentarily static with

$$\tilde{\Phi}_{0}=A\exp\left[-\left(\frac{\tilde{r}-\tilde{r}_{0}}{\sigma}\right)^{2}\right]\,Y_{l}(\theta),$$

(69)

and $\tilde{\Phi}_{,t}=0$ at $t=0$ implies

$$\tilde{\Pi}_{0}=\frac{2\tilde{r}}{1+\tilde{r}^{2}}\tilde{\Phi}_{0,\tilde{r}}.$$

(70)

We superpose two modes with $l=2$ and $l=3$ and choose the parameters $\tilde{r}_{0}=0.3$ and $\sigma=0.07$. The amplitude for both modes is taken to be $A=12$ for $n=3,\,p=5$ and $A=200$ for $n=5,\,p=3$, which is smaller than but close to the critical amplitude beyond which blow-up occurs in the focusing case ($\mu=-1$). The numerical resolution is taken to be $N_{\tilde{r}}=4000,\,N_{\theta}=16$.

Figure 2 demonstrates that the energy balance

$$E(\tilde{t})-F(0,\tilde{t})=E(0)$$

(71)

is well satisfied numerically, i.e., the sum of the numerically computed terms on the left-hand side is indeed nearly constant. For an evolution with a defocusing nonlinearity ($\mu=1$), the initial energy is higher than for a focusing nonlinearity ($\mu=-1$) due to the different sign of the potential energy term in (4.2).