{'fraction_non_alphanumeric': 0.04139389666633095, 'fraction_numerical': 0.01597053364602967, 'mean_word_length': 4.4308485700890765, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A spectral method algorithm for numerical simulations of gravitational fields

A numerical study of the Einstein field equations, based on the 3+1 foliation of the spacetime, is presented. A pseudo-spectral technique has been employed for simulations in vacuum, within two different formalisms, namely the Arnowitt-Deser-Misner (ADM) and the conformal Baumgarte-Shapiro-Shibata-Nakamura (BSSN) approach. The numerical code is based on the Fourier decomposition, accompanied by different filtering techniques. The role of the dealiasing, as well as the influence of the filter type, has been investigated. The algorithms have been stabilized via a novel procedure that controls self-consistently the regularity of the solutions. The accuracy of the model has been validated through standard testbeds, revealing that the filtered pseudo-spectral technique is among the most accurate approaches. Finally, the procedure has been stressed via black hole dynamics and a new strategy, based on hyperviscous dissipation that suppresses spurious boundary problems, has been proposed. The model represents a valid tool of investigation, particularly suitable for the inspection of small scale nonlinear phenomena in gravitational dynamics.


Introduction
In past decades, great effort has been devoted to the numerical study of gravitation. More recently, numerical simulations of the Einstein field equations have been improved significantly. Many of these studies have successfully described the inspiraling of compact objects in their final coalescence phase (and the relative emission of gravitational waves), becoming therefore an extraordinary guidance for the observational campaign of LIGO and VIRGO [1][2][3][4][5][6][7][8][9][10]. Several works have been dedicated to the dynamics of Kerr-type black holes, as well as the understanding of the plasma dynamics in their accretion discs. These simulations helped to interpret novel observations, including the celebrated M87 Event Horizon Telescope image of a supermassive black hole [11][12][13][14][15][16][17].
Because of the intrinsic nonlinearity and complexity of the governing equations, several works have been dedicated to the study of both the stability and accuracy of the numerical solutions [18][19][20][21]. Different forms of the governing equations have been proposed and different numerical strategies have been developed, in a variety of initial conditions known as gravitational testbeds [22]. Many factors affect the dynamics and the accuracy of the solutions, including the formulation, the choice of the gauge conditions, as well as the boundary conditions, the accuracy of the space-time integration, the use of external dissipation and so on [23,24].
In this paper, we adopt two different 3 + 1 formalisms, based on the standard Arnowitt-Deser-Misner (ADM) decomposition [25,26] and the conformal Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation [27][28][29][30][31]. We propose a very accurate numerical strategy, in vacuum conditions, based on classical spectral methods. In analogy with high-precision simulations of fluid and plasma dynamics, our numerical code is based upon a Fourier pseudo-spectral method, in concert with filtering techniques. These anti-aliasing techniques are important in direct numerical simulations, since they stabilize the numerical evolution, establishing a direct connection between pseudo-spectral methods and exact Galerkin representations [32]. The main advantage of this approach, apart from its simplicity, is the accuracy of the solutions due to the spectral projection, while its limitation relies on the imposition of periodicity at the boundaries. Nevertheless, similar spectral decomposition can be locally used in the flux-reconstruction of finite volume methods, when the 3+1 Einstein field equations are represented in a quasi-hyperbolic form [33,34].
Spectral methods are among the most accurate and performing techniques for numerical relativity. Among all, the SpEC code represents a first, clear example of highly performing spectral algorithm applied to gravitation [35,36], where the basic model relies on a first-order, symmetric-hyperbolic, generalized harmonic formulation. Alternatively, in the framework of the BSSN representation, another good example, based on spectral representation of the spherical harmonics, is given by the SGRID code [37]. Both of the above are multi-domain, pseudo-spectral codes particularly suitable for the study of binary systems. These important, inspiring studies suggest that the spectral decomposition might be both accurate and stable. On this path we proceed here, concentrating on an alternative pseudo-spectral technique and proposing new strategies that prevent numerical pathologies, typical of spectral methods simulations.
Apart of the advantage given by the accuracy of the spectral decomposition, the periodic representation can be a valid (and easy) choice for the study of small scale gravitational dynamics, especially in cases where the fields are strongly nonlinear and form local homogeneous patterns, as it can be in the vicinity of the merging-region between compact objects, as well as in the coalescence/collapse of multiple massive bodies [38][39][40]. Note that, in the case of chaotic patterns with a finite characteristic scale [41], periodic conditions can be a valid strategy for the study of small scale dynamics, as commonly done in simulations of turbulent astrophysical plasmas [42].
In the present work, the accuracy of our numerical model has been validated through a series of standard tests, like the gauge wave, the robust stability test, the linear gravitational waves and the Gowdy spacetime, as suggested in previous pioneering works [22,[43][44][45]. The code successfully passed all the numerical tests. Apart from the classical tests, we also proposed a new class of initial data, that satisfies the constraint equations and leads to the formation of standing waves of the metric. In order to minimize the aliasing instability due to nonlinear terms, we propose different antialiasing filters, showing how smoother filters might lead to more stable solutions. The algorithms have been stabilized via an adaptive time-refinement, based on a method that self-consistently controls the regularity of the solutions, at each time-step. Finally, we stress the code via black hole dynamics. In order to suppress spurious boundary problems, we propose a new strategy based on hyper-viscous dissipation. The work is organized as follows. In section 2 we discuss two different 3+1 models of general relativity, namely the ADM and the BSSN, while in section 3 we explain our numerical technique. Here we present the decomposition, the spacetime filters and the time-step control. The classical numerical testbeds are presented in section 4, while in section 5 we introduce a novel technique that suppresses boundary effects, tested via gravitational wavepackets and a head-on collision of black holes. Finally, in section 6, we summarize the main achievements of our work and discuss our conclusions.

Spacetime foliation in general relativity
We use the same spectral algorithm for two different formulations, namely the ADM and the BSSN set of equations. In particular, we build an algorithm based on dealiased spectral methods and compare the above approaches through all the classical testbeds. As follows, we briefly summarize these two models of general relativity.

The ADM formulation
Gravitational fields obey to the celebrated Einstein field equation, which in vacuum can be written as be satisfied and these are known as the Hamiltonian and the momentum constraints, namely: In the above relations K ij is the extrinsic curvature and K = γ ij K ij is its trace, while R is the scalar curvature and D i is the covariant derivative. Any initial condition must obey equations (2)-(3) and this is known as the initial data problem [47][48][49]. Associated to the above constraints, the evolution equations for the spatial metric and the extrinsic curvature are given by where L β is the Lie derivative with respect to the shift vector β k . Finally, we need the evolution equation for the lapse and the shift, choosing therefore the slicing conditions. In the framework of the Bona-Massó formalism [50,51], we choose with f (α) > 0 (but arbitrary) [52,53]. In the case in which f (α) = 1, the above equation describes the harmonic slicing for the lapse [54]. In the ADM code we will set β i = 0 therefore we do not need an equation for the shift. We will introduce such an equation for β i , in the next subsection, within the BSSN formalism.

The BSSN formulation
Several formulations alternative to the ADM have been proposed in the literature [55,56]. A popular extension of the ADM is based on the work by Baumgarte and Shapiro [27] and by Shibata and Nakamura [28]. A model known as the BSSN approach. In order to introduce the BSSN formulation, let consider first a conformal rescaling of the spatial metric of the form where ψ is a conformal factor that can in principle be chosen in a number of different ways. Here we follow the approach suggested by Campanelli et al. [57], where the conformal factor is written as ψ −4 = χ, choosing χ such that the conformal metric γ ij has unit determinant, namely χ = γ −1/3 , where γ is the determinant of the physical spatial metric. In the BSSN, the extrinsic curvature is subjected to the same conformal transformation described by equation (7) and its trace K is evolved as an independent variable. Instead of K ij , the new variables are K = γ ij K ij and A ij = χ K ij − 1 3 γ ij K -a choice that implies tr A ij = 0. Finally, a new field is usually introduced, namely the Christoffel symbols of the conformal metric Γ i = γ jk Γ i jk . This is a crucial point since it has been proven that this choice further stabilizes numerical algorithms [22]. By using the above definitions and change of variables, the system of BSSN equations reads: In the above set, "T F " indicates the trace-free parts of a tensor object and D i is the covariant derivative associated with the physical three-metric γ ij . The Ricci tensor associated with the physical metric is separated into two contributions, is the Ricci tensor related to the conformal metric and denotes the additional term that depends on χ. In addition to equation (6), we use the following slicing conditions [58]: Here the factor 3/4 is somewhat arbitrary, but leads to good numerical results and the parameter η is a positive constant [59]. In all the evolution equations above, Γ i is replaced by −∂ j γ ij wherever it is not differentiated. In addition to the evolution equations, the BSSN variables must also obey to the constraints in equations (2)-(3), that now become The specific choice of evolution variables introduces five additional constraints, namely: During the numerical simulations, we enforce the algebraic constraints in equations (17) and (18), as discussed in [20]. The remaining constraints, H, M i and G i are not actively enforced and are used to monitor the accuracy of the numerical solutions. For more details on these additional constraints, see for example Ref. [22]. As a summary, we solve equations (8)- (12) for the new variables, by decomposing the Ricci tensor with (13)- (14) and closing the system for the lapse α and the shift β i via the slicing equations (6), (15) and (16). Finally, it is worth mentioning also the generalized harmonic formulation, for example used by spectral or other efficient finite difference codes [60]. This latter formulation is shown to be symmetric hyperbolic and admits well posed numerical schemes, being in fact the first method to achieve stable evolutions of more than a single orbit, as discussed in the work by Pretorious [61].

The numerical technique
In order to solve the gravitational fields equations, using both the ADM and the BSSN formalism, we make use of spectral techniques [62][63][64][65][66]. More specifically, in our code we use the pseudo-spectral method, which is a technique characterized by a recurrent switch, back-and-forth, between the physical (real) and the Fourier (complex) space. The rule is quite simple: derivatives are computed in the spectral space while products are computed in the physical space (in order to avoid convolutions.) This iterative transformation, as we will see in the next subsection, needs to be regularized via a dealiasing procedure. Classical spectral and pseudo-spectral methods are widely used for complex problems that involve strong couplings and nonlinearities. The main advantage of these techniques relies on the fact that each field and all its derivatives are represented at the collocation points in a very accurate way. This leads to very precise solutions, with an error on the computation of derivatives that is on the order of the machine truncation error.
In order to compute any spatial derivative, imposing periodic boundary conditions in all the three spatial dimensions, our code makes use of standard Fast Fourier Transforms (FFTs). More specifically, each field f (x, t) is projected on a 3D spatial box, over a lattice with a number of meshes N = (N x , N y , N z ). We decompose the ADM and the BSSN fields as: with f N being the fields of the 3+1 foliation (i.e. γ ij , K, K ij and so on.), at the grid points.
As usual, f k (t) ∈ C are the Fourier coefficients and k = (k x , k y , k z ) is the k-vector. In each Cartesian direction, k = 2π L 0 m, L 0 is the box length, m = 0, ±1, ±2, . . . ± N k with N k = N/2 being the Nyquist mode and N is the number of points along this direction.
Even if pseudo-spectral methods are among the most precise numerical techniques, particularly suitable for solving small scale complex problems such as wave dynamics, turbulence and so on, great care must be dedicated to the computation of nonlinear terms. Nonlinearity characterizes gravitation, as in the case of hydrodynamics and plasma-dynamics (although in these cases the level of nonlinearity is certainly lower, as we shall discuss later.) In the spectral methods representation, the nonlinear terms become a convolution and there are several transform-based techniques for evaluating it efficiently [62,67,68]. Numerical problems might arise from the so-called aliasing instabilities [69], since any product among fields creates new k-components (spectral leakage). In discrete Fourier transforms, alias happens when the size N k of sampled Fourier modes is less than the effective size of modes. In such cases, due to the cyclic assumption in discrete Fourier transform, a Fourier mode with wavenumber out of the size-range is aliased to another wavenumber in the domain. In a time-evolving system, this aliasing creates numerical instabilities. The importance of eliminating the aliasing error (dealiasing) has been studied since pioneering works by Orszag et al. [70]. As follows, we briefly mention this problem, well-known in hydrodynamics, and propose a similar "cure" to this numerical instability in the case of gravitational dynamics.
Let consider the simplest possible nonlinearity, i.e. the product of two functions f (x) and g(x) in a 1D, periodic domain (the 3D extension is straightforward). The truncated Fourier series with a set of {N k } modes are where f p and g q are the complex Fourier coefficients, as defined in (20). The product of the above two functions Q(x) = f (x)g(x) becomes a convolution product in Fourier space, namely where the last sum is the convolution on all the possible modes p and q of the set {N k } that satisfies the resonance k = p + q. It can be clearly seen that the product (21) generates high order harmonics with respect to the Fourier series truncated at N k . They will contribute to the well-known aliasing error [71]. The dealiasing technique consists of making the maximum number of Fourier modes k * in the above centered convolution free of aliasing, by extending the vector on both sides by zeros, to a larger vector of size N k (zero padding). Analogously, it is easy to demonstrate that the full dealiasing can be achieved by defining a k * in such a way that all coefficients with p, q, ... > k * are zero. For a quadratic nonlinearity like in equation (21), it has been demonstrated that is sufficient to filter out modes with k > k * = 2N k /3 = N/3. This fully eliminates the aliasing instability [62] and is often called the 2/3 rule, since the number of dealiased modes is 2/3 of the total number of modes. Compared to the most widely inspected cases, such as the Burgers equations and Navier-Stokes, gravitation involves high-order products and derivatives. In particular, from a quick analysis of the ADM equations, it is clear that one has products of order p = 4, or more. The nonlinearity is not completely well-defined because of contravariant tensors which complicate the degree of the nonlinearity, involving inverse computation of  (23) and (24). The figure shows a section of the Fourier space, where all the |k| ≥ k * are truncated. The Heaviside filter is a unit step-function, while the smooth filter goes from 0 to 1 more gently. In the example it has been used a = 10 to emphasize the smoothness of this filter. In all our simulations we set a = 20.
the metric. In these cases it is still possible to suppress the aliasing instability by finding an accurate truncation [69], which establishes a direct comparison between spectral methods and a full Galerking truncation. The above quadratic dealiasing described before can be extended to simple polynomial nonlinear terms of order p, by setting k * = N 1+p , therefore reducing the number of degrees of freedom in the Fourier space. In this work we will inspect the role of the aliasing truncation and also the shape of the filter [72,73]. The technique is very simple to implement: on the final product it is enough to set Q k = 0 for k > k ⋆ . In the case of the Einstein field equations, it is important to consider the very high nonlinearity of the system, where one has products of the type: The above quantity has a nonlinearity of order 5, which corresponds to a challenging convolution in the Fourier space. It is easy to envision a process in which the above products produce immediately high-order harmonics and therefore a pronounced aliasing instability. In order to suppress the above numerical pathology, different values of k * have been chosen, depending on the difficulty of the simulation and on the initial data. Generally, by filtering high Fourier modes, the price to pay is the loss of effective resolution (information). However, by suppressing these high-k's couplings, the codes become more stable and accurate. Apart of the mobile cutoff at k * , in this paper we propose two different shapes of the spectral filter, described as follows. For each representation, we define where Φ k * (k) is the spectral anti-aliasing filter. The first type, is the classical spherical truncation technique, namely a "Heaviside" filter (H) defined as: A second smooth filter (S), usually adopted for treating discontinuities and singularity formation [73,74], is given by where ξ = | k|/k * and a is a free parameter. It's clear that in the limit a → ∞, the two filters become identical. A cartoon of the two filters is represented in figure 1. Generally, in our numerical experiments, we will use k * = ∞ (no-filter), k * = N/2 (grid-size), k * = N/3 (typical quadratic nonlinearities), N/4 and so on. For the smooth filters Φ S k * (k) in equation (24), we set a = 20. The main challenge with using spectral methods for numerical relativity is the fact that they require sufficiently smooth solutions. A spectral decomposition can provide derivative evaluations with exponential accuracy (that is, higher than any polynomial power of the grid spacing) only if the data being represented is C ∞ . When derivatives are continuous at some lower order, then the error might rise too much. The filtering techniques might help to suppress the above problem. However, the accuracy of the method needs to be verified via a convergence analysis, especially when evolving nonsmooth solutions. These numerical tests are reported in the Appendix.
Note that even flux-conservative, Godunov-type schemes, although they formally solve the ideal equations, retain a considerable amount of artificial dissipation in the flux-reconstruction scheme [75,76]. The dissipation of these space-filters is sometimes difficult to quantify. In our much simple pseudo-spectral approach, the algorithm is stabilized by definite spectral kernels, that act only at very small scales, over a controlled and limited band-width of modes.
The above pseudo-spectral technique, used for the products and the spatial derivatives, has been associated to a simple second-order Runge-Kutta method for the time integration [77]. At each time discretization n, the function is advanced in time by using where f n ≡ f (x, y, z, t n ) is any of the gravitational fields at the time t n = n∆t, F is the time-derivative defined by the 3+1 governing equations, and ∆t is the timestep, chosen to be sufficiently small. It is worth saying that second-order evolution might be too dissipative, especially in very stressful conditions such as the inspiraling of binary systems, while it might still manifest a good "convergence". In the Appendix, we report on the adaptation of the current method to fourth-order Runge-Kutta, which is surely more precise, but is more numerically expensive. Moreover, in a pseudo-spectral algorithm, where the number of allocated arrays is large, memory-saving becomes a challenge. Hereafter, we will therefore show the main results for the second-order scheme, which has the advantage of being simpler, more dissipative for the highest frequencies and adaptable to other algorithms such as the implicit Crank-Nicholson (presented in Section 5). In our integration scheme, the time-step can be changed during the evolution, with a simple predictor scheme, as described in the following subsection.

A control on the stability of simulations
In numerical methods, the Courant-Friedrichs-Lewy (CFL) condition is a necessary stability condition for the solution of certain partial differential equations problems [78]. In the simplest case of the motion of a wave traveling through a discrete spatial grid with speed v, the CFL condition imposes ∆t < C∆x/v, namely a condition on the characteristic time of integration. In the above condition, C is the Courant number and is in general a number smaller than one. Numerical simulations of purely hyperbolic systems must satisfy this condition. The case of numerical relativity is more complicated because of the strong nonlinear couplings and because of the low level of hyperbolicity of the 3+1 models (the strictly-hyperbolic nature of 3+1 models is still an open issue [79]). The above complications might lead to stability conditions that are more restrictive than a simple CFL. A general criterion for the determination of the time-step integration ∆t is therefore necessary in order to have a good balance between time-efficiency of the integration and avoidance of nonphysical solutions. Here we re-elaborate the CFL idea, in a more general sense, as follows.
For any ADM or BSSN dynamic field f , one can estimate a global characteristic time as where δf = (f − f ) 2 is the rms of the field and . . . indicates a volume average. The above quantity can be interpreted as a typical dynamical time of f , evaluated as the ratio between the level of fluctuations and its typical time-derivative, similarly to plasma characteristic times in turbulence [80]. The characteristic time in equation (26) is defined at any given time of the foliation and does not involve any complicated computation in explicit time-integration schemes. For any dynamical field f , one has ∂f ∂t ≃ ∆f ∆t which, from the time discretization in equation (25), corresponds to F t n + ∆t 2 , f n + ∆t 2 F (t n , f n ) . From the latter approximation, at the intermediate step, for example, one can estimate the dynamical time in the equation (26). Note that one can also build local versions of this time, or can measure the minimum time by finding the maximum of the denominator.
In analogy with the CFL condition, it is possible to use the characteristic time T f for guessing an appropriate ∆t. If the system undergoes fast evolution, T f will drop off and the algorithm can self-adjust ∆t. This might indicate a "precursor" of a simulation blow-up, as we shall see in our numerical experiments. Throughout this work, we will introduce the Running Stability-Check (RSC), which is a method that consists in a continuous adjustment of the Runge-Kutta time-step, based on the measure of T f in equation (26), such that Here C is similar to the Courant number. The method simply corresponds to a adaptive time-refinement, typical of many schemes such as the leapfrog integration, except here is constrained with the minimum characteristic time given by (26). This strategy is very similar to the adaptive stepsize control for Runge-Kutta, described in standard numerical textbooks [77]. The general purpose of these classical methods is the achievement of some predetermined accuracy in the solution, with minimum computational effort. Here we re-elaborate this technique, but instead of establishing a desired error, we define a global time T f for the system and keep the integration step always smaller than this dynamical time. Note also that the above technique might provide a general condition not only in the case of Runge-Kutta methods, but in any explicit time-integration scheme.

Standard numerical testbeds
In this section we will perform direct numerical simulations in order to validate our numerical procedure, via standard gravitational testbeds [22,43,44,47]. For each initial data, we will check the accuracy by inspecting the conservation of several quantities, for both the ADM and BSSN formalism. Regarding the spatial integration, we will show how to stabilize the pseudo-spectral code by varying our dealiasing filter in the projection (22), by changing both the k * and the filter shape. For the time integration, in each test, will apply the RSC by monitoring the characteristic times in (26). After these typical, fundamental tests, we will show a new possible initial data that leads to stationary, large amplitude, nonlinear waves. These waves bounce with growing amplitude and therefore represent another good test for numerical relativity. We will then move to more stressful conditions, such as the Gowdy spacetime (both forward and backward in time) and the Schwarzschild equilibria. All the simulations performed are summarized in table 1, where we report all the main parameters used, such the type of the initial condition, the amplitude of the perturbations, the number of points of the spatial grid, the k * of anti-aliasing filter, the type of filter, the formalism adopted, and so on. In the last column, for the reader's convenience, we also report the stability of the run.

The gauge wave
The gauge wave test is based on a gauge transformation of the Minkowski spacetime.
As suggested by Ref. [22], the metric is given by   where describes a sinusoidal modulation of amplitude A, propagating along the x-axis. Since derivatives are zero in the y and z directions, the problem is essentially 1D. The metric in equation (28) implies β i = 0, and K xx = −∂ t γ xx /2α. For the spatial metric γ ij one gets while the only nontrivial component of the extrinsic curvature is It is easy to demonstrate that the above fields satisfy the initial data constraints in equations (2) and (3). In order to stress the stability of our new code, we solved the ADM equations (4)- (6), with both small-and large-amplitude initial conditions [43]. We chose a spatial domain x ∈ [0, 1] with a number of meshes N x = 64. The initial time-step is ∆t = 5 × 10 −3 and, according to the literature, we used the harmonic slicing, namely f (α) = 1 in equation (4). For the first run, RUN 1 in table 1, we impose A = 10 −2 and for the numerical technique we choose the filtering mode k * = N x /3, with the smooth filter described by equation (24). In order to test the RSC method described in section 3.1, the run was performed twice: with a fixed time-step (RUN 1 ) and with an adaptive-timerefinement suggested by the RSC method (RUN 2 ). Figure 2 (a) shows the evolution of the characteristic times given by equation (26). Note that only the characteristic times associated to the non-null fields have been reported, namely T Kxx , T γxx and T α . The full (black) line represents the time-step ∆t of the simulation. One can immediately see that when a dynamical time T f becomes on the order of the time-step, instabilities arise and the code crashes. In the second test, reported in the panel (b) of the same figure, we repeated the experiment with the RSC switched on. In particular, we set the adaptive time refinement in equation (27) with C = 1/6, and the time-step has been halved any time that the condition is violated. The run is therefore stabilized and carried out until 250 crossing times. This simulation is reported as RUN 2 in table 1. Similar behavior has been obtained for other initial data (not shown here).
For the large-amplitude regime of the gauge-wave test, we have chosen A = 0.9 and A = 0.96, with a spatial grid along the x-axis of N x = 128 points and ∆t = 10 −3 . The other parameters are the same as in the small amplitude case, including the slicing gauge. Again we have performed two experiments: first without anti-aliasing filter (RUN 3 ) and second with a k * = N/3 smooth filter (RUN 4 ), defined in equation (24). A comparison between the waveforms of the extrinsic curvature trace K, in both cases, has been represented in the top-row of figure 3. As it can be seen, there is a net improvement with the anti-aliasing filter. One can observe, overall, an excellent agreement between the exact and the numerical solutions. In the unfiltered case numerical instabilities arise, as can be seen in panel (a). The analogous filtered simulation is reported in (b).
A second, more challenging test has been performed, with a larger perturbation, i. e. A = 0.96. We found that now a stable solution can be achieved even with k * = N/2.5. In order to test the different anti-aliasing filters we performed the test twice: RUN 5 uses the Heaviside filter defined in equation (23) while RUN 6 adopts the smooth filter (24). The panel (d) of figure 3 shows the L 2 errors of the x−momentum constraint, where L 2 norm is given by where dΩ |γ| is the volume element. The smooth filter again stabilizes the code. In figure 3 (c) we report the waveform of K at time t = 10, which is in very good agreement with the exact solution. It is important to emphasize that the last set of runs, where A = 0.96, is considered to be a very strong stress for codes, as discussed in [43]. For this testebed, we obtained similar results for the BSSN code. Our conformal model follows the approach by Zlochower et al. [20]. We used A = 10 −2 (RUN 7 ) and A = 10 −1 (RUN 8 ), a spatial grid with N x = 128 mesh points and a time-step ∆t = 5 × 10 −4 . The harmonic gauge in equation (6) has been used. In order to better stabilize the simulation, a smooth filter with k * = N/3 has been chosen, for both amplitudes, similarly to the ADM cases. The left panel of figure 4 shows the L 2 errors of the Hamiltonian for the two amplitudes, while in the right panel we compare the wave-form of γ xx , at t = 500. At this time, the norm of the Hamiltonian violation ||H|| 2 is about 0.1% of A, for the case with A = 10 −2 . The run shows an excellent agreement between numerical and exact solutions.
In agreement with the existing literature (see for example [21]), we found that generally the BSSN formalism is less able to handle gauge waves compared to the ADM  [20,81], although the Hamiltonian error shown here is at least one order of magnitude smaller than the one presented in the above works. This good performance is evidently due to the precision of the spectral representation, accompanied by the RSC and the dealiasing filters.

The robust stability test
The robust stability test is a classical initial condition for numerical relativity, particularly stressful from the numerical point of view, that makes use of random, small perturbations in the initial data. This testbed is able to initiate numerical instabilities that might exponentially grow in time, arising from small scales and propagating back ‡ This difference is mainly due to truncation and roundoff error propagation. The BSSN formalism involves a higher number of products than the ADM. Just as an example, by using double precision (that we use throughout the paper), even the initial Hamiltonian is already ∼ 10 −16 . We checked that by switching to quadruple precision the errors reduces to H ∼ 10 −32 . to larger structures. These growing modes are quite dangerous since they might be confused with physical mechanisms. It is therefore important to inspect these numerical instabilities in order to prevent numerical artifacts.
The initial configuration consists of small random perturbations to the Minkowski space, where ds 2 = −dt 2 + dx 2 + dy 2 + dz 2 , with the unperturbed metric η ij . The numerical noise is distributed among every field. Since we base our numerical algorithm on the use of FFT's, we vary several resolutions, keeping the number of meshes as a power of 2. In order to follow the same configuration of Ref. [22], we chose a spatial domain x, y, z ∈ [0, 1], a number of grid points N x = 50ρ n , where the mesh size is therefore ∆x = 1/(50ρ n ) and the initial time-step is ∆t = 0.01/ρ n . Here we vary ρ n = 64n/50, with n = 1, 2, 4. To speed up calculations, we use only four grid points in the y and z directions. We solve the equations by using the harmonic gauge, i.e. f (α) = 1 in the slicing equation (6). In the ADM tests, the initial conditions are given by where ε f are random independent perturbations with ε f ∈ −10 −10 /ρ 2 n , 10 −10 /ρ 2 n .
In our robust stability tests we vary the resolution (n), the type of the dealising filter (Heavyside or smooth), and we vary from fixed time-step to the RSC. The goodness of our numerics has been monitored by looking at the evolution of the constraints in time, such as the L 2 -norm of the Hamiltonian constraint.
In figure 5 (a) we report the Hamiltonian error, for different resolutions, namely for n = 1, 2, 4. No anti-aliasing filters have been used for the runs in this panel. As it can be seen, the code becomes suddenly unstable. These runs are reported in table 1, labeled as RUN 9 -RUN 11 . In order to show the benefit of the anti-aliasing filters and the robustness of the code, a second set of tests has been performed, by using filters at k * = N/4, but without RSC. A time-step of ∆t = 10 −3 and n = 1 has been chosen. The L 2 -norms of the Hamiltonian constraint are shown for these cases in figure 5 (b). The simulations show that, with the "Heaviside" filter in (23), the code becomes now stable until t ∼ 1000, while is even more stable using the smooth filter. These runs are reported in table 1, as RUN 12 and RUN 13 . We complete the robust stability test campaign with the BSSN approach and finally introduce the RSC condition. The parameters are the same of the ADM previous tests. For the slicing condition we used the harmonic and hyperbolic Gamma-driver slicing in equations (15) and (16). The tests are carried out for t = 10 4 crossing times, at most. First we compare the two filters described by equations (23)-(24), but keeping a fixed time-step of integration. In the bottom left panel of figure 5, we show the ||H|| 2 for the two filters. Note that using the Heaviside filter the code crashes at t ∼ 2000 crossing times (RUN 14 ), while using the smooth version it remains more stable (RUN 15 ). In panel (d), finally, we show the same tests with the RSC, by using an adaptive time refinement with C = 1/4, as described in section 3.1. In this case, the full algorithm stabilizes the simulation until t = 10 4 crossing times (RUN 16 -RUN 17 ). Stability of the run has been achieved with the combined use of filters and adaptive time-integration.

New initial data: the standing waves
Here we propose a very simple, 1D, initial set of data that satisfies equations (2)- (3) and that can be used to test codes in numerical relativity, with various amplitudes. This initial condition consists of a simple initial data for the extrinsic curvature only, on an initial flat metric, that suddenly produces some non-propagating, oscillatory modes. The physical pattern consists of a "bouncing metric". This solution to the constraint conditions is not supposed to have an immediate physical application, as in many of the other testbeds, but only for testing purposes. However, it might have some applications and similarities with cosmological models of inflation. For this regard, see for example Adamek et al. [82].
The initial spatial metric γ ij obeys Minkowski flat space, which means zero curvature and where Ricci tensor vanishes (R = 0). In this case, the Hamiltonian and the momenta constraints reduce to respectively. We now chose an initial perturbation only for the extrinsic curvature K ij , with a Minkowski flat initial metric. With this choice, one has: and K ij = 0 otherwise, which satisfies the set of constraints in equation (32). Unlike the previous case of gauge waves, this condition does not lead to a wave propagation in space but to standing waves of the metric. One can stress the code by increasing the amplitude of the extrinsic curvature perturbation A, as follows.
In figure 6 we report the time evolution of the metric, for several ADM tests. For RUN 18 , a small amplitude A = 10 −1 and no anti-aliasing filter has been used. We chose a spatial domain with x ∈ [0, 1] and a spatial grid N x = 64. The time-step is ∆t = 10 −3 and the harmonic slicing [f (α) = 1] has been chosen. The simulation has been carried out until 100 crossing times. We have carried out a second test (RUN 19 ) with same parameters but choosing a smooth filter with k * = N x /3, as reported in figure 6 (b). It is evident that the filter improves the stability, in fact this run is carried out until 1000 crossing times. The top panels of figure 6 show that the Hamiltonian constraint remain constant in both cases, but the x−momentum constraint (for the symmetry of the problem the y−momentum and z-momentum are zero) grows without filter and remains constant with k * = N/3. In order to test a more nonlinear case, a second, high-amplitude condition has been used, where we chose A = 8, a spatial grid with N x = 128 points and a k * = N/3 filter. All the other parameters have been chosen as in the previous tests, as summarized in table 1. The higher amplitude of A now induces asymmetry in the standing wave. The metric starts to bounce and becomes asymmetric. Again, this test has been performed twice: a first test (RUN 20 ) using the standard truncated filter described by equation (23) and a second one (RUN 21 ) with the filter in equation (24). The time-evolution of L 2 norm of x−momentum constraint, computed using the definition in equation (29), is shown in the panel (c) of figure 6, for both cases. It is clear that the smooth filter works better than the Heaviside filter.
We ran the same tests for the BSSN code (not shown), setting A = 10 −1 , a spatial grid with N x = 64 points, a time-step ∆t = 10 −4 and the "1 + log" slicing (i.e. f (α) = 2/α in equation (6) and β i = ∂ t β i = 0). The tests are carried out until 100 crossing times. Similarly as in the previous ADM tests, the smooth filter better stabilizes the code.

The gravitational wave test
In this subsection we test our two codes, the ADM and the BSSN spectral algorithms, via linear gravitational waves. In particular we check whether they are able to retain accuracy, for long spacetime travels, by comparing amplitudes and phases with the wellknown analytical solution. This solution to the linearized Einstein field equation can be written as where b = A sin 2π(x−t) d and d = 1 is the linear size of the propagation domain. The wave amplitude is small, in order to satisfy the linear regime conditions. The metric in equation (34) represents the transverse-traceless gauge in which the wave amplitudes are purely spatial, with null trace and transverse to the propagation direction. It describes a wave propagating along the x axis, with the polarization aligned with the y and z axes [83]. This metric is written in Gauss coordinates, that is α = 1 and β i = 0. Using this metric, the extrinsic curvature reduces to K ij = −∂ t γ ij /(2α), so that the nontrivial components reduce to K yy = 1 2 ∂ t b and K zz = − 1 2 ∂ t b. Since this wave propagates along the x-axis and all derivatives are zero in the y and z directions, the problem is essentially 1D. We chose a spatial domain x ∈ [0, 1] on a spatial grid with N x = 64, with an initial ∆t = 10 −3 and the harmonic slicing. We use a very small amplitude, namely A = 10 −8 , in order to satisfy the linearity of the evolution [43]. The evolution is carried out up to t = 1000 and a smooth filter, described by equation (24), with k * = N/4, has been used. Figure 7 (a) represents the time evolution of the L 2 norm of ADM constraints. The errors are well bounded. In the panel (b) we compare the wave-form of K zz at the final time, i.e. at t = 1000, showing overall an excellent agreement between numerical and exact solutions. This run is reported in table 1 as RUN 22 . We performed an analogous test with the BSSN formalism, choosing N x = 128, ∆t = 5 × 10 −5 , with the harmonic slicing and the "Gamma-driver" shift condition. The evolution is carried out again for a time of t = 1000 and the smooth filter with k * = N x /8 has been used. This run is labeled in table 1 as RUN 23 . The bottom left panel of figure 7 shows the L 2 Norm of the Hamiltonian, and the bottom right panel of the same figure shows the excellent agreement with the analytical solution for the waveform of K zz , at the final time. It is important to stress that, in the BSSN runs, the choice of a stronger dealiasing filter has been used in order to obtain a longer time, stable and accurate simulation. We also point out that here we used a finer grid, so that the simulation corresponds essentially to the ADM runs. More points are now needed since the BSSN is less accurate. As discussed in previous sections, this is related to the fact that BSSN is a more elaborated formalism, where a larger number of products leads to higher truncation errors and therefore to less precise solutions [22]. In summary, the gravitational wave test therefore is well solved via both algorithms, the smooth filter stabilizes the code for very long times, giving a very good representation of the solution.

The Gowdy spacetime
All the tests described so far considered perturbations to the flat metric. In this section we deal with a genuinely curved exact solution, known as the polarized Gowdy spacetime [20,81]. The Gowdy waves are vacuum cosmological models used to test codes in a strong field context and present a serious challenge for numerical relativity. These particular solutions of the vacuum Einstein equations on the 3-torus describe an expanding/collapsing universe, containing plane polarized gravitational waves where Here λ and P are functions of z and t only. The time coordinate t is chosen such that time increases as the universe expands. The metric is singular at t = 0 which corresponds to the cosmological singularity. We will carry out our tests in both time directions, i.e. in the collapsing and expanding dynamics. With the metric described by (35), the Einstein equations lead to [22,84] P (z, t) = J 0 (2πt) cos(2πz), (36) λ(z, t) = −2πJ 0 (2πt) J 1 (2πt) cos 2 (2πz) + 2π 2 t 2 J 2 0 (2πt) + J 2 1 (2πt) where J 0 and J 1 are Bessel functions. Equations (36) and (37) yields a diagonal spatial metric of the type To complete the initial data conditions, the extrinsic curvature is obtained as a time derivative of the metric. Since the metric (35) depends only on the z−direction, the spatial evolution is 1D. We chose a spatial domain z ∈ [0, 1], with N z = 50ρ n , ∆z = 1/(50ρ n ) (see before), and a ∆t = 5 × 10 −4 /ρ n . In all the Gowdy tests, we used a k * = N z /5, smooth, anti-aliasing filter.
In the expanding form, the time coordinate t increases during the evolution of the simulation, P approaches zero asymptotically and λ increases linearly. Since the metric is singular at t = 0, we set the initial data from the exact solution at t = 1 and proceed forward in time. Due to the exponential growth, such evolution may crash rather soon but it will test the accuracy of a code in a rather harsh situation. In figure 8 (a) we report the L 2 norm of the Hamiltonian, for the three resolutions n = 1, 2, 4 (RUN 24 -RUN 26 ). The errors are bounded in all the tests and our results are in good agreement with the basic literature [85,86]. Similar results have been found for the ADM simulation.
In the collapsing version, the time coordinate in the Gowdy metric can be transformed so that the initial singularity is approached asymptotically, backwards. The new time coordinate τ is defined by τ = 1 c ln(t/k). In other terms, we replace t → ke cτ and the time-step ∆τ is chosen to be negative. We choose a particular value of t 0 such that the initial slice is far from the cosmological singularity and, following the standard approach [22], we chose the twentieth zero of the Bessel function, so that J 0 (2πt 0 ) = 0. Finally we get c ∼ 2 × 10 −3 and k ∼ 9.6 × 10 13 . In the panel (b) of figure 8, we present a comparison between the ADM (RUN 27 ) and BSSN (RUN 28 ) runs, for the collapsing direction. The ADM code shows a smaller error on the Hamiltonian error, but both schemes proceed until 1000 crossing times (τ = −1000). This evolution is less challenging than the expanding case, since the lapse function is essentially an exponential in τ so that the spacetime is becoming less dynamical. On the other hand, the value of the physical metric γ zz at t = t 0 is of the order of 10 3 . Overall, as it can be seen, also in this more challenging test, the pseudo-spectral code is able to handle the numerical evolution, with violation errors that are comparable (or smaller) than in the existing literature [22].

The Schwarzschild black hole
In this classical test we consider the stability of an isolated, Schwarzschild black hole (BH), in a fully 3D spatial geometry. The Schwarzschild metric is used as initial data to test the ability of the code to evolve BH spacetimes within the puncture approach [57,[87][88][89]. The metric in isotropic coordinates can be written as [22,58]   (b) Conformal factor χ = ψ −4 , at different times of the BSSN simulation, with the same parameters as RUN 32 , but higher resolution 128 3 . A smooth filter k * = N/2 has been used.
where r is the isotropic radius, m is the mass of the black hole and ψ = 1 + m 2r is the conformal factor. The initial spatial metric and the extrinsic curvature become respectively. The lapse is initially set to one and the shift is vanishing everywhere. For a Schwarzschild black hole both the linear momentum and the dimensionless spin are set to zero [90]. We set the puncture mass m = 1.0M and the runs are carried out until t = 1000M, corresponding to 1000 crossing times, unless the code becomes unstable. The spatial domain is x, y, z ∈ [0, 300], with an isotropic spatial grid N = 64 3 and a time-step ∆t = 10 −1 . The 1 + log gauge condition is now used [i.e. f (α) = 2/α in equation (6) and β i = ∂ t β i = 0], for both the ADM and BSSN tests. In order to handle Schwarzschild equilibrium, since the puncture limits differentiability of the solution, we shifted the position of the singularity by half grid-step, namely ∆/2, for each Cartesian direction. In this regard, it is worth saying that exact spectral projections always allow to obtain the field in the continuum via a simple phase-shift of the Fourier coefficients [91].
In the panel (a) of figure 9 we report the evolution of the Hamiltonian error ||H|| 2 , for the ADM (RUN 29 and RUN 30 ) and the BSSN tests (RUN 31 and RUN 32 ), with and without the smooth dealiasing filter. Here we found that k * = N/2 is enough to stabilize the simulation. As expected, and as well known in the literature [46,58], the ADM formalism is not well suited for the simulations of singular objects such as BH's. Even if the ADM "beats" the BSSN formalism in several numerical challenges (see previous sections) regarding the precision, the BSSN remains one of the most stable numerical approaches for massive stars dynamics.
In figure 9 (b) we report the conformal factor χ = ψ −4 for a simulation with the same parameters of RUN 32 but with higher resolution (128 3 ). We report a 1D cut through the center of the domain. The simulation shows slightly higher non-conservation (probably due to the choice of the time-step), but it is extremely more precise, as we shall discuss later in the convergence tests of the Appendix. For this version of the algorithm, namely the filtered BSSN, the simulation is stable (due to the constraint propagation properties of BSSN), accurate and convergent.
As is clear from figure 9, the BSSN is able to simulate black holes dynamics, accurately, for long times. Even with very small dealising filters, with k * = N/2, the BSSN retains the black hole shape and keeps a very small error on the constraints. However, the small violation of the periodicity, together with the high accuracy of the spectral code, might produce wiggles that propagate from the discontinuous boundary. In the next section, we will deal with these effects and cure this pathology by introducing a new numerical strategy, where we still keep the high precision of the pseudo-spectral strategy.

The implicit hyperviscous boundary method
We study black hole dynamics via a head-on collision of two massive objects. We propose a new method in order to treat the boundary conditions with our pseudo-spectral, periodic code, based on the use of absorbing hyperviscous boundaries. As demonstrated in previous sections, our pseudo-spectral technique, together with an appropriate choice of the dealiasing and the time-step, is able to handle all the classical numerical testbeds. The main advantage of the pseudo-spectral code is of being relatively fast (since it stems on FFTs) and incredibly precise (since it relies on an exact decomposition). It would be hard to beat the precision of such methods with finite difference or finite volume/elements techniques, at a given resolution [70]. On the other hand, it is evident that the main weakness of the method relies on the boundary conditions that need to be periodic for this particular spectral decomposition. Periodic boundaries are not always a good choice for general relativity problems, as for example in the case of compact object mergers. Here we present an additional algorithm that helps to overcome this limitation, which is based on the use of dissipating boundaries, combined with implicit time integration. This Implicit Hyperviscous Boundary method (IHB) has been tested with gravitational waves and compact objects dynamics. The method is relatively easyto-implement and efficient.
The following IHB method has been inspired by works in hydrodynamics and magnetohydrodynamics, with applications to turbulence and to stellar dynamo. It is indeed possible to use Cartesian geometries in order to simulate spherical domains, using for example volume penalization methods [92,93], or by using hyperviscous forces [94], starting from a given distance from the center of the Cartesian lattice. In this way it is possible, for example, to embed a spherical domain inside a cube, by dissipating the regions close to the corners, with a radial envelope. Here we re-elaborate this idea and apply it to numerical relativity, in a pretty straightforward way.
Hyperviscoscous terms of the type ν n ∇ n f , with ν n a numerical coefficient, are able to dissipate quickly ripples and numerical artifacts, but they might reduce dramatically time-step since they involve high-order derivatives and fast diffusion. For the method presented here show results with fourth-order dissipation (but the technique can be extended to any order.) We make use of a technique for the hyperviscous terms based on the implicit, second-order Crank-Nicholson method [95]. Hereafter we will apply this technique only to the BSSN code, which is a stable method for solving problems such as the binary mergers and black holes dynamics.
For any variable of the BSSN method, we advance in time the solution f {ideal} (x, t) by using the second order Runge-Kutta described in equation (25). Simultaneously, we evolve in time an accompanying, twin field f {H} (x, t) which obeys to the same BSSN equation but is also subject to hyperviscous dissipation, being therefore highly damped.
We advance in time this hyperviscous field f {H} (x, t) by using a Crank-Nicolson, semiimplicit method (CN). The CN integration can be easily matched with the Runge-Kutta method and it is numerically very stable with respect to viscous dynamics. For any BSSN hyper-dissipated variable we have an evolution equation where N indicates the BSSN evolution equation, ν 4 is the hyper-viscous coefficient. In the case of the "ideal" BSSN set of equation, ν 4 = 0 and taking the Fourier transform, from the typical Runge-Kutta 2-steps method one gets At this point, for any damped field, it is easy to get the time-splitting, second order CN formula, which is stable and particularly easy to implement in pseudo-spectral codes. For all the BSSN variables, we advance in time equation (39) for both the ideal fields f (ν 4 = 0). By taking the inverse transform now one has both fields f {ideal} (x, t + ∆t) and f {H} (x, t + ∆t) and then has to deal with the solutions match.
The scope is to suppress spurious effects only at the boundaries, therefore we interpolate between the two solutions gradually (or in a Heaviside way) going towards the borders, only in a region close to the boundaries. With this procedure we retain the pure, ideal dynamics of the BSSN formalism in the center of the box. Overall, the price to pay for the IHB is that we have to double the time of integration, slowing down the computation. On the other hand, the method gives clear benefits on the stability and the goodness of the solutions, as we shall see in the next subsections, where we present two tests.

Gravitational wavepacket absorption
In order to test the IHB method, we check whether the algorithm is able to absorb fluctuations that move outward, from the center of the domain toward the boundaries. At the initial time, we build a gravitational wavepacket, namely a windowing of the gravitational wave at the center of the domain via a Gaussian-shaped filter, as shown in figure 10-(a). As in the linearized gravitational wave test (see section 4.4), the metric is given by equation (34). We initialize the perturbation b(x, t = 0) as a wavepacket that propagates along the x−axis, with a Gaussian window of the type exp[−( x−x 0 σ ) 2 ]. We set an amplitude of A = 10 −5 , the wavenumber k = 40 and the width of the window σ = 0.05 is small enough in order to localize the waves in the center, namely In order to match the two solutions f {ideal} (x, t) and f {H} (x, t) in the IHB technique, we chose a linear matching when distance approaches to the borders. In particular, at the end of each time-step, we impose for the BSSN variables In the above matching-condition we chose λ = L 0 /4, namely one quarter of the box. In this test we set the hyperviscous coefficient ν 4 = 10 −7 . This test is labeled as RUN 33 in table 1. With our IHB approach, when the packet is far from the boundaries, the BSSN evolution equation f {ideal} (x, t) dominates in equation (40) -all the field equations evolve in the standard way. Otherwise, when the waves approaches to the boundaries, the hyperviscous terms of the type ν 4 ∇ 4 f dominate, and the wavepacket is totally absorbed.
The time history of the propagating gravitational wave is reported in figure 10. As it can be seen, the modulation moves coherently from the center of the domain, toward L 0 = 1. When the wavepacket enter the region of the viscous boundary, at x = 3 4 L 0 , it starts to "feel" the hyperviscous damping in equation (38). Before approaching to the borders, the wave is totally absorbed. Similar results (not shown here) have been obtained with other matching conditions, discussed in the next subsection. The method can be easily exported to 3D gravitational dynamics, in order to suppress boundary effects, as we will see in the next, final experiment.

Head-on collision of black holes
In this last section we test the IHB technique via one of the most stressful conditions, namely the head-on collision of two black holes. In this test, the initial conditions consist of two non-rotating, massive compact objects. We test this black hole crash only via the most numerically stable approach, namely the BSSN spectral code, with a smooth filter in equation (24) and with the RSC active.
The parameters of the initial conditions are the following. The three-dimensional computational domain is x, y, z ∈ [0, 25] and flat Minkowski spacetime is imposed as a boundary condition everywhere. We chose two identical black holes with masses m 1 = m 2 = 0.5M, and C 1 = {10.5, 12.5, 12.5} and C 2 = {14.5, 12.5, 12.5}, where the C j are the locations of the two punctures. These parameters correspond to an initial separation of the BHs a little larger to that of an approximate innermost stable circular orbit (ISCO), as determined in [96]. The two punctures consist of two Schwarzschild solutions with initially zero linear momentum and spin, as described in section 4.6.
As suggested by [43], the lapse is initially set to: where r * d = (r 4 + 10 −24 ) 1 4 and r is the coordinate distance of a grid point from the puncture. We present two types of runs for the head on collision. We evolved a first simulation (RUN 34 ) without the IHB method. We carried out the simulation until the code crashed, that is a few times after the collision. The evolution of the contour surfaces of the lapse α are reported in the top row of figure 11. The panels (a), (b) and (c) show the x − y section of the 3D domain, at z = L 0 /2, i.e. at the middle of the z−domain, at times t = 1, 13 and 25M, respectively. As it can be seen, after the BH's collision, spurious harmonics develop and grow in time. These numerical disturbances are related to gravitational fluctuations, generated at the central collision-region, that propagate away and interact at the periodic borders. This non-physical interference is also amplified by the violation of the mean gravitational profile of the black holes at the border, which eventually leads the simulation to crash at t ∼ 30.
In order to stabilize the spectral algorithm, a second simulation (RUN 35 ) has been carried out by switching-on the IHB technique. In order to quickly truncate noise that periodically travels through the periodic boundaries, we match the two solutions, f {ideal} (x, t) and f {H} (x, t), with a rigid condition where x 0 is set at the center of the 3D spatial domain, and λ = 3 8 L 0 . In the bottom row of figure 11 we represent the contour surfaces of the lapse [panels (d), (e) and (f)], at the same times of the previous run. The simulation has been noticeably improved, thanks to the diffusive boundaries that absorb the spurious periodic effects, as well as the out-propagating gravitational fluctuations. The line contours reveal a very smooth and well-behaved profile of the final (merged) black hole, confirming that dealiased pseudo-spectral methods can represent a good numerical strategy for the simulation of compact object dynamics.
Finally, we compare the violation of the Hamiltonian constraints for the two runs (with and without IHB), as reported in figure 12. As it can be seen, the run with the viscous boundaries is much more stable, manifesting also lower violation of the constraints. This is due to the fact that the boundary ripples have been absorbed and the method prevents the contamination of these wiggles in the central ideal part of the domain, where the merger occurs. A similar strategy can be used for a variety of studies, including the inspiring binaries and the multiple black holes systems [97][98][99]. There are other techniques able to handle black hole dynamics, as in the SpEC code [35,36], where the method imposes a large infrastructure able to handle the singularity. The method presented here is different and simpler than the above (more refined) strategy. In the case of massive objects, we base our numerical approach on (1) a simple grid-shift of the black hole (in the case of the static Schwarzschild), (2) the use of filtering/dealiasing techniques that retain the differentiability of the solution and (3) the use of hyperviscous boundaries which absorb perturbations.

Conclusions
Pioneering works on spectral methods applied to numerical relativity, as in the case of the SpEC and the SGRID codes, suggested that the spectral approach is very efficient. These codes, together with more common approaches such as BAM [100], the Einstein Toolkit [101,102], and more exotic codes such as Dendro-GR [103], represent today a valid tool of investigation for the understanding of gravitational dynamics. In this work we pursued the path of spectral algorithms, by examining the role of numerical instabilities and the performances of filtering methods.
We have studied the evolution of gravitational fields by solving numerically the Einstein field equations, in vacuum, with two different approaches, namely the ADM and the BSSN formalism. Our numerical method is based on a pseudo-spectral technique, where we compute spatial derivatives via simple and efficient Fast Fourier Transforms. A detailed analysis of the aliasing instabilities has been presented and two types of antialiasing filters have been used in order to suppress numerical artifacts that are due to the intrinsic nonlinear nature of the governing equations. We presented a new technique that improves the stability of the simulations, named the Running Stability Check, which monitors the strength of the time derivatives, providing an optimal time-step. The code shows a very low violation of the constraints, due to the precision of standard spectral methods, which are high-order accurate. We have found that the filter type and the cutoff wavenumber are crucial for the evolution of the gravitational equations, with very stable and accurate solutions, especially for smooth filters with cutoffs k * < N/3.
The stabilized, pseudo-spectral code can handle even more difficult gravitational problems, by solving BH equilibria and the dynamics of binary systems. In this regard, we did overcome the typical limitations of periodicity, which is needed for such pseudospectral, FTT-based methods. Our strategy, borrowed from fluid and plasma dynamics, is able to suppress spurious boundary effects. We matched the ideal solution in the center of the domain, based on the pseudo-spectral solution of the BSSN system of equations, with a semi-implicit, second order Crank-Nicholson technique at the boundaries, where we added hyper-viscous diffusion. This strategy has been successfully tested via a simple gravitational wavepacket absorption, as well as more catastrophic initial data, such as the head-on collision of black holes.
In summary, the code has been successfully tested via the classical numerical relativity testbeds. Via these typical initial data, we have validated the goodness and the robustness of our numerical strategy. Future works will concentrate on the approach of the present algorithm to the dynamics of in-spiraling binary systems, as well as to the dynamics of many-body problems. Finally we plan to export the method to the solution of the gravitational fields in presence of matter as well its coupling with electromagnetic fields, in the framework of general relativistic magnetohydrodynamics.

Appendix A.
Appendix A. 1

. The time integration scheme
The most commonly used methods to integrate in time systems of coupled partial differential equations are the well-known Runge-Kutta (RK) schemes. Several classifications of these methods can be done, according to their convergence order, the number of stages or their explicit/implicit structure. Low-order RK methods have the advantage of simplicity and the relatively low need of memory allocation. This is surely a good point for numerical relativity, especially in the BSSN formalism, where the number of tensors proliferates compared to the ADM approach. However, low-order RK, such as the second-order method (RK2) are quite dissipative, causing a departure (in time) from exact solutions (when these are known). Some of these accuracy problems are partially solved by using higher-order, such as the popular fourth-order RK4 §.
In this brief subsection of the Appendix we compare the accuracy of both RK2 and RK4, in the case of our pseudo-spectral filtered method. The RK2 has been described in Section 3, by equation (25). Analogously, for the RK4 one gets f n+1 = f n + ∆t 6 [F 1 + 2F 2 + 2F 3 + F 4 ] , (A.1) where F 1 = F (f n ), F 2 = F (f n + ∆t 2 F 1 ), F 3 = F (f n + ∆t 2 F 2 ) and F 4 = F (f n + ∆tF 3 ). In figure A1 we show the Hamiltonian error ||H|| 2 , for the most stable runs shown in table 1. For the sake of simplicity, we show only few testbeds, namely: the large amplitude Gauge wave (Run 8 ), robust stability test (Run 17 ), the gravitational wave test (Run 23 ) and the Schwarzschild black hole (Run 32 ). In all the cases, the RK4 has always a smaller violation error, as expected.
In table A1 we report the difference between the simulated profile of the field χ and its analytical solution. In particular, we computed the L 2 norm of the error ∆ a χ = χ RK# − χ analytical , (A.2) for both the second (#=2) and the fourth (#=4) order. As it can be seen, the RK4 is more accurate, although it requires extra computational time and memory. In future § Alternatively, explicit (e.g. Adams-Bashford) or implicit (Adams-Moulton) linear multi-step methods improve the stability properties, albeit at the expense in increasing the complexity of the algorithm.  works, a combination of high-order time-integration schemes with the IHB will be presented, coupling alternative techniques with implicit diffusion for the boundaries.

Appendix A.2. Convergence tests
In this second Appendix we present the convergence tests at different resolutions N, for the filtered BSSN algorithm. We start with the Schwarzschild black hole, discussed in subsection 4.6, by varying the resolution from very-low (N = 32 3 ) to moderate resolution (N = 256 3 ). For this set of runs we used the same parameters as for RUN 32 (see table  1) and by shifting the puncture position in between mesh-nodes, as described before. In order to quantify accuracy, for each resolution we computed the difference between a profile at the end of the simulation and its analytical solution, as in equation (A.2), where now ∆ a χ = χ N −χ analytical and χ N is the conformal factor for the simulation at a given resolution N. This norm can be computed for any variable but here we report results for χ (other fields give similar results, not shown). In figure A2 (a) we show this discrepancy, for different numerical resolutions, at the end of the Schwarzschild simulations (t = 100). As it can be seen (and as expected) the error diminishes with increasing resolution. To quantify the convergence, we computed the L 2 norm of this error, as reported in table A2. The global error slightly diminishes, going from low to high resolution.
As a final set of convergence tests, we repeated the head-on collision of black holes, discussed in the subsection 5.2. In particular, we performed RUN 35 , going from 128 3 up to moderately-high resolution, with N = 512 3 mesh points. We compared the simulations just after the collision of the BHs, at t = 22, reported in panel (b) of figure  A2. This sequence of runs suggests that the solution is well described within our typical resolutions, for N ≥ 128 3 . In particular, the profiles are very similar, although the highest resolutions reveal a more singular and accurate profile. After the merging, the filter, together with the good differentiability properties of the BSSN formalism [57], is able to handle the final singularity at the location of the puncture, which typically is expected to cause problems for spectral codes. To quantify the level of convergence here, we computed the error ∆ N χ = χ N 1 − χ N 2 , for neighboring spectral resolutions N 1 and N 2 . In order to perform these direct measurements, we computed the difference only on the coarser grid and we used for all the simulations the same time-step, constrained by the 512 3 run. We obtained ||χ 256 − χ 128 || 2 = 3.3 × 10 −3 and ||χ 512 − χ 256 || 2 = 9.7 × 10 −4 . This further confirms that the solution is accurate and convergent.