The d’Alembert solution in hyperboloidal foliations

. We explicitly construct the analogue of the d’Alembert solution to the 1+1 wave equation in an hyperboloidal setting. This hyperboloidal d’Alembert solution is used, in turn, to gain intuition into the behaviour of solutions to the wave equation in a hyperboloidal foliation and to explain an apparently anomalous permanent displacement of the solution in numerical simulations discussed in the literature


Introduction
A hyperboloidal hypersurface in a spacetime (hyperboloid for short) is a spacelike hypersurface which intersects null infinity at cuts having the topology of the 2-sphere -see e.g.[1,2].Prime examples of hyperboloids are constant mean curvature (CMC) hypersurfaces in an asymptotically flat spacetime -although, of course, these are not the only possibilities.Hyperboloidal foliations were considered in the seminal work by H. Friedrich on the semiglobal stability of the Minkowski spacetime [3,4].More recently, hyperboloidal foliations have been considered as a natural way of providing gauges for the numerical evolution of the equations of linearised gravity in the context of the study of perturbations of black hole spacetimes -see e.g.[5,6,7].As discussed in [8], the use of hyperboloidal foliations in the study of wave equations avoids the use of unnatural boundary conditions by providing a built-in outflow boundary condition.Alternative approaches to boundary conditions in the context of numerical simulations in General Relativity can be found in e.g.[9,10,11,12,13].
An important step when working with a new type of gauge conditions, is to develop an intuition on the behaviour of some basic solutions.In the case of the wave equation in Cartesian coordinates, the basic intuition is provided by the d'Alembert solution (in the 1+1 case), the Poisson-Hadamard's solution (in the 1+2-dimensional case) and the Kirkhoff solution in the 1+3-dimensional setting.These formulae allow to compute the solutions to the wave equation in the whole of space in terms of the prescribed initial data -see e.g.[14,15].
In the present article we restrict ourselves, for conciseness, to the discussion of the 1+1 dimensional wave equation.The more physically relevant Kirkhoff solution for the 1+3-dimensional wave equation can be recovered from the d'Alembert solution via the method of spherical decomposition means -see e.g.[14].Moreover, the 1+1dimensional wave equations (with a potential) naturally arise from a decomposition in spherical harmonics of the solutions to the 1+3-dimensional wave equation.
Key insights provided by d'Alembert's solution into the behaviour of solutions to the 1+1-dimensional wave equation include the fact that the support of the solution on hypersurfaces of constant time (i.e. the regions where the solution is non-zero) propagates with finite speed.Moreover, if the solution had initially compact support and the initial value of its time derivative vanishes then after the front of the wave has passed, the solution returns to the trivial value (zero) -i.e. the solution shows no permanent displacements.
The phenomenology described in the previous paragraph is now to be contrasted with what is observed in the numerical evaluation of solutions to the 1+1 dimensional wave equation in the hyperboloidal setting [16,17,18].Initial configurations with initial data of essentially compact support and vanishing time derivative do not return to the zero value after the wave front has passed.Instead, the solution settles to a constant value as it will be demonstrated by our numerical simulations implementing numerical methods developed in previous work [18,19,20,21].In the following we refer to this behaviour as permanent displacement.When contrasted with the Cartesian case this behaviour may appear, in first instance, puzzling.In this article we construct the analogue, for hyperboloidal foliations, of the d'Alembert solution and use it to explain the above mentioned behaviour in the numerical solutions.Further insights into the behaviour of solutions to the wave equation in the hyperboloidal wave equation can be obtained by considering various choices of initial data.Permanent displacements can also be observed in solutions to the Cartesian wave equation by choosing appropriately the initial value of the time derivative of the scalar field -in this particular case this piece of the data appears in a term involving an integral.The novelty in the hyperboloidal setting is that the permanent displacement is produced by integral terms involving the actual value of the scalar field.Our hyperboloidal D'Alembert solution allows, in principle, to fine-tune data so as to minimise this effect.
Solutions in higher spatial dimensions, in both Cartesian and hyperboloidal settings, can be analysed by means of the method of spherical means [14].This method relies on d'Alembert's formula to construct the explicit solutions to the equations.As such, analogous permanent displacements can be observed in the auxiliary variables to which d'Alembert's method is applied.However, given that the solutions in higher spatial dimensions decay (no decay occurs in 1+1 dimensions) the displacements are harder to observe in the simulations.This plot shows two pulses of half the size of the original bump moving in opposite directions.The support of the solution remains, at all times, compact.Bottom: contour plot of the d'Alembert solution for the choice f = 0 and g a bump.This choice of initial conditions generates a permanent displacement in the value of the scalar field in a compact domain.

The 1 + 1-dimensional wave equation in Cartesian coordinates
In this section we make some remarks regarding the properties of solutions to the wave equation in the 1 + 1 and the spherically symmetric 1 + 3 settings.

The d'Alembert solution
We start by looking at the initial value problem for the 1 + 1-dimensional wave equation on the real line.That is, we want to study the problem where f ∈ C 2 (R) and g ∈ C 1 (R) are two specifiable functions over R. Furthermore, we note here and in the rest of the manuscript we work with geometric units and thus c = 1.
It is well-known that the general solution to the 1 + 1-dimensional wave equation is of the form where F , G are two twice differentiable functions of a single argument.For the initial value problem (1a)-(1c) one can express the functions F and G in terms of the initial value problem data f , g.This gives rise to d'Alembert's solution: Intuition regarding formula (3) can be obtained from looking at a number of particular cases.Setting g = 0 and choosing f to have, say the shape of a bump, gives the well-know picture of the initial bump splitting into two smaller bumps of half the size each one spreading in the opposite direction with velocity 1.As we are in a 1 + 1 setting the solution does not decay in time.If the function f is of compact support then for t > 0 the solution has also a compact support.However, the support will now consist of two disconnected parts -see Figure 1, top.
An alternative situation, rarely discussed, is to set f = 0 and set g to be, again, a bump function.A particularly convenient choice of a bump function is which can be easily integrated.In this case d'Alembert's formula gives the solution Observe that for a fixed x, one has that ϕ(t, x) → 1 as t → ∞.Consequently, the support of ϕ grows as t increases.This situation is illustrated in Figure 1, bottom.
Remark.In view of the linearity of the wave equation a generic solution is a combination of the two effects discussed previously.

The wave equation in the hyperboloidal setting
Following [8] we introduce coordinates τ and ρ via where the height function h is given by i Figure 3: 10 in the height function h as given by equation ( 4).The hyperboloids are hypersurfaces of constant mean curvature given by K = 3/ √ S -see [1].The Penrose diagram is quantitatively correct in the sense that it is the plot of the height function composed with the transformation implementing the conformal embedding of the 1+1 dimensional Minkowski spacetime.
where S is a positive number and It can be verified that the hypersurfaces of constant coordinate τ are hyperboloids -i.e.spacelike hypersurfaces reaching to future null infinity; see Figure 2. The specific choice S = 1 will be used in the numerical experiments discussed in this article -see Figure 3.
The endpoints of the hyperboloid correspond to the sets for which ρ = ±S.In particular, the 1 + 1 wave equation takes the form where now ϕ = ϕ(τ, ρ).The above equation evaluated at ρ = ±S takes the form From the above factorisation it follows that one of the characteristic speeds of the system vanishes (associated to the operator ∂ τ ) so that the wave equation behaves, at the boundary, as an advection equation (the operator ∂ τ ∓ ∂ ρ ) propagating information out of the system.Accordingly, the boundaries are at ρ = ±S are of an outflow nature and, thus, no boundary conditions are required.We will consider the initial value problem for equation ( 5) on the hyperboloid with τ = 0.For this we prescribe, as usual where f and g are now thought of as functions of ρ ∈ [−S, S] -not to be confused with the functions with the same name used in (1b)-(1c).
Remark 1.In the following we will show how to adapt d'Alembert's method to obtain a formula for the solution in terms of the initial conditions.

Construction of the solution in the hyperboidal setting
As already discussed, the general solution to the 1+1 wave equation on the Minkowski spacetime is given by formula (2).Making use of the transformation formulae between the Cartesian coordinates (t, x) and the hyperboloidal coordinates (τ, ρ), one has that where we have set so that the general solution to the wave equation ( 5) is given by Observe that ζ → ∞ as ρ → S and χ → ∞ as ρ → −S while ζ(−S) = χ(S) = 0.
A calculation then readily shows that where ′ denotes differentiation with respect to the argument of the function.
Differentiating equation (10a) with respect to ρ and taking into account the initial conditions (7a)-(7b) one obtains the system of equations The above is a linear system for Its solution can be readily be found to be given by Now, re-expressing (11a) in terms of ζ as given by ( 8) one has that Integration of this equation yields Using integration by parts one can remove the derivative in the second integral so that S + s ds.
Now, using (9) one finds, after some manipulations, that S + s ds.
Combining the above expressions using the replacements This is the main result in this article.It provides the solution to the wave equation in terms of the initial values of ϕ and ∂ τ ϕ on the hyperboloid given by the condition τ = 0.
Remark 2. It can be readily be verified using the relation . Moreover, given that equation ( 12) is manifestly of the form (9) allows to verify that it provides the unique solution to the initial value problem for equation (5) with data given by (7a)-(7b).
Remark 3. Equation ( 12) is consistent with the fact that the wave equation ( 5) requires no boundary conditions at ρ = ±S -i.e. that one has outflow boundary conditions.In other words, equation (12) does not see the boundaries.

Negative permanent displacements
A peculiarity of expression (12) is the integral involving f -in the third line of the equation.For ease of reference let with the understanding that χ and ζ are functions of (τ, ρ) via the transformation (8).
As τ ≥ 0, it follows that the lower integration limit in the above integral is always positive for fixed χ if τ is sufficiently large.
The effect of I(τ, ρ) can be better appreciated by choosing data such that g(ρ) = 0.In this case, if f is chosen to be in the form of a non-negative "bump" (of say, compact support) then, for fixed ρ, provides a negative, monotonously decreasing contribution to the value of ϕ(τ, ρ).In particular, one has that lim The integral is finite if f is of compact support.Moreover, one has that its value is independent of the value of ρ.In fact, under the above assumptions on f one has that which vanishes for a function f of support strictly contained in [−S, S].Expression (13) explains the results observed in the numerical simulations of [16,17,18].That is, it shows that the solution exhibits a permanent displacement for late times.A particular example of this phenomenon will be discussed in the next section.

Numerical experiments
We will now verify our theoretical results numerically by solving equation (1a) in both the original (t, x) Cartesian coordinate chart and the hyperboloidal slice as defined by  4).Initially the gaussian pulse has maximum amplitude (highlighted in red) it then splits with equal amplitudes as the simulation runs its course.Outflow is observed, and by the end of the simulation only numerical residue is detected as expected (highlighted by the colour purple).
equation (5).Our numerical framework is based on the method-of-lines recipe, where the partial differential equation is solved as a system of coupled ordinary differential equations (ODEs).We start by introducing the notation where the matrix evolution operator is defined as and finally reduce the system to the reduced ODE In this framework, the fields ϕ(t, x) and ψ(t, x) are evaluated on a discrete spatial grid x = {x i } N i=0 such that ϕ(t, x) → ϕ(t) and ψ(t, x) → ψ(t).The components ϕ(t, x i ) ≡ ϕ i (t) and ψ(t, x i ) ≡ ψ i (t) of the vectors ϕ(t) and ψ(t) are the fields evaluated on the grid points.Generically as we demonstrated in [18,19,20,21] we will be solving a system of coupled 2N+2 time-dependent ODEs given as where χı = χ(x ı ), ιı = ι(x ı ), Ṽı = Ṽ(x ı ), εı = ε(x ı ), ρı = ρ(x ı ) and I ıȷ is the N+1 identity matrix, obtained in Mathematica for example as, I = IdentityMatrix[N+1].For the Left: fractional energy error of the numerical solution up to when the pulse is inside the computational domain.High accuracy is observed and energy is conserved both exactly and numerically [18,19,20].Right: numerical error between the d'Alembert solution given by eq.( 3) and with the numerical recipe of eq.(19a -21b) in the numerical domain x ∈ [−7, 7] and t ∈ [0, 25.4742] .
wave equation in the (t, x) chart we have, χı = χ(x ı ) → 1 and ιı = Ṽı = εı = ρı = 0.For the hyperboloidal chart defining equation ( 5) we have For spatial discretisation in the spatial domain x ∈ [a, b] where a ≤ x i ≤ b we use Chebyshev collocation nodes ranging from 0 < i < N and follow the recipe given in Section 3.2.1 of [19].Particularly, in this work we use Mathematica Version 13.2 to obtain the n-th order differentiation matrices with the command, The time-component is resolved numerically by using the implicit-turned-explicit IMTEX Hermite time-integration schemes [18,19,20,21].In this work one opts for the fourth-order Hermite IMTEX scheme, given as For details and comparisons with higher-order IMTEX Hermite schemes, purely implicit IM Hermite and purely explicit EX Runge-Kutta schemes we refer the reader to [18,19,20].5) on a transformed hyperboloidal coordinate chart t → t(τ, ρ), x → x(ρ) numerically solved through the implementation of IMTEX Hermite time-integration scheme of order-4.As expected, outflow is observed, the initial pulse splits, travelling in opposite directions with equal field magnitudes until the solution leaves the numerical domain.However, unlike Figure 4, we now observe a settling down to a negative constant, at around ϕ(τ, ρ) ≈ −0.0832767, consistent with our theoretical predictions highlighted by equations ( 3)-( 13) (last plot on the right).

Numerical solution to the wave equation in the (t, x) coordinate chart with outflow boundary conditions
Having selected our numerical framework, and based on the discussion in previous chapters, we will first study the behaviour of our numerical solution in the (t, x) coordinate chart with outflow boundary conditions defined as, For initial data we pick, in accordance to equations (1b, 1c), We check the accuracy of our numerical solution by computing the energy, which should be conserved while the pulse remains inside the numerical computational domain  7. Initially the gaussian pulse has maximum amplitude (highlighted in red) it then splits with equal amplitudes as the simulation runs its course.However, unlike in Figure 5, the simulation settles down to a negative constant (highlighted by the colour purple) at around ϕ(τ, ρ) ≈ −0.0832767.[18,20].The energy is calculated as where here we use a Clenshaw-Curtis quadrature to integrate x ∈ [a, b].The fractional energy error is calculated as and the numerical error as, Our numerical domain spans x ∈ [−7, 7] and we use N=452 Chebyshev collocation nodes.Evolving our numerical solution we observe in Figure 4 outflow, as expected, for the entirety of the numerical simulation running time t ≈ [0., 25.4742].These results are further substantiated by the contour plot of Figure (5).By the end of the simulation i.e at around t = 25.4742all that is left, is noise from the numerical framework used.To gauge the numerical error associated with the implementation of the numerical scheme, we calculate the L-infinity norm, l ∞ of difference between the exact and numerical solutions as defined in eq.(24).As given by the plot on the right of Figure 6 the solution remains accurate throughout the entire simulation.The accuracy of our simulations is further verified by computing the fractional error of the energy as defined in eq.( 23).As we use IMTEX Hermite integration schemes we are able to conserve energy with near machine precision as showed in the left plot of Figure (6).Table 1.Numerical solution ϕ(τ, ρ) evaluated at relatively late times at around τ = 20, 000 verifying the negative permanent displacement predicted by eqs.(12)(13) for different S values.The numerical error given by the l ∞ is calculated as per eq.( 24) where the exact solution is given by eqs.(12)(13).
4.2.Numerical solution to the wave equation in the (t → t(τ, ρ), x → x(ρ)) hyperboloidal coordinate chart Our motivation is though to understand how the conversion of the Cartesian chart (t, x) to the hyperboloidal chart described by equation ( 5) i.e (t → t(τ, ρ), x → x(ρ)) affects our numerical results.To that effect we implement the same numerical recipe as described above, however, now there will be no need to implement boundary conditions, as equation ( 5) automatically enforces outflow.For initial data we make the following choice where here Ξ = 1/30.For the sake of consistency, in the computation of the numerical solution to ϕ in the hyperboloidal chart (τ, ρ), we choose N = 452 Chebyshev collocation nodes on a compact spatial domain spanning ρ ∈ [−S, S] where S = 1.Initially our simulation runs from τ ∈ [0, 12.1306] with an optimal time-step of ∆τ ≈ 0.0000242613.
As in the left plot of Figure 6, here, in the left plot of Figure 9, we too observe conservation of energy to a fractional error near machine precision.However, as demonstrated in Figure 7, even though the initial behaviour is the same as that observed in the first two plots (from left to right) of Figure 4, we now observe, Figure 7, a settling down to a negative constant, at around ϕ(τ, ρ) ≈ −0.0832767, reflecting the negative permanent displacement predicted by our theoretical results highlighted in equations ( 12)-( 13).This is further made clear by the contour plot given in Figure 8 describing the numerical evolution of ϕ(τ, ρ).To further assess the accuracy of our results we run our simulation for several values of S = 1, 4, 7, 10.Our interest is to verify that the value of the negative permanent displacement matches the theoretical prediction in eq.( 13), and hence it is necessary to significantly increase the simulation running time to τ = [0, 20, 000].As a compromise with speed, we also significantly inflate the time-step size of our simulation  5).Left: fractional energy error of the numerical solution up to when the pulse is inside the computational domain.High accuracy is observed and energy is conserved both exactly and numerically [18,19,20].Right: numerical error between the numerical and exact solution predicted by eqs.( 12)-( 13) for different values of S = 1, 4, 7, 10 at very late times τ = [0, 20, 000].As a compromise with simulation running times and unlike in Figure 6, or the left plot here depicted, we significantly increased our time-step size to ∆τ = 0.02.Our results are consistent with Figure 6.
to ∆τ = 0.02.As it is observed in Table 1, the simulation remained accurate for all values investigated provided one the initial gaussian pulse remained sufficiently compact, which can be easily achieved by calibrating the constant value Ξ.Both our theoretical and numerical results, as given by the plot on the right of Figure 9, are consistent with the accuracy of results obtained when solving the wave equation in the (t, x).
It is important to note, that this behaviour is not exclusive to the particular chart given in [8] but to any hyperboloidal chart.More generally, one would expect an analoguous effect in any coordinate transformation that modifies the time coordinate.For completion, we demonstrate this in Appendix A by implementing our numerical recipe to two different hyperboloidal slices.
Finally, the right plot in Figure 9 shows the same fractional energy error as that obtained for the numerical solution in the (t, x) chart, showing the merit of an hyperboloidal slice implementation.Not only does it further simplify our numerical implementation, it ensures numerical simulations remain accurate while automatically enforcing outflow behaviour.This advantage is of particular use when applied to problems such as those arising when modelling Extreme-Mass-Ratio-Inspirals (EMRIs) where accurate boundary conditions choices are less obvious and usually compromise the accuracy of the numerical solutions as discussed in [19,21].In Table 8 of [19] we give a brief outline of previous attempts within the community and it is noteworthy to highlight the merit of hyperboloidal approaches over traditional boundary conditions implementations on standard coordinate charts.

Conclusions
In this brief note we have obtained the hyperboloidal analogue of d'Alembert's solution to the wave equation in 1+1 dimensions.Our analysis has concentrated on a particular type of hyperboloidal foliation given in [8].The analogue expressions can be easily obtained for other families of hyperboloids such as those found in [22].Moreover, the results can also be adapted to incorporate the spherically 3+1 dimensional wave equation -for this one can make use of the method of spherical means to express the solution of the 3+1 equation in terms of a solution to 1+1 problems.A detailed discussion of this solution goes beyond the illustrative purposes of this article.
Our result brings to the fore the key role of explicit exact solutions in developing intuition into the nature of solutions to partial differential equations in non-standard physical and geometric settings.over σ ∈ [0, 1], a time-step size of ∆τ ≈ 0.00001897247 over a time interval of τ = [0, 18.974].As we can see from the first two plots, Top and Middle of Figure A1 a negative permanent displacement is observed at late times of ϕ(τ, σ) ≈ −0.0474046 and ϕ(τ, σ) ≈ −0.0672598 for slice A and slice B, respectively.In the latter, and unlike in the other two slices investigated in this work, the splitting of the gaussian pulse is not symmetric, this is due the fact the slicing was derived for curved spacetime and designed to improve the regularity of Schwarzschild (-like) cases [22].As with Figures 6,9 energy conservation is observed.‡ Checking the numerical error would require the derivation of an analogous hyperboloidal d'Alembert equation as showed in eqs.(12)(13) and is out of scope of this work.‡ We thank Rodrigo Panosso Macedo for having suggested and numerically independently verified our results with these slices [17].Bottom: fractional energy error of the numerical solutions of Slice A and Slice B up to when the pulse is inside the computational domain.High accuracy is observed and energy is conserved both exactly and numerically [18,19,20] and is consistent with the result observed in Figures 6,9.

Figure 1 .
Figure1.Top: contour plot of the d'Alembert solution corresponding to a choice of f with the shape of a bump (i.e. of compact support) and g = 0.This plot shows two pulses of half the size of the original bump moving in opposite directions.The support of the solution remains, at all times, compact.Bottom: contour plot of the d'Alembert solution for the choice f = 0 and g a bump.This choice of initial conditions generates a permanent displacement in the value of the scalar field in a compact domain.

Figure 4 .
Figure 4. Numerical evolution of the 1+1D wave equation defined by equation (1a) solved numerically through the implementation of IMTEX Hermite time-integration scheme of order-4 with open boundary conditions as defined in equations (20a)-(20b).As expected outflow is observed, the initial pulse splits, travelling in opposite directions with equal field magnitudes until the solution fully leaves the numerical domain (last plot on the right).

Figure 5 .
Figure 5. Contour plot further corroborating observations in Fig.(4).Initially the gaussian pulse has maximum amplitude (highlighted in red) it then splits with equal amplitudes as the simulation runs its course.Outflow is observed, and by the end of the simulation only numerical residue is detected as expected (highlighted by the colour purple).

Figure 6 .
Figure 6.Numerical solution to the wave equation in the (t, x) coordinate chart.Left: fractional energy error of the numerical solution up to when the pulse is inside the computational domain.High accuracy is observed and energy is conserved both exactly and numerically[18,19,20].Right: numerical error between the d'Alembert solution given by eq.(3) and with the numerical recipe of eq.(19a -21b) in the numerical domain x ∈ [−7, 7] and t ∈ [0,25.4742]

Figure 7 .
Figure 7. Numerical evolution of the 1+1D wave equation defined by equation (5) on a transformed hyperboloidal coordinate chart t → t(τ, ρ), x → x(ρ) numerically solved through the implementation of IMTEX Hermite time-integration scheme of order-4.As expected, outflow is observed, the initial pulse splits, travelling in opposite directions with equal field magnitudes until the solution leaves the numerical domain.However, unlike Figure4, we now observe a settling down to a negative constant, at around ϕ(τ, ρ) ≈ −0.0832767, consistent with our theoretical predictions highlighted by equations (3)-(13) (last plot on the right).

Figure 8 .
Figure 8. Contour plot further corroborating observations in Figure7.Initially the gaussian pulse has maximum amplitude (highlighted in red) it then splits with equal amplitudes as the simulation runs its course.However, unlike in Figure5, the simulation settles down to a negative constant (highlighted by the colour purple) at around ϕ(τ, ρ) ≈ −0.0832767.

Figure 9 .
Figure 9. Numerical solution to the wave equation in the (τ, ρ) coordinate chart described by equation (5).Left: fractional energy error of the numerical solution up to when the pulse is inside the computational domain.High accuracy is observed and energy is conserved both exactly and numerically[18,19,20].Right: numerical error between the numerical and exact solution predicted by eqs.(12)-(13) for different values of S = 1, 4, 7, 10 at very late times τ = [0, 20, 000].As a compromise with simulation running times and unlike in Figure6, or the left plot here depicted, we significantly increased our time-step size to ∆τ = 0.02.Our results are consistent with Figure6.

Figure A1 .
Figure A1.Numerical solution to the wave equation in two hyperboloidal (τ, σ) coordinate charts.Top, Middle: contour plots showing the negative permanent displacement of two different slices, slice A and slice B of a value of ϕ(τ, ≈ −0.0474046 and ϕ(τ, σ) ≈ −0.0672598, respectively.Bottom: fractional energy error of the numerical solutions of Slice A and Slice B up to when the pulse is inside the computational domain.High accuracy is observed and energy is conserved both exactly and numerically[18,19,20] and is consistent with the result observed in Figures6,9.