We explicitly construct the analogue of the dAlembert solution to the 11 wave equation in an hyperboloidal setting. This hyperboloidal dAlembert 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.
['Wave equation', 'Exact solutions', 'Hyperboloids', "D'Alembert solution", 'General Relativity and Quantum Cosmology', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Mathematical Physics', 'Physics - Computational Physics']
The dAlembert solution in hyperboloidal foliations
{'Juan A. Valiente Kroon & Lidia Gomes da Silva': "School of Mathematical Sciences, Queen Mary, University of London, Mile End Road, London E1 4NS, UK  \nE-mail: j.a.valiente-kroon@qmul.ac.uk, lidiajoana@pm.me  \n10 September 2024  \nAbstract. 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.", '1. 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].  \nAn 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].  \nIn 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.  \nKey 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.  \nThe 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.  \nSolutions 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.  \nFigure 1. 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.  \n<!-- image -->", '2. 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.', "2.1. 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  \n∂ 2 t ϕ -c 2 ∂ 2 x ϕ = 0 , t ≥ 0 , x ∈ R , (1a)  \nϕ (0 , x ) = f ( x ) , (1b)  \n∂ t ϕ (0 , x ) = g ( x ) , (1c)  \nwhere 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.  \nIt is well-known that the general solution to the 1 + 1-dimensional wave equation is of the form  \nϕ ( t, x ) = F ( x -t ) + G ( x + t ) , (2)  \nwhere 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 :  \nϕ ( t, x ) = 1 2 ( f ( x -t ) + f ( x + t ) ) + 1 2 ∫ x + t x -t g ( s )d s. (3)  \nIntuition 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.  \nAn 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  \ng ( x ) = 1 1 + x 2  \nwhich can be easily integrated. In this case d'Alembert's formula gives the solution  \nϕ ( t, x ) = 1 2 ( arctan( x + t ) -arctan( x -t ) ) .  \nObserve 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.  \nRemark. In view of the linearity of the wave equation a generic solution is a combination of the two effects discussed previously.", '3. The wave equation in the hyperboloidal setting': "Following [8] we introduce coordinates τ and ρ via  \nτ ≡ t -h ( x ) , x ≡ ρ Ω ,  \nwhere the height function h is given by  \nh ( x ) ≡ √ S 2 + x 2 (4)  \nFigure 3:  \nFigure 5: Ω( ρ ) = 1 2 ( 1 -ρ 2 S 2 ) .  \n<!-- image -->  \ni  \n-  \nFigure 4: t Figure 2. Examples of hyperboloids in the Minkowski spacetime corresponding, from top to bottom, to the choices of S = 7, S = 1, S = 1 2 , S = 1 4 and S = 1 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.  \nwhere S is a positive number and  \nx  \n3 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.  \nThe endpoints of the hyperboloid correspond to the sets for which ρ = ± S . In particular, the 1 + 1 wave equation takes the form  \n∂ 2 τ ϕ + 2 ρ S ∂ τ ∂ ρ ϕ -Ω 2 ∂ 2 ρ ϕ + 2 S Ω S 2 + ρ 2 ∂ τ ϕ + (3 S 2 + ρ 2 ) ρ Ω S 2 ( S 2 + ρ 2 ) ∂ ρ ϕ = 0 , (5)  \nwhere now ϕ = ϕ ( τ, ρ ). The above equation evaluated at ρ = ± S takes the form  \n∂ τ ( ∂ τ ∓ 2 ∂ ρ ) ϕ = 0 . (6)  \nFrom 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.  \ni  \n-  \ni  \n0  \nI  \n-  \nI  \n-  \nFigure 6:  \n<!-- image -->  \nFigure 7:  \nFigure 5: 3 4 Figure 3. The hyperbolidal foliation corresponding to the choice S = 1 for the height function h given by equation (4). Left: diagram in Cartesian coordinates. Right: Penrose diagram. The hyperboloids correspond, from bottom to top, to the times τ = { 0 , 1 , 2 , 3 } . The 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.  \nWe will consider the initial value problem for equation (5) on the hyperboloid with τ = 0. For this we prescribe, as usual  \nϕ (0 , ρ ) = f ( ρ ) , (7a)  \n∂ τ ϕ (0 , ρ ) = g ( ρ ) , (7b)  \nwhere 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).  \nRemark 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.", '3.1. 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  \nx + t = ζ + τ, x -t = χ -τ,  \nwhere we have set  \nζ ≡ S ( ρ + S ) S -ρ , χ ≡ S ( ρ -S ) S + ρ , (8)  \nso that the general solution to the wave equation (5) is given by  \nϕ ( τ, ρ ) = F ( ζ + τ ) + G ( χ -τ ) . (9)  \n3  \ni  \n0  \nObserve that ζ →∞ as ρ → S and χ →∞ as ρ →-S while ζ ( -S ) = χ ( S ) = 0.  \nA calculation then readily shows that  \nϕ (0 , ρ ) = F ( S ( ρ + S ) S -ρ ) + G ( S ( ρ -S ) S + ρ ) , (10a)  \n∂ τ ϕ (0 , ρ ) = F ' ( S ( ρ + S ) S -ρ ) -G ' ( S ( ρ -S ) S + ρ ) , (10b)  \nwhere ' 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  \nF ' ( S ( ρ + S ) S -ρ ) -G ' ( S ( ρ -S ) S + ρ ) = g ( ρ ) , 2 S 2 ( S -ρ ) 2 F ' ( S ( ρ + S ) S -ρ ) + 2 S 2 ( S + ρ ) 2 G ' ( S ( ρ -S ) S + ρ ) = f ' ( ρ ) .  \nThe above is a linear system for  \nF ' ( S ( ρ + S ) S -ρ ) , G ' ( S ( ρ -S ) S + ρ ) .  \nIts solution can be readily be found to be given by  \nF ' ( S ( ρ + S ) S -ρ ) = ( S 2 -ρ 2 ) 2 4 S 2 ( S 2 + ρ 2 ) f ' ( ρ ) + ( S -ρ ) 2 2( S 2 + ρ 2 ) g ( ρ ) , (11a)  \nG ' ( S ( ρ -S ) S + ρ ) = ( S 2 -ρ 2 ) 2 4 S 2 ( S 2 + ρ 2 ) f ' ( ρ ) -( S -ρ ) 2 2( S 2 + ρ 2 ) g ( ρ ) . (11b)  \nNow, re-expressing (11a) in terms of ζ as given by (8) one has that  \nF ' ( ζ ) = S 2 S 2 + ζ 2 g ( S ( ζ -S ) S + ζ ) + 2 S 2 ζ 2 ( S + ζ ) 2 ( S 2 + ζ 2 ) f ' ( S ( ζ -S ) S + ζ ) .  \nIntegration of this equation yields  \nF ( ζ ) -F (0) = ∫ ζ 0 S 2 S 2 + s 2 g ( S ( s -S ) S + s ) d s + ∫ ζ 0 2 Ss 2 ( S + s ) 2 ( S 2 + s 2 ) f ' ( S ( s -S ) S + s ) d s.  \nUsing integration by parts one can remove the derivative in the second integral so that  \nF ( ζ ) = F (0) + ∫ ζ 0 S 2 S 2 + s 2 g ( S ( s -S ) S + s ) d s + ζ 2 S 2 + ζ 2 f ( S ( ζ -S ) S + ζ ) -2 S 2 ∫ ζ 0 s ( S 2 + s 2 ) 2 f ( S ( s -S ) S + s ) d s.  \nNow, using (9) one finds, after some manipulations, that  \nG ( χ ) = χ 2 S 2 + χ 2 f ( S ( S + χ ) S -χ ) -F (0) -∫ -S 2 /χ 0 S 2 S 2 + s 2 g ( S ( s -S ) S + s ) d s +2 S 2 ∫ -S 2 /χ 0 s ( S 2 + s 2 ) 2 f ( S ( s -S ) S + s ) d s.  \nCombining the above expressions using the replacements  \nζ ↦→ ζ + τ, χ ↦→ χ -τ,  \none obtains the hyperboloidal d'Alembert formula  \nϕ ( τ, ρ ) = ( ζ + τ ) 2 S 2 +( ζ + τ ) 2 f ( S ( ζ + τ -S ) S + ζ + τ ) + ( χ -τ ) 2 S 2 +( χ -τ ) 2 f ( S ( S + χ -τ ) S -χ + τ ) -2 S 2 ∫ ζ + τ -S 2 / ( χ -τ ) s ( S 2 + s 2 ) 2 f ( S ( s -S ) s + S ) d s + ∫ ζ + τ -S 2 / ( χ -τ ) S 2 S 2 + s 2 g ( S ( s -S ) s + S ) d s. (12)  \nThis 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.  \nRemark 2. It can be readily be verified using the relation  \nζ = -S 2 χ  \nthat indeed ϕ (0 , ρ ) = f ( ρ ) and ∂ τ ϕ (0 , ρ ) = g ( ρ ). 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).  \nRemark 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.", '3.2. 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  \nI ( τ, ρ ) ≡ -2 S 2 ∫ ζ + τ -S 2 / ( χ -τ ) s ( S 2 + s 2 ) 2 f ( S ( s -S ) s + S ) d s,  \nwith 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.  \nThe 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  \nlim τ →∞ I ( τ, ρ ) = -2 S 2 ∫ ∞ 0 s ( S 2 + s 2 ) 2 f ( S ( s -S ) s + S ) d s < 0 . (13)  \nFigure 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).  \n<!-- image -->  \nThe 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  \nϕ ( ∞ ) ≡ lim τ →∞ I ( τ, ρ ) < 0 ,  \nas  \nlim τ →∞ ( ( ζ + τ ) 2 S 2 +( ζ + τ ) 2 f ( S ( ζ + τ -S ) S + ζ + τ ) + ( χ -τ ) 2 S 2 +( χ -τ ) 2 f ( S ( S + χ -τ ) S -χ + τ )) = f ( S ) + f ( -S ) ,  \nwhich 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.", '4. 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  \nFigure 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).  \n<!-- image -->  \nequation (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  \n∂ t U = L U , U = ( ϕ ψ ) , (14)  \nwhere the matrix evolution operator is defined as  \nL = ( 0 1 L 1 L 2 ) , (15)  \nand finally reduce the system to the reduced ODE  \nd U dt = LU . (16)  \nIn 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  \nd dt ( ϕ ı ψ ı ) = N ∑ ȷ =0 ( 0 I ıȷ ˜ χ ı D (2) ıȷ +˜ ι ı D (1) ıȷ + ˜ V ı I ıȷ ˜ ε ı D (1) ıȷ + ˜ ϱ ı I ıȷ )( ϕ ȷ ψ ȷ ) , ı = 0 , 1 , ..., N (17)  \nwhere ˜ χ ı = ˜ χ ( x ı ) , ˜ ι ı = ˜ ι ( x ı ) , ˜ V ı = ˜ V ( x ı ) , ˜ ε ı = ˜ ε ( x ı ) , ˜ ϱ ı = ˜ ϱ ( x ı ) and I ıȷ is the N+1 identity matrix, obtained in Mathematica for example as, I = IdentityMatrix[N+1] . For the  \n<!-- image -->  \nFigure 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]  \n<!-- image -->  \n.  \nwave equation in the ( t, x ) chart we have, ˜ χ ı = ˜ χ ( x ı ) → 1 and ˜ ι ı = ˜ V ı = ˜ ε ı = ˜ ϱ ı = 0. For the hyperboloidal chart defining equation (5) we have  \n˜ ε ( ρ ) = 2 ρ S , (18a)  \n˜ ϱ ( ρ ) = 2 S Ω S 2 + ρ 2 , (18b)  \n˜ χ ( ρ ) = Ω 2 , (18c)  \n˜ V ( ρ ) = 0 , (18d)  \n˜ ι ( ρ ) = (3 S 2 + ρ 2 ) ρ Ω S 2 ( S 2 + ρ 2 ) . (18e)  \nFor 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,", 'Dn = NDSolve\'FiniteDifferenceDerivative[Derivative[i], x, DifferenceOrder -> "Pseudospectral"]@"DifferentiationMatrix"].': '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  \nU n +1 = U n +(∆ t L ) · HFH4 · U n , (19a)  \nHFH4 = ( I -∆ t L 2 · ( I -∆ t L 6 )) -1 . (19b)  \nFor 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].  \nFigure 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 order4. 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).  \n<!-- image -->', '4.1. 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,  \n∂ t ϕ ( t, x 0 ) = c ∂ x ϕ ( t, x ) , (20a)  \n∂ t ϕ ( t, x N ) = -c ∂ x ϕ ( t, x ) . (20b)  \nFor initial data we pick, in accordance to equations (1b, 1c),  \nf ( x ) = e -16 x 2 , (21a)  \ng ( x ) = 0 . (21b)  \nWe check the accuracy of our numerical solution by computing the energy, which should be conserved while the pulse remains inside the numerical computational domain  \nFigure 8. Contour plot further corroborating observations in Figure 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.  \n<!-- image -->', '[18, 20]. The energy is calculated as': 'E = 1 2 ∫ b a ( ∥ ϕ ( t, x ) ∥ 2 + ∥ ψ ( t, x ) ∥ 2 ) dx, (22)  \nwhere here we use a Clenshaw-Curtis quadrature to integrate x ∈ [ a, b ]. The fractional energy error is calculated as  \nδE E 0 = E ( t f ) -E t 0 E t 0 , (23)  \nand the numerical error as,  \n|| ∆ ϕ ( t, x ) || = ϕ exact ( t, x ) -ϕ numerical ( t, x ) . (24)  \nOur 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 . 4742 all 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).  \nTable 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  \nϕ (0 , ρ ) = exp ( -ρ 2 2 Ξ 2 ) , (25a)  \n∂ τ ϕ (0 , ρ ) = 0 , (25b)  \nwhere 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.  \nAs 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  \n<!-- image -->  \nFigure 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 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.  \n<!-- image -->  \nto ∆ τ = 0 . 02. As it is observed in Table 1, the simulation remained accurate for all values investigated provided one ensured 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 ).  \nIt 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.  \nFinally, 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.', '5. 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 the 1+1 problems. A detailed discussion of this solution goes beyond the illustrative purposes of this article.  \nOur 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.", 'Acknowledgements': "We would like to thank Rodrigo Panosso Macedo and Nelson Eir'o for independently checking our numerical results in 2020-2021. We would also like to thank Anil Zenginoglu, Peter Diener, Alex Va˜no-Vi˜nuales, Enno Harms and Sebastiano Bernuzzi for discussions and previous work. JAVK acknowledges the financial support of EPSRC grant 'Geometric scattering methods for the conformal Einstein field equations', EP/X012417/1.", 'Appendix A. Negative permanent displacement detection in other hyperboloidal coordinate charts numerical solutions to the wave equation': "For completion we now implement two different hyperboloidal slices as studied by [22] see also [23, 24, 25, 26], defined by the coordinate chart transformation, ¯ t → ( ¯ t ( τ, σ ) , ¯ x → ¯ x ( σ )). Following the notation of [22] we define,  \n¯ t = τ -h ( σ ) , (A.1)  \n¯ x = g ( σ ) , (A.2)  \nwhere generically as shown in [22],  \n¯ t = t λ , ¯ x = r ∗ λ , (A.3)  \nwith λ to be specified further for each hyperboloidal slice we will study. r ∗ is the tortoise coordinate given as dr/dr ⋆ = f ( r ) where f ( r ) = (1 -2 M/r ), in an extended domain of r ⋆ = [ ∞ , ∞ ]. Furthermore, h ( σ ) is the height function implementing the hyperboloidal slicing and g ( σ ) is the spatial compactification bringing ¯ x ∈ [ -∞ , ∞ ] down to a compact domain of σ = [ a, b ].  \nThe first slice we will study, which we will here call Slice A uses compactified hyperboloids given by Biz'on-Mach coordinates and maps R to ] -1 , 1[ via,  \n¯ t = τ -h ( σ ) = τ -1 2 ln (1 -σ 2 ) , (A.4)  \n¯ x = arctanh( σ ) , (A.5)  \nwith λ = 1 and ¯ x ∈ [ -∞ , ∞ ] being compactified to the domain of [ a, b ] where a = -1 , b = 1. The evolution operators L 1 , L 2 as defined in equations (15-17) now transform to,  \n˜ ε ( σ ) = -2 σ, (A.6a)  \n˜ ϱ ( σ ) = -1 , (A.6b)  \n˜ χ ( σ ) = (1 -σ 2 ) , (A.6c)  \n˜ V ( σ ) = 0 , (A.6d)  \n˜ ι ( σ ) = -2 σ. (A.6e)  \nFor initial data, we inject the following gaussian pulse  \nϕ (0 , σ ) = exp ( -a ( -1 2 log (1 -σ 2 ) -arctan ( σ ) ) 2 ) , (A.7)  \n∂ τ ϕ (0 , σ ) = 0 , (A.8)  \nwhere a = 350. The second slice, hereafter slice B, commonly referred to as minimal gauge is attained via the transformations,  \n¯ t = τ -h ( σ ) = τ -1 2 ( ln σ +ln(1 -σ ) -1 σ ) , (A.9)  \n¯ x = 1 2 ( ln (1 -σ ) + 1 σ -ln σ ) , (A.10)  \nwhere here λ = 4 M and ¯ x = r ∗ /λ , and thus ¯ x ∈ [ -∞ , ∞ ] is mapped to the compactified domain of [ a, b ] where a = 0 , b = 1.  \nFor the second slice, slice B, commonly referred to as the minimal gauge , the evolution operators are exactly as defined by equations (27)-(31) of [21] with the difference that in this work ˜ V ( σ ) = 0. For initial data for slice B, we pick,  \nϕ (0 , σ ) = exp ( -( σ -0 . 5) 2 2 Ξ ) , (A.11a)  \n∂ τ ϕ (0 , σ ) = 0 , (A.11b)  \nwhere Ξ = 1 / 3000. For the numerical simulation of slice A, we use N = 256 Chebyshev collocation nodes on spatial domain spanning over σ ∈ [ -1 , 1], a time-step size of ∆ τ = 0 . 0000379449 over a time interval of τ = [0 , 18 . 974]. In the numerical simulation of slice B, we use N = 256 Chebyshev collocation nodes on spatial domain spanning  \nover σ ∈ [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.  \nFigure 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 Figures 6,9.  \n<!-- image -->", 'References': "- [1] J. A. Valiente Kroon. Conformal Methods in General Relativity . Cambridge University Press, 2016.\n- [2] A. Zenginoglu. Hyperbolic times in minkowski space. In arXiv:2404.01528 , 2024.\n- [3] H. Friedrich. Cauchy problems for the conformal vacuum field equations in general relativity. Comm. Math. Phys. , 91:445, 1983.\n- [4] H. Friedrich. On the existence of n-geodesically complete or future complete solutions of Einstein's field equations with smooth asymptotic structure. Comm. Math. Phys. , 107:587, 1986.\n- [5] A. Zenginoglu. A geometric framework for black hole perturbations. Phys. Rev. D , 83:127502, 2011.\n- [6] A. Zenginoglu and L. E. Kidder. Hyperboloidal evolution of test fields in three spatial dimensions. Phys. Rev. D , 81:124010, 2010.\n- [7] S. Gautam, A. Va no Vi˜nuales, D. Hilditch, and S. Bose. Summation by parts and truncation error matching on hyperboloidal slices. Phys. Rev. D , 103:084045, 2021.\n- [8] A. Zenginoglu. Hyperboloidal layers for hyperbolic equations on unbounded domains. J. Comput. Phys. , 230:2286, 2011.\n- [9] T. Hagstrom and S. Lau. Radiation boundary conditions for maxwell's equations' a review of accurate time-domain formulations. J. Comp. Math. , 25:305, 2007.\n- [10] S. Lau. Radiation boundary kernels for time-domain wave propagation on black holes: thoery and numerical mathods. J. Comp. Phys. , 199:376, 2004.\n- [11] S. Lau. Analytic structure of radiation boundary kernels for blackhole perturbations. J. Math. Phys. , 46:102503, 2005.\n- [12] S. Lau. Fast evaluation of far-field signals for time-domain wave propagation. J. Sci. Comp. , 64:647, 2015.\n- [13] M. J. Grote and J. B. Keller. On nonreflecting boundary conditions. J. Comp. Phys. , 122:231, 1995.\n- [14] L. C. Evans. Partial Differential Equations . American Mathematical Society, 1998.\n- [15] F. John. Partial differential equations . Springer, 1991.\n- [16] Enno Harms. Numerical solution of the 2+1 teukolsky equation on a hyperboloidal foliation of the kerr spacetime. Master's thesis, Friedrich-Schiller-Universitat Jena, Physikalisch-Astronomische Fakultat, Deutschland, 2012.\n- [17] Rodrigo Panosso Macedo. Fully-spectral minimal gauge time-domain code. Private Communication , 2021.\n- [18] L. J. Gomes Da Silva. Numerical Algorithms for the modelling of EMRIs in the time domain . PhD dissertation, QMUL University, School of Mathematical Sciences, 2023.\n- [19] Lidia J Gomes Da Silva. Discotex: Discontinuous collocation and implicit-turned-explicit (imtex) integration symplectic, symmetric numerical algorithms with higher order jumps for differential equations with numerical black hole perturbation theory applications. arXiv preprint arXiv:2401.08758 , 2024.\n- [20] Michael F O'Boyle, Charalampos Markakis, Lidia J Gomes Da Silva, Nelson Eir'o, Rodrigo Panosso Macedo, and Juan A Valiente Kroon. Conservative evolution of black hole perturbations with time-symmetric numerical methods. arXiv preprint arXiv:2210.02550 , 2022.\n- [21] Lidia J Gomes Da Silva, Rodrigo Panosso Macedo, Jonathan E Thompson, Juan A Valiente Kroon, Leanne Durkan, and Oliver Long. Hyperboloidal discontinuous time-symmetric numerical algorithm with higher order jumps for gravitational self-force computations in the time domain. arXiv preprint arXiv:2306.13153 , 2023.\n- [22] Jos'e Luis Jaramillo, Rodrigo Panosso Macedo, and Lamis Al Sheikh. Pseudospectrum and black hole quasinormal mode instability. Physical Review X , 11(3):031003, 2021.\n- [23] M. Ansorg and R. P. Macedo. Axisymmetric fully spectral code for hyperbolic equations. J. Comp. Phys. , 276:357, 2014.  \n- [24] M. Ansorg and R. P. Macedo. Spectral decomposition of black-hole perturbations on hyperboloidal slices. Phys. Rev. D , 93:124016, 2016.\n- [25] R. P. Macedo, J. L. Jaramillo, and M. Ansorg. Axisymmetric fully spectral code for hyperbolic equations. Phys. Rev. D , 98:124005, 2018.\n- [26] L. Al Sheik. Scattering resonances and Pseudospectrum : stability and completeness aspects in optical and gravitational systems . PhD thesis, Institut de Math'ematiques de Bourgogne, Dijon, France, 2022."}