{'fraction_non_alphanumeric': 0.03975046480539635, 'fraction_numerical': 0.012184815844870238, 'mean_word_length': 4.293371289141928, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A spectral solver for evolution problems with spatial S3-topology

We introduce a single patch collocation method in order to compute solutions of initial value problems of partial differential equations whose spatial domains are 3-spheres. Besides the main ideas, we discuss issues related to our implementation and analyze numerical test applications. Our main interest lies in cosmological solutions of Einstein's field equations. Motivated by this, we also elaborate on problems of our approach for general tensorial evolution equations when certain symmetries are assumed. We restrict to U(1)- and Gowdy symmetry here.


Introduction
Numerical studies of initial value problems of partial differential equations on certain spatial domains have a long history both in basic research and in applied science, see for instance [6] and references therein. In particular, for applications in geometrical classical theories in physics, as for instance general relativity, Maxwell theory, but also Ricci flow (introduced e.g. in [39]), there are interesting applications where the spatial domain is a 3-sphere. In general relativity, spatial S 3 -topology plays a particularly important role for the standard model of cosmology based on the spatially homogeneous and isotropic Friedmann-Robertson-Walker solutions (see for instance [22,40]). Beyond those simple and, at least for certain matter fields, well-understood models with high symmetry, there exist several outstanding open problems; of particular outstanding interest and motivation for this work here are the strong cosmic censorship and BKL conjecture in Gowdy vacuum solutions of Einstein's field equations for spatial S 3 -or S 1 ×S 2 -topologies [21,25,9,18,4,36].
It is not straight forward to deal with "non-trivial" spatial topologies, such as S 3 , numerically. Recall the well-known problems occurring at the coordinate axis in the case of standard cylindrical coordinates (ρ, φ, z) in R 3 . These coordinates degenerate at the "axis" given by ρ = 0. In typical equations, this has the consequence that derivatives with respect to the azimuthal angle φ always come together with a factor 1/ρ. If f is a smooth function on R 3 , then a term like 1/ρ ∂ φ f is well-behaved at the axis. However, when an equation with such terms is solved numerically, the "formal singularity" 1/ρ can cause numerical instabilities. Applications with axial symmetry have a great history in all over science. Just to name examples of numerical studies of axisymmetric problems in general relativity, we list [19,8,34]. In this paper, we will not be interested in axial symmetry, however, it is instructive to keep it in mind for the following reason. Let us assume that we cover a dense subset of S 3 with one coordinate patch, for instance the Euler coordinates introduced below. Then it turns out that this leads, loosely speaking, to two "axes" on S 3 , each of which with similar properties as the axis for cylindrical coordinates on R 3 . Hence, we call this the "axis problem" also in the case of spatial S 3 -topology.
Alternatively to this approach with one singular coordinate patch on the spatial domain, one can try a multipatch technique. The idea is to cover the spatial domain with several local regular coordinate patches. In general relativity, examples of implementations of the multipatch technique are [38,11]; of particular interest for our work here are implementations based on spectral methods in [32,35,12]. In any case, the multipatch technique does not seem advantageous for the applications which we have in mind. First, its implementation is difficult, since one must find a stable and efficient way of guaranteeing the necessary communication between the local patches. Second, we are interested in cosmological solutions with symmetries, and in order to take full advantage of those, it is often a good idea to adapt the coordinates to the symmetries, even though this can mean that one has to deal with singular coordinate systems. We have decided to develop a single patch code based on the collocation method 1 for the spatial discretization. It is general experience that such methods typically yield high accuracy [6]. Furthermore, a spectral single patch approach seems very natural from our geometric point of view, which we introduce in this paper. For the time discretization, we use the method of lines with a Runge Kutta integrator. Such a discretization technique for various spatial domains has a long history in computational physics, see for instance the references and examples in [6]. An alternative approach, which uses spectral discretization both in space and in time, has been reported on in [23]. However, to our knowledge, the case of spatial S 3 -topology has not been studied yet.
Our aim is to find numerical solutions of systems of first order quasi-linear evolution PDEs, written schematically as Here, the unknown f is a vector, the terms B are vector valued, and A i and ∇ i represent a collection of matrices and spatial derivative operators, as explained in more detail later. The spatial domain, represented by the abstract coordinates x, is S 3 . We assume in the following, without further notice, that the initial value problem for such a system is well-posed, i.e. that for any choice of initial data in a given regularity class, there exists a unique solution f locally in time, which depends continuously on the initial data in a well-defined manner. In our applications, ignoring the "axis problem" above for the moment, the operators ∇ i can be thought of as the spatial partial derivatives, the quantities A i as symmetric matrices, and all coefficients depend smoothly on their arguments. In this case, the system is symmetric hyperbolic, and well-posedness follows [26,28]. In particular, when the initial data is smooth, then the solution is smooth until it breaks down. We remark, that in principle our approach applies to other forms of hyperbolicity for Eq. (1), or even parabolic systems. We will often use the following equivalent geometric language. Namely, we say that we look for solutions of Eq. (1) on R × S 3 , where time t is interpreted as the canonical coordinate on the factor R of the manifold R × S 3 , and the t = const-hypersurfaces have S 3 -topology. In this paper, we will often be concerned with the case when f in Eq. (1) represents components of smooth tensor fields. Then we call Eq. (1) "tensorial equation". This paper is organized as follows: Section 2 is devoted to the description and discussion of our numerical technique. We show our geometric point of view which leads, in a natural way, to a single patch collocation discretization in space. We demonstrate that it allows a straight forward treatment of the "axis problem" mentioned above. Since we are interested, in particular, in tensorial equations, we discuss several related issues in Section 2.3. Namely, in order to express the tensor fields as collections of smooth functions we have decided to work with smooth global frames on S 3 . The formulation of tensorial equations in terms of smooth global frames on S 3 leads to certain problems in the presence of symmetries. Two particular classes of symmetries of our interest are discussed. For these sections, but indeed at many places of this paper, some background in differential geometry would be helpful, which can be obtained for instance from [22,29]. In Section 2.4, we elaborate on our numerical infrastructure in general, and we discuss certain issues related to numerical stability which are present for spatial S 3 -topology. In Section 3, we proceed as follows. First, we introduce the mathematical and physical background of our test application. Then we test the code in that setting and elaborate on numerical errors, stability and performance. Finally, in Section 4, we summarize and conclude. We comment briefly on possible problems when some of our special assumptions are dropped, point to open issues and list aspects for future work.
2 Geometric ideas and numerical implementation

Euler coordinates and the Euler map
In the following, we consider the 3-sphere S 3 as the submanifold of R 4 given by The Euler coordinates of S 3 , also appearing as Euler angle parametrization, etc. in the literature, is the coordinate patch given by in terms of the coordinate functions χ ∈ [0, π], λ 1 , λ 2 ∈ [0, 2π[. We will rather use the coordinates {χ, ρ 1 , ρ 2 } determined by The Euler coordinates smoothly cover the dense subset of S 3 given, when the points χ = 0, π are taken away. We expect that other choices of coordinates with similar properties are also appropriate. Although the motivation for choosing the Euler coordinates stems from Gowdy symmetry, as becomes clearer later, they are robust enough for more general cases. Certainly, the relations in Eq. (3) are well-defined for all χ, ρ 1 , ρ 2 ∈ R. This is also true when we consider χ, ρ 1 , ρ 2 ∈ (R mod 4π). Geometrically, the Euler coordinates given by Eq. (3) can thus be interpreted as a smooth map from the 3-torus to S 3 , which we call the Euler map The Euler map Φ is even a diffeomorphism when we restrict it to e.g. χ ∈]0, π[ and restrict the image correspondingly. But at the points χ = 0, π, the inverse is not welldefined. Note, that in the whole paper we will often make the canonical identification of the isomorphic groups U(1), S 1 , (R mod 2π) and (R mod 4π), and henceforth not distinguish between them. However at this point, our definition of the map requires to stand on 4π-periodicity at least for the coordinate χ, but we come to the standard 2π-periodicity in a moment. Let f be a smooth function on S 3 . In the following we considerf := f • Φ, where Φ is the Euler map defined in Eq. (4). Hencef is a smooth function on T 3 . For simplicity, we write f instead off , and often make no difference between the "original function f on S 3 " and the "corresponding functionf on T 3 ". If necessary, in order to avoid confusions, we sometimes say that "f is a smooth function on T 3 originating in a smooth function on S 3 ".
Motivated by this simple geometric relation given by the Euler map, our approach is the following one; the details are worked out in the subsequent sections. Our aim is to solve partial differential equations with a spatial domain S 3 . Let us suppose that all coefficient functions and unknown functions, which we want to solve for, in the equation are smooth functions on S 3 . Furthermore, suppose that all derivative operators stem from smooth globally defined vector fields on S 3 . Since Φ becomes a diffeomorphism, when we restrict it to a dense subset of S 3 , all these functions and vector fields correspond in a unique manner to smooth functions and vector fields on the corresponding subset of T 3 . Hence, we can solve the equation as if it its spatial domain were T 3 . However, because Φ is not a diffeomorphism globally, we have introduced formal singularities to the equations, analogous to the singularities at the axis of cylindrical coordinates on R 3 discussed above. One main conclusion in the following sections is that the hypothesis, that all these quantities on T 3 originate in smooth quantities on S 3 , allow to regularize the formally singular behavior at those points. By all this, we will successfully make "S 3 periodic in all three spatial directions", and can henceforth use spectral methods based on the standard Fourier basis for the spatial discretization of the equations.
In all of what follows, we will assume, for simplicity, that all functions involved do not depend on the coordinate ρ 2 . We call such functions U(1)-symmetric, and later we interpret this symmetry geometrically. The generalization of the following ideas is, however, straight forward. This assumption has the following nice property. Let f be a smooth U(1)-symmetric function on S 3 . Then f • Φ is 2π-periodic (instead of 4πperiodic) in χ and ρ 1 . This is so, because we can use the symmetry in ρ 2 to switch the signs of all terms in Eq. (4) consistently as needed. Hence, in the following, we only need to deal with the standard 2π-periodicity.
Note, that analogous spectral approaches for equations with spatial domains diffeomorphic to S 2 have been implemented before, see for instance [2,6] and references therein.

Analysis of Fourier series of given smooth functions on S 3
Let f be a given smooth U(1)-symmetric function on S 3 . The corresponding function f on T 3 is smooth, and hence has a representation in terms of a Fourier series. We come back to the question of convergence of such a series later. At this stage, we want to assume that this Fourier series is finite, and, without loss of generality, with N summands both in χ and ρ 1 . Since f is U(1)-symmetric function, it is 2π-periodic in both χ and ρ 1 , as argued before. Let f be real valued. Then it must be of the form Here, F 0 (χ) is a real valued function, and F p (χ) is allowed to be complex for all p ≥ 1. In [3], we demonstrated that the properties of the Euler map under U(1)-symmetry implies for some, in general, complex coefficients f n,p ; only the coefficients for p = 0 must be real. Some of the factors in these expressions are chosen for later convenience. Note, that in [3], the coordinates (χ, ρ 1 , ρ 2 ) are defined slightly differently. For example, the function χ here must be substituted by 2χ to compare to the expressions in [3]. In any case, we further find that for even p > 0, there are the "compatibility conditions" These originate in the fact that the functions F p (χ) must vanish at the degenerated places χ = 0, π for all p > 0. For odd p, this is implied by the expression Eq. (6a) automatically, and hence there are no compatibility conditions. The corresponding expressions are more complicated, but analogous, when we give up U(1)-symmetry. Now, we want to study what happens to the Fourier series of f , when it is differentiated along a smooth tangent vector field on S 3 . By this, we mean that the abstract derivative operator Here, our convention is that x α represents three abstract spatial coordinates, which neither need be the Cartesian coordinates, nor the Euler coordinates used before. However, we will restrict to Euler coordinates in the following. Furthermore, we assume Einstein's summation convention. The coefficients V α are just functions on S 3 which are not necessarily smooth when we consider the Euler coordinates. Of particular importance will be the following tangent vector fields, whose origin we explain later and which, with respect to the Euler coordinates, take the form The factors 2 are chosen for consistency with our discussion in [3]. Now, as we explain later, any smooth vector field V on S 3 can be written as a linear combination V = V a Y a with V 1 , V 2 , V 3 smooth functions on S 3 . Hence, under the assumption that all differential operators in the equations stem from smooth vector fields on S 3 , it follows for U(1)-symmetry, that all "formally singular" differential operators in our equations must be of the form −F (χ, ρ 1 ) cot χ∂ ρ 1 , with F some smooth function on S 3 . Without the assumption of U(1)-symmetry, there can additionally be singular operators of the form (F (χ, ρ 1 , ρ 2 )/ sin χ)∂ ρ 2 . This shows that the formally singular terms here are of the same type as in the "axis problem". The differences to the case of cylindrical coordinates on R 3 are twofold. First, in the case of S 3 , we have two such "axes" simultaneously at χ = 0 and π. Second, the axis itself is not topologically a line here, but a closed circle. In the following, we restrict our attention to the operator relevant to U(1)-symmetry − cot χ∂ ρ 1 .
So, let f be as before, and g := ∂ ρ 1 f . Since g is again a smooth U(1)-symmetric function on S 3 with finite Fourier series, the analogue of Eqs. (5) and (6) holds with F p (χ) substituted schematically by the function G p (χ). Then, for χ = 0, π, we get, for p > 0 odd. The computations leading to this result are described in [3]. This result means that, as soon as the Fourier coefficients of ∂ ρ 1 f are known, the Fourier coefficients of the complete "formally singular term" − cot χ ∂ ρ 1 f can be computed. Now, let us consider the general case of a given smooth U(1)-symmetric function f on S 3 . The associated function on T 3 has an infinite Fourier representation in a general, with rapidly decreasing coefficients f n,p . This is a standard result from Fourier analysis, which can be found in [37,7]. This last property means that the modules of the Fourier coefficients is bounded by a uniform constant times any negative integer power of the two summation indices n and p. This property is often referred to as "exponential convergence". It is a general fact under our conditions that the Fourier series converges pointwise absolutely and even uniformly. We find straight forwardly that the Fourier series of f must be of the form given by Eqs. (6), setting N → ∞. The infinite series of the compatibility condition Eq. (6b) converges because the coefficients are rapidly decreasing. Now consider the inverse question. Let a function be given on T 3 , of the form of Eqs. (5) and (6), for N = ∞ with rapidly decreasing coefficients f n,p . The standard theory implies that the series converges pointwise absolutely and uniformly to a smooth function f on T 3 . However, does f originate in a smooth function on S 3 ? In general the answer is no, because the compatibility conditions Eq. (6b) are just necessary, but not sufficient for smoothness. Indeed, it is sufficient that for any p, the function F p (χ) is a smooth 2π-periodic function on R, which has a zero of order p at χ = 0 and π. In particular, each function F p (χ)e ipρ 1 then originates in a smooth function on S 3 . The argument for proving that this implies that f is a smooth function on S 3 , uses the theorem about Fourier series on S 3 in [37]. Namely one can show that under these assumptions, f as a function on S 3 can be represented as an infinite series of spin-weighted spherical harmonics with again rapidly decreasing coefficients.
Consider the derivative g = ∂ ρ 1 f . In particular, the formula for − cot χ g in Eqs. (8) also holds in the limit N → ∞, and the series expressions there converge at least pointwise at all χ = 0, π. The function − cot χ ∂ ρ 1 f is a smooth function on T 3 because each F p (χ) in Eq. (5) is a smooth 2π-periodic function in χ with appropriate zeros at the "singular places" χ = 0, π. This is a nice result because it shows that Eq. (8) is meaningful at the singular locations, and hence allows to evaluate the formally singular term − cot χ ∂ ρ 1 f explicitly there, even in the limit N → ∞. However, we remark that − cot χ ∂ ρ 1 f does not originate in a smooth function on S 3 , because Eqs. (8) is not consistent with Eq. (6). This is not a problem because the formally singular operator is only a part of a differential operator defined by a smooth vector field on S 3 . Indeed, the result, when this "full" differential operator is applied to a smooth function, yields a smooth function on S 3 .

Symmetry and related issues for tensorial equations
Before we discuss, how these results can be applied in practice, let us first consider some consequences for tensorial equations.

Smooth frames on S 3
Recall, that one of our main assumptions is that, at any given instance of time, all unknowns and coefficients in the equations are smooth functions, and that all differential operators are determined by smooth globally defined vector fields on S 3 . However, in order to turn tensorial equations into partial differential equations for smooth scalar functions, we need to introduce smooth frames on S 3 . We would like to mention that an alternative way of treating tensorial equations on S 3 in the case of Gowdy symmetry, see below, can be found in [18] for the case of spatial S 1 × S 2 -topology.
Let us recall some well-known facts. Let S 3 be given as in Eq. (2). Assume the standard matrix representation of the Lie group SU(2) [37]. The map is a diffeomorphism, which can be used to transport the group structure of SU(2) to S 3 . Hence, both SU(2) and S 3 can be considered as identical Lie groups via the map Ψ. Thus, from the standard SU(2) group multiplication, we can define left and right translation maps, so that L u and R u are diffeomorphisms S 3 → S 3 for each point u ∈ S 3 . Those maps can be employed to construct smooth global frames. First one chooses a basis of the tangent space at the unit element e of the group. We choose the Pauli matrices with non-standard normalizatioñ considered as elements of T e (SU (2)). Then, one uses the push forward of L u or R u to transport this basis smoothly to any other point u ∈ SU(2) (2)-right invariant and both are smooth global frame fields on S 3 . It is straight forward to check that they satisfy where ǫ abc is the totally antisymmetric symbol with ǫ 123 = 1. Here, the brackets denote the Lie bracket. For {Y a }, we have already written the explicit expressions with respect to the Euler coordinates in Eqs. (7). For {Z a }, we have On R × S 3 with a time function t, we assume that each t = const-hypersurface is diffeomorphic to S 3 with Euler coordinates {χ, ρ 1 , ρ 2 }. Hence, on each of these surfaces, the expressions Eqs. (7) and (10) define the fields {Y a } and {Z a }. Geometrically, we thus have for all a = 1, 2, 3 Now we write an arbitrary globally defined smooth frame {e i } on R × S 3 as follows. By {e i } we mean the collection of 4 frame fields {e 0 , e 1 , e 2 , e 3 }. We set and write, where (e b a ) is a smooth 3 × 3-matrix valued function with non-vanishing determinant on S 3 . Our conventions for frames is that the index 0 always corresponds to the "time frame vector" e 0 and a, b, . . . = 1, 2, 3 correspond to the "spatial frame vectors". When we write indices i, j, . . . = 0, . . . , 3, we mean both the time and spatial frame vectors. For a tensor field S, say for example a covariant 2-tensor, we write S ij := S(e i , e j ), S ab := S(e a , e b ) etc. We stress, that it is important to understand, that writing e i , does not mean the ith component of a vector e, but rather the ith vector field of the frame {e i }.

Symmetry reductions of tensorial equations
Let S be an arbitrary smooth tensor field on R×S 3 . We say that S is ξ-invariant, provided L ξ S = 0 everywhere. Here L ξ denotes the Lie derivative along ξ. The coefficients S ij of S, with respect to an arbitrary frame {e i }, are constant along ξ, if, and in general only if, [ξ, e i ] = 0. Now, suppose that the functions S ij of such a tensor are among the unknowns of the system of partial differential equations which we would like to solve. Often, we would like to exploit the symmetry of the unknowns and reduce the equations to some simpler form. If the unknown functions are constant along ξ, such a reduction can be done directly, if ξ has the meaning of a spatial coordinate vector field. Hence, in the following, we will be interested in frames such that [ξ, e i ] = 0; in this case, we say that the frame is ξ-invariant.
Let us consider special cases of interest for us. We have already introduced the notion of U(1)-symmetry for functions before. For general smooth tensor fields S on R × S 3 , we define it by the requirement that S is Z 3 -invariant. One can define U(1)-symmetry more geometrically, but for the purpose of this paper, our definition is sufficient. The integral curves of Z 3 = ∂ ρ 2 are circles. Hence the symmetry group is isomorphic to U(1), which motivates the name. Now, it is straight forward to construct Z 3 -invariant frames {e i }. Namely, for the ansatz Eqs. (11), this requirement is equivalent to ∂ ρ 2 (e b a ) = 0, due to Eqs. (9). The consequence is, that both the frame matrix (e b a ) and the frame components of all tensor fields, which are Z 3 -invariant, are constant along ρ 2 , and hence are U(1)-symmetric functions. Provided, our equations are formulated with respect to functions with this property only, the spatial domain reduces to that of the coordinates {χ, ρ 1 }. We refer to this as the 2 + 1-reduced equations.
Another symmetry assumption of interest is Gowdy symmetry with two spatial symmetry vector fields. We say, that a tensor field is Gowdy symmetric, if it is U(1)symmetric and additionally Y 3 -invariant. Again, this definition can be made more geometric. In any case, let us suppose that a frame field {e i }, obeying Eqs. (11) and [Z 3 , e i ] = 0, is given as before. Now, it turns out, as argued in [3], that the assumption that {e i } is a smooth globally defined Z 3 -invariant frame on S 3 , does not allow that it is Y 3 -invariant in addition. Hence, although we can find a frame such that the frame components of arbitrary Gowdy symmetric tensor fields are constant along ρ 2 , it is not possible to achieve that in addition they are constant along ρ 1 ; recall Y 3 = ∂ ρ 1 . Thus, even if we assumed Gowdy symmetry, the spatial domain of our equations would not reduce directly to that of the coordinate χ, i.e. not to 1 + 1 dimensions. This difficulty is a consequence of our assumption that the frame is smooth globally on the manifold; if we allowed the frame to become singular at some places, then the situation would be different. However, we would get other problems due to the additional singularities.
Nevertheless, even under these assumptions, it is possible to perform the following "indirect" reduction of the equations to 1+1 dimensions in the case of Gowdy symmetry. For the frame components of any smooth tensor field S, for example in the case of a covariant 2-tensor, which is Y 3 -invariant, we find with . Under our assumptions for the frame Eqs. (11), we see that For the spatial components, one has The matrix (f a c ) is defined here as the inverse of the matrix (e c a ). As soon as we fix the transport of the frame in time, we can compute the time derivative of the matrix (T a ′ a ); we do this for a particular example in Section 3.1. In general, we can expect that these evolution equations for the matrix (T a ′ a ) are non-trivial. In any case, the idea for the "indirect" reduction to 1 + 1 in the Gowdy case is the following: first, substitute the Y 3 -derivatives of all tensor field components in the equations by means of Eq. (12). Second, append the evolution equations for the matrix (T a ′ a ) to the system of equations and, third, evaluate the equations only at ρ 1 = 0. Then, with respect to the Euler coordinates, all unknowns only depend on t and χ, and the evolution system is closed. Note, however, that it depends on the properties of the evolution equations of (T a ′ a ), whether the resulting system of evolution equations yields a well-posed initial value problem. In the example, which we discuss later, this is the case.
We remark that under the assumption of U(1)-symmetry, all results obtained here also hold for spatial S 1 × S 2 -topology. We do not elaborate on this further; a discussion is given in [5].

Discretization and our numerical infrastructure
In order to compute approximate solutions of our system of partial differential equations Eq. (1) by means of a computer, we need to discretize the equations and the unknowns. Our analysis before based on Fourier series suggests spectral discretization [7,6] in space with the standard trigonometric basis. We follow most of the conventions in [6]. In order to keep the presentation as short as possible here, we do not write down formulas wherever we follow the standard conventions. In particular, we use the collocation method. For the spatial grid in any of the spatial dimensions, referred to as x, we set where N is the number of grid points in the chosen spatial direction. The quantity µ ∈ [0, 1[ is a shift quantity. For the standard collocation method which we use, the quantity N must be odd. For simplicity we use the so-called partial summation algorithm [6] for computing the discrete Fourier transforms (DFT) so far; however, we plan to switch to the Fast Fourier Transform algorithm (FFT) [10,33], in order to optimize performance for high spatial resolutions.
This discretization of the equations and unknowns in space yields a system of ordinary differential equation (ODE) in time for the spectral coefficients of the unknowns, or equivalently for the values of all unknown functions at the spatial grid points. This system of ODEs is called the semi-discrete system. In order to solve this numerically, one must discretize time as well. For this, we have implemented a couple of Runge Kutta (RK) variations described in [33]; namely, first, the non-adaptive 4th order RK scheme, second, the 4th order "double-step-adaption" RK scheme and, third, the adaptive 5th order "embedded" RK scheme. For those schemes, the time adaption is always "global in space", namely, at a given time the maximal estimated error at all spatial points is taken. The parameter η is the desired accuracy, according to Eq. (16.2.7) in [33], where it is called ∆ 0 . The lower its value, the stronger is the tendency of the adaption scheme to decrease the time step h. For practical reasons, we also define a minimal time step h min , so that the adaption scheme is prevented from reaching unpractically small values of h.
A sophisticated discussion of errors and convergence in such discretization approaches is given in [7]. We will not elaborate on this, in general, very complicated problem analytically. Instead, we will investigate errors and convergence in our test applications empirically in Section 3.2. The positive experience, which people have gathered over many years of research with the collocation method, is summarized in Boyd's empirical "assumption of equal errors" [6], which we decided to rely on in our numerical work.
It is clear that many classes of problems require adaptive techniques for the spatial resolution. One particular effect for underresolved numerical solutions obtained by the collocation method is aliasing. We have not yet implemented any of the explicit antialiasing recipes, given for instance in [6]. Instead, we use the following simple global spatial adaption technique so far. At each time step of the numerical evolution, the program computes the Fourier transform of one representative unknown function; which one is chosen requires some experiments. In our applications, where symmetry implies that only one spatial direction is significant, it is then sufficient to do the following. From the spectral coefficients of this unknown, the code determines the "power" of the upper third of the frequency spectrum with respect to the significant spatial direction, divided by the total power. It is straight forward to generalized this to more general situations. By "power" we mean the sum over the squares of the modulus of the Fourier coefficients. We call the result of this computation the "adaption norm" Norm (adapt) . Besides adaption itself, this "norm" can also be thought of as a measure for the aliasing error. When, during the numerical evolution, Norm (adapt) exceeds a prescribed threshold value, the code stops automatically, interpolates all quantities to a higher spatial resolution and continues the run. In each of these adaption steps, we have found it to be reasonable to almost double the spatial resolution. In any case, note that this is a primitive adaption method, since it is "global in space". In particular, for solutions, which develop sharp localized features, a local adaption method in space would be more desirable. This is a future work project. However, even Gowdy spikes have been treated with our much simpler method in [3].
Let us, furthermore, mention the possibility of the Intel Fortran compiler [24] on Intel CPUs, which we have worked on exclusively up to now, to switch from the standard machine supported number representation called "double precision" with round-off errors of the order 10 −16 to software emulated "quad precision" with round-off errors of the order 10 −32 . Indeed, this possibility is exploited in our applications, as is discussed later.

Practical issues for spatial S 3 -topology
In this paper, we discuss two codes for spatial S 3 -topology which both use the infrastructure above, but with slight differences. For one of the codes, we assume U(1)-symmetry with a choice of orthonormal frame reducing the equations to two spatial dimensions; we call this code the 2 + 1-code. For the other code, which we call 1 + 1-code, we assume Gowdy symmetry, and suppose that the indirect reduction to one spatial dimension described before leads to a well-posed initial value formulation. The particular equations which we implement and study in both cases are discussed in Section 3.1.
In order to apply the results, which we obtained for the properties of Fourier series of smooth U(1)-symmetric functions on S 3 in Section 2.2 in our discretization approach, we assume that for any choice of resolution, the solution of the corresponding discretized equations originates in smooth functions on S 3 . Let us suppose, that our evolution equations have the property, that if the initial data of the continuum problem is smooth, then the corresponding solution of the continuum equations is smooth. For instance, this is the case for symmetric hyperbolic systems. Then, if the initial data is approximated by functions, which originate in smooth functions on S 3 , and, if all our assumptions about the coefficients and derivative operators in the equations before hold, then the solution of the discretized originates in smooth functions on S 3 for any resolution. This is at least true, as far as we can neglect errors caused by aliasing and the finite number representation in our computer. Let us assume just for a moment, that these errors can be neglected. In particular, Eq. (8) can be applied to the 2 + 1-code directly. For the 1 + 1-code, we need one further observation, since the equations are only evaluated at ρ 1 = 0, and hence there is no information about even and odd Fourier modes with respect to ρ 1 . Namely, due to Eq. (6a), all coefficients associated with cos-modes with respect to χ must correspond to an even mode with respect to ρ 1 . Analogously, all coefficients associated with sin-modes with respect to χ must correspond to an odd mode with respect to ρ 1 . This information is sufficient to use Eqs. (8) as in the 2 + 1-case. Now, for even p > 0, there are two ways of computing the coefficients b N r,p and c N r,p of the formally singular terms. Namely, due to the compatibility conditions Eqs. (6b), we can write both For odd p, there is only one way of writing these coefficients. Although each pair is equivalent in exact computations, there can be a difference numerically. We refer to the first way of computing these coefficients as "up-to-down", since we need the information of all high frequencies to compute the low frequency coefficients recursively. The second variant is called "down-to-up", since the information from low frequency coefficients is used to compute the high frequency coefficients recursively. A priori, both ways have the potential to amplify numerical instabilities. In particular, although the solution originates in smooth functions on S 3 at one time of the evolution, this together with round-off errors and aliasing can cause a drift, so that the form given by Eqs. (6) are violated eventually. We have not yet built the special structure of the Fourier series into our numerical infrastructure. Thus, it is possible, that such errors accumulate, such that, after some time of evolution, the numerical solution does not represent a smooth solution on S 3 anymore. Indeed, we found in our numerical experiments with the 2 + 1-code in [3], that, without precautions, the numerical solution typically drifts away strongly for both the up-to-down and down-to-up method. We were not able to pin-point the problem. However, for the down-to-up method, it turns out to be sufficient, after each time step, to set all those Fourier coefficients to zero explicitly, which are supposed to vanish according to Eq. (6a). With this manipulation, the numerical evolution becomes stable. In particular, Eq. (6b) stays satisfied within reasonable error limits, and the code is convergent and able to reproduce explicitly known solutions. The up-to-down method, however, we were not able to stabilize.
We can expect that similar practical issues exist for the 1 + 1-code. Here, however, we must proceed slightly differently, since the form given by Eqs. (6) cannot be enforced explicitly. In order to control the smoothness of the numerical solution nevertheless, the idea is to control the unknowns directly at the coordinate singularities χ = 0, π in terms of "boundary conditions" 2 . Let the frame {e i } be Z 3 -invariant as before, and let S be one of the unknown Z 3 -invariant tensor fields. Since Y 3 = ±Z 3 on the symmetry axes, one obtains Y 3 (S ij ) = 0 there. Exploiting this information by means of Eq. (12), implies a homogeneous linear algebraic "boundary system", which yields the "boundary conditions" for S For our particular test case later, we solve the boundary system in Section 3.1.3. Let us assume for the moment that we have solved this system. In general, we would like to have the possibility of either letting the numerical evolution proceed freely and just monitoring, how well those boundary conditions are satisfied, or, if necessary, we would like to enforce the boundary conditions. The latter means that we set the values of the unknowns at χ = 0, π explicitly to the values implied by the boundary conditions. In order to make this possible, we modify the spectral conventions slightly as follows, so that both boundary points χ = 0 and χ = π correspond to grid points. Let f be some unknown function, which we discretize as (a n cos nχ + b n sin nχ).
In the standard collocation approach, which we use for the 2 + 1-code, the number of grid points N , according to Eq. (13), is odd, and M = N . Recall that the discrete Fourier transform [6] is the linear map from the values of f at the grid points (f (x 0 ), . . . , f (x N −1 )) to the Fourier coefficients (a 0 , b 1 , a 1 , . . . , b (M −1)/2 , a (M −1)/2 ), which is bijective for these choices of N and M . This is true for any choice of shift µ. For the 1 + 1-code now, we choose µ = 0, any even number N , and set M = N + 1. In this case x 0 = 0 and x N/2 = π. One finds easily that for this, the standard discrete Fourier transform is the map Hence, the map has the standard properties except for the highest frequency. The fact, that the discrete Fourier transform always yields zero for the highest sin-mode can be understood easily, because the value of sin nχ is always zero for n = (M − 1)/2 at χ = k 2π N . The main point is now that this map is nevertheless invertible. For spectral differentiation, we ignore the frequency n = (M + 1)/2 completely. In practice, it is expected that this is not problematic, since the highest frequencies are typically insignificant. Now, in our numerical experiments, which have so far restricted to the down-to-up method, we find that the 1 + 1-code with this spectral infrastructure is very stable and convergent. This is true even without enforcing the boundary conditions at all and hence without any explicit control of the smoothness in regimes where the solution is relatively smooth. Recall that this is not so for the 2 + 1-code. Furthermore, the violations in the boundary conditions typically converge to zero with increasing resolution. However, if the simulation approaches a non-smooth regime of the solution, it seems often necessary to enforce the boundary conditions; this is discussed for our test applications.

Physical and mathematical background
Our aim is to compute cosmological solutions of Einstein's theory of relativity; in particular we are interested in the strong cosmic censorship conjecture in Gowdy vacuum solutions of Einstein's field equations for spatial S 3 -or S 1 ×S 2 -topologies [21,25,9,18,4,36]. All our discussions assume vacuum and a cosmological constant λ, so that Einstein's field equations (EFE) in geometric units c = 1, G = 1/(8π) read Hereg µν is the spacetime metric, which is the fundamental unknown encoding the information about the gravitational field. Its Ricci tensorR [22] is a 2nd order quasi-linear expression in the metric. We will always assume four spacetime dimensions, that the signature of the metric is Lorentzian (−, +, +, +), and that Cauchy surfaces, i.e. the "surfaces of constant time", are diffeomorphic to S 3 . Furthermore, we suppose λ > 0. In [30,31], Penrose introduced his notion of conformal compactifications. The idea is to rescale the physical metricg by means of a conformal factor Ω, which is a smooth strictly positive function on the spacetime manifoldM . This yields the so called conformal metric g := Ω 2g .
Now, loosely speaking, if it is possible to attach those points toM , which are the limit points of vanishing Ω, so that the new manifold M is smooth and the metric g can be extended as a smooth metric on M , then we say that the original spacetime has a smooth conformal compactification. The references above, but in particular [15], give further necessary technical requirements to make this loose statement rigorous. Under those conditions, the set Ω = 0 is a smooth surface in M , called conformal boundary J . Physically it represent "infinity". In [15], it is shown, that conformal boundaries must be spacelike hypersurfaces with respect to the conformal metric for all solutions of Eq. (14) with λ > 0. One calls such solutions "future asymptotically de-Sitter" (FAdS) [13,1], if its conformal boundary has a smooth non-empty future component J + ; there is the analogous concept for the past time direction. In particular, the de-Sitter solution [22] is FAdS. Under these conditions, J + represents the infinite timelike future ofM . Some of the asymptotic geometric properties of FAdS solutions are discussed in [3]. Friedrich introduced his conformal field equations (CFE), as reviewed for instance in [15], in order to deal with conformally compactified solutions of Einstein's field equations. In these conformally invariant equations, the fundamental unknown is the conformal metric g and the conformal factor Ω related to the physical metricg. The non-trivial property of these equations is, that they are, first, equivalent to Einstein's field equations wherever Ω > 0, and, second, yield regular hyperbolic evolution equations even where Ω = 0. Under the assumptions above, the CFE allow us to formulate what we call "J + -Cauchy problem" [13]. The idea is to prescribe data for the CFE on the hypersurface J + , including its manifold structure, subject to certain constraints implied by the CFE. These data can then be integrated into the past by means of evolution equations implied by the CFE. Friedrich proved that the J + -Cauchy problem is well-posed, and that the unique FAdS solution corresponding to a given choice of smooth data on J + is smooth, as long as it can be extended into the past. It is remarkable that this Cauchy problem allows to control the future asymptotics of the solutions explicitly by the choice of the data on J + . Concerning the past behavior of the solution corresponding to a given choice of data on J + , however, there is only limited understanding and a-priori control, because of the complexity of the field equations. In this paper, we will give no details on the constraints on J + , and say only briefly what the relevant initial data components are, since we do not want to introduce all necessary geometric concepts now. However, we write down a special class of solutions of the constraints in Section 3.1.2. We refer to [13,3], where the details have been carried out.
We decided to use the so-called general conformal field equations, which are the CFE in conformal Gauss gauge [14,15]. In our applications, we specialize the gauge even further to what we call Levi-Civita conformal Gauss gauge [3]. In this paper here, we will discuss neither the physical properties, nor the possibly bad implications of this choice of gauge, but just refer to [3,4]. In any case, assuming, without loss of generality, λ = 3, and having fixed the residual gauge initial data, as described in [3], the implied set of evolution equations is for the unknowns The unknowns are the spatial components e b a of a smooth frame field {e i } as in Eq. (11b), with e 0 = ∂ t , which is orthonormal with respect to the conformal metric, the spatial frame components of the second fundamental form χ ab of the t = const-hypersurfaces with respect to e 0 , the spatial connection coefficients Γ b a c , given by Γ b a c e b = ∇ ea e c −χ ac e 0 where ∇ is the Levi-Civita covariant derivative operator of the conformal metric, the spatial frame components of the Schouton tensor L ab , which is related to the Ricci tensor of the conformal metric by and the spatial frame components of the electric and magnetic parts of the rescaled conformal Weyl tensor E ab and B ab [15,16], defined with respect to e 0 . In this special conformal Gauss gauge, the timelike frame field e 0 is hypersurface orthogonal, i.e. (χ ab ) is a symmetric matrix. In order to avoid confusions, we point out that, in principle, the conformal factor Ω is part of the unknowns in Friedrich's formulation of the CFE. However, for vacuum with arbitrary λ, it is possible to integrate its evolution equation in any conformal Gauss gauge explicitly [14], so that Ω takes the explicit form Eq. (15g) in our gauge. We note, furthermore, that, since (E ab ) and (B ab ) are tracefree by definition, we can get rid of one of the components for each of the two. Our simple minded choice is the 33-component by E 33 = −E 11 − E 22 ; the same for the magnetic part. The evolution equations Eqs. (15e) and (15f) of E ab and B ab are derived from the Bianchi system [15]. In our gauge, the constraint equations implied by the Bianchi system take the form Here, ǫ abc is the totally antisymmetric symbol with ǫ 123 = 1, and indices are shifted by means of the conformal metric. The other constraints of the full system above are equally important, but are ignored for the presentation here. Further discussions of that evolution system and the quantities involved can be found in the references above. Note that in Eqs. (15e), (15f) and (16), the fields {e a } are henceforth considered as differential operators, using Eq. (11b) and writing the fields {Y a } as differential operators in the Euler coordinate basis as in Eqs. (7). Seen as a system of partial differential equations, the system Eqs. (15) is symmetric hyperbolic, and hence the initial value problem is well-posed. Note that in this gauge, our initial hypersurface J + corresponds to t = 0. The past conformal boundary, if it exists, corresponds to t = 2. Hence, our time coordinate runs backwards with respect to physical time.
These equations hold without any symmetry assumptions. In the following we will assume that all unknown fields involved are Gowdy symmetric. For the 1 + 1-code, we need to derive the evolution equations of the matrix (T a ′ a ). The fact that the conformal metric g is Y 3 -invariant, implies, after straight forward computations, that the matrix (T ab ) := (T a ′ a g a ′ b ) is antisymmetric. Our choice of frame transport is parallel transport with respect to the conformal metric. This implies, after some algebra, that Hence, since the original system of equations is symmetric hyperbolic, also the "indirect reduction" to 1 + 1-dimensions is symmetric hyperbolic. So, the initial value problem for these equations is also well-posed.

A class of initial data
As initial data on J + , we use the "Berger data", which are solutions of the constraints derived for U(1)-and Gowdy symmetry in [3]. Those data are close to data of the λ-Taub-NUT solutions and hence are particularly interesting for the strong cosmic censorship conjecture [4]. Here, we restrict to Gowdy symmetry. Under the conventions above, these data take the form (e b a ) = diag(1, 1, a 3 ), (17a) The induced conformal 3-metric of J + is a Berger sphere with a free parameter a 3 > 0. The only inhomogeneous, i.e. space dependent part of the initial data is given by the components E ab . For our definition of the functions w np , consult [3]; we just note that, with respect to the Euler coordinates, we have For all these data, one finds In total this family of solutions of the constraints has three free parameters a 3 > 0, E 0 ∈ R and C 2 ∈ R. We remark that the reason for the strange names of these parameters is consistency with our notation in [3].

Boundary control for the 1 + 1-code
In Section 2.4.2, we have motivated our boundary control approach for the 1 + 1-code. Because the analysis depends strongly on the particular equations and choice of frame transport, it was not possible to give a further discussion there in full generality. Hence, let us continue here for our special choice of equations and frame transport. Due to what was said before, we have for all t, for our choice of initial data. In this case, we say, that the frame is "boundary adapted". Now, we introduce the new fields and similar for the magnetic part B ab , so that the boundary system, introduced in Section 2.4.2, yields the following conditions at χ = 0 and χ = π, whereas E * 11 and B * 11 are free. For all other symmetric invariant 2-tensor fields, for instance the 2nd fundamental form, we get analogous relations, but, in addition, the 3-3-components are free, if the tensor is not tracefree. Since the behavior of the connection coefficients Γ b a c can be derived on the symmetry axes as well, which are the only nontensorial objects in our set of unknowns, we obtain a complete set of boundary conditions for all the unknowns. The quantity Norm (BC) is now defined as the sum of the actual absolute numerical boundary values of all those unknowns, which are supposed to be zero there according to these results. Monitoring Norm (BC) in a numerical computation, yields information on how well the boundary conditions are satisfied.
In order to implement the 1 + 1-code numerically, we write the unknowns in terms of the new electric and magnetic fields defined in Eq. (18). Actually, it would be better to introduce the analogous combinations of fields for the other unknowns, but this has not yet been done. Hence, so far, the code lacks a clean way of enforcing e.g. the boundary condition χ 11 − χ 22 = 0 at χ = 0 and π. To circumvent this problem temporarily, we have decided to work with a "partial enforcement" scheme, which, at a given time of the evolution, enforces all boundary conditions except for those of this type. In addition, we monitor the quantity Norm (BC) , and so far this treatment has turned out to be sufficient.

Numerical results for the test application
In [3], we have performed a couple of tests with the 2 + 1-code, discussed the findings and drew conclusions about the numerical method. Here, we rather focus on the 1 + 1code, and show so far unpublished tests and discussions in Section 3.2.1. Afterwards in Section 3.2.2, we also compare a simulation done with the 2 + 1-code and the 1 + 1-code directly.
We just note that we have not made systematic investigations of the CFL condition for our codes yet.

Analysis of computations with the 1 + 1-code
For our numerical test case here, we choose a 3 = 0.7, C 2 = 0.1 and E 0 = 0 in Eqs. (17) as initial data parameters, corresponding to the "large inhomogeneity case" in [3] and to one of the simulations presented in [4]. The associated solution turns out to develop a singularity, and hence can be seen as an interesting test case for our code. The evolution of a spatial norm of the curvature invariant called Kretschmann scalar is shown in Fig. 1. All the results we show here were done without the automatic spatial adaption approach described in Section 2.4.1, because, in order to study convergence, it seems more useful to control and adapt the spatial resolution manually. The adaption norm, computed with respect to E 13 , was used only for estimating the aliasing error. The time integration was done with the adaptive 5th order embedded RK scheme with control parameters η , where the Ricci tensorR ij of the physical metricg ij is evaluated algebraically from the conformal Schouton tensor L ij and derivatives of the conformal factor Ω. The indices in this expression are defined with respect to the physical orthonormal frame given bỹ e i = Ωe i . The norm is computed by summing over the L 1 -norms of each component. Now, we will distinguish two phases of the evolution for these initial data, in which different aspects and effects are important: the early evolution close to J + for t between 0.0 and 0.69, and the late evolution close to the singularity, which we find at t ≈ 0.69520493. In order to avoid confusions, we recall that the terms "early" and "late" are always understood with respect to the time coordinate t, which, however, runs backwards with respect to the physical time.
At early times, it is achieved easily that the spatial discretization error is not significant, until some later time when small spatial structure starts to form more rapidly. One hint that this is true, as we do not show here, is that Norm (adapt) is more or less constant over a long time period, and small. Another hint is that the behaviors of Norm (constr) and Norm (einstein) are not strongly influenced by the spatial resolution. See Fig. 2 and where N represents the number of spatial grid points, which is constant in this early regime. Indeed, the higher N , the larger is the initial value Norm (constr) due to higher round-off errors for computing spatial derivatives. This is not visible for Norm (einstein) , since this quantity is defined purely algebraically in the unknowns. In Fig. 2, we see that Norm (constr) grows less, the higher the time resolution is, i.e., in particular, the smaller the parameters η and eventually also h min are. However, we always observe at least weak approximately exponential growth. In Fig. 3, we see a similar behavior for Norm (einstein) . We do not show here that there is neither a particular growth of Norm (constr) , nor of Norm (einstein) , at the symmetry axes. Rather, the maximal growth takes place, where the curvature increases most strongly. This can be seen as a confirmation that our treatment of the coordinate singularities works well, cf. [3]. Note that there is an optimal time resolution, in the sense that, if we choose a higher resolution, the constraint error and Norm (einstein) are actually increased caused by higher round-off errors. Although we do not show any plots, we want to mention, that we have experimented with "quad precision". For the 1 + 1-code, this yields reasonable performance and has several consequences. First, the initial data for the constraint violations are decreased by many orders of magnitude, since those are determined primarily by the precision of the numerical number representation. By choosing appropriate resolutions, we find that the constraint violations and Norm (einstein) can then be kept several orders of magnitude smaller than in the double precision case during the whole run. However, they always show exponential growth, which suggests that this is the typical behavior of the constraint propagation in our system of equations. Furthermore, quad precision allows us to work in a regime in which discretization errors are much larger than roundoff errors, and hence it is easier to interpret convergence tests. Fig. 4 shows the behavior of Norm (BC) . The errors at the boundaries are small and stable, despite of some weak growth, which is expected in situations close to a singularity. In these tests, this error does not decrease for higher resolutions, which is a hint that round-off errors have a significant effect. With quad precision, we were able to confirm that the violations of In the runs underlying these plots, we adapted the spatial resolution several times during the runs manually; the number N in the figures is the final spatial resolution in each case. Each manual spatial adaption step is visible in the plots as a jump, because for different spatial resolutions, the numerical values of the norms slightly change. In all these runs, we choose η = 10 −13 , and hence the time steps h decrease so strongly that h = h min at that time when the runs were stopped. In this late time regime, the errors are dominated by spatial discretization errors, because the solution has the property that spatial structures shrink without bound. This can be seen by looking at the late time plot of the adaption norm in Fig. 5. It shows, how strongly the demand for spatial resolution grows with time, but also, that it is possible to gain control by increasing the resolution at least temporarily. However, the demand for spatial resolution increases very strongly with time, and it turns into a difficult numerical issue to keep track of that eventually. In Fig. 6, we demonstrate, how the choice of spatial resolution influences the propagation of the constraint violations, and that this quantity converges to a weakly exponentially growing for sufficiently high spatial resolutions. This is a promising result, and shows that the constraint propagation is more or less under control, as long as it is possible to increase the resolution in practice. Fig. 7 shows the violations of the boundary conditions. They turn out to be very small and under control for sufficiently large resolutions. The higher the final spatial resolution is, the smaller these violations appear. Finally, in Fig. 8 we show Norm (einstein) . We do not observe a very strong difference between the various resolutions; indeed, in order to make the differences visible at all, the time axis represents an even smaller time neighborhood now. It is unexpected that at very late times this norm is not necessarily smaller the higher the resolution is, and it has to be investigated whether this is a problem. Next, we want to compare directly the numerical results obtained with the 2 + 1-code and the 1 + 1-code for a Gowdy symmetric solution of our equations. We restrict to a singularity free case, since we do not want to be spoiled by the lack of spatial resolution, in particular for the 2+1-code. The following initial data parameters are chosen: a 3 = 0.93, E 0 = 0 and C 2 = 0.5, which correspond to the "regular case" in [3]. The corresponding solution is smooth for all 0 ≤ t ≤ 2 and hence develops a smooth past conformal boundary at t = 2. Fig. 9 shows the evolution of the curvature invariant Kretschmann scalar. We show neither convergence plots, nor the error quantities before, since the situation is very similar to the early phase of the singular solution in the previous section. Instead, we report only on one run done with one fixed resolution: the size of the time step is h = 5 · 10 −4 and the number of spatial points is N = 40 for the 1 + 1-code and N χ = 41, N ρ 1 = 21 for the 2 + 1-code. We use the non-adaptive 4th order RK time integrator, and the automatic spatial adaption has been switched off here as well. The run with the 1 + 1-code was done without enforcing the boundary conditions now, cf. Fig. 10. We see that this error stays very small and is stable. In order to illustrate, how well the results obtained with the 1 + 1-and the 2 + 1-code coincide, let us define the following norm Consider Fig. 11 where this norm is plotted vs. time. We find very good agreement between the two codes. Some deviations can be expected, since in the 2+ 1-code, Gowdy symmetry is valid only approximately; see more comments in [3]. So far, we have made no effort to explain the oscillatory behavior in both figures, but we conjecture them to be caused by aliasing, since its amplitude becomes smaller, the higher the spatial resolution is. We have also compared more variables at other grid points and found similar results.

Summary and conclusions
The purpose of this paper is to introduce our numerical approach for the study of evolution equations with spatial S 3 -topology. First, we have given details on a geometric point of view in order to find a natural discretization. Second, we have discussed issues related to the implementation and, third, analyzed test applications. By all this we have demonstrated the feasibility of this approach. Indeed, these techniques have been applied elsewhere, for instance in [3,4,5].
Although our main interest lies in the field of general relativity, we believe that the applicability of the method is more general. For the presentation of this paper, we made a couple of special choices which need to be overcome in order to use the method in more general situations. First, our analysis is based on U(1)-symmetry so far. We believe that it is straight forward to generalize in principle. Furthermore, we decided to formulate the equations in terms of smooth global frames on S 3 . This had the consequence that only singular terms of special type are present, and with this knowledge, we were able to work out the Fourier series of the formally singular terms explicitly in Section 2.2. However, in general, as soon as we have a-priori knowledge about the structure of the formally singular terms in the equations, and if their type is not completely different than ours, a similar regularization is possible. Hence, from this point of view, one should be able to modify our method to cases, where other kinds of frames or spatial coordinates are used. In our application, we made a special choice of coordinate gauge and frame transport. This implied that, in particular, the matrix (T a ′ a ) is constant in time. If the original evolution equations imply a well-posed initial value problem, then, in this gauge, also the evolution equations for Gowdy symmetry reduced to 1 + 1 dimensions have this property. For other gauges, however, this might not be the case. A further issue turns up when the frame transport does not keep the orthonormal frame boundary adapted in the sense of Section 3.1.3. Then a different boundary treatment than ours might become necessary. As a last point of this certainly incomplete list, we mention that we have not made experiments with other time integrators yet. Of particular interest might be implicit schemes, for instance, to treat parabolic evolution equations.
We believe that our current numerical infrastructure has not yet been pushed to its limits, but is also not yet really optimized for high spatial resolutions. For instance, we still do not use FFT but only the partial summation scheme. Also it may be true, that there is a more optimal trade-off between accuracy and efficiency for other time integrators than the Runge Kutta schemes of our choice; comments on this can be found in [6]. Furthermore, it might make sense to think about parallelization of the code. This should be straight forward with some publicly available FFT libraries, for instance [17]. We expect, that the 1 + 1-code, in its current form, is limited to a few thousand spatial grid points. We want to stress that this turned out to be sufficient in our runs with T 3topology in [3]. There, we were able to reproduce "spiky features" numerically. Hence, we are optimistic that also such difficult phenomena can be studied with our numerical infrastructure in 1 + 1 dimensions. Dependent on the kind of the application, one may nevertheless doubt, if spectral methods are suitable at all. In general, it is fair to say that, although these methods are highly accurate for lower resolutions, they might be too inefficient for high resolutions. Thus it makes sense, to also investigate into other methods, for instance finite differencing methods, maybe with multipatch [38,11,27], in order to make systematic comparison studies. Further automatic local spatial adaption, usually called adaptive mesh refinement, would be desirable; for instance cosmological solutions with spikes have been investigated with such techniques in [20], however not for S 3 -topology.
In summary, due to the results of the tests here and in [3,4,5], we believe that our numerical approach has promising potentials, which have not been fully exhausted yet, for many kinds of applications with spatial S 3 -topology and beyond. For instance, we find it particularly appealing that our infrastructure can be used directly also for problems with spatial T 3 -topology and, at least for Gowdy symmetry, with S 1 × S 2topology [5]. For the applications, which we have in mind, future research will show, how much the method must be modified, or whether completely new approaches will become necessary, in order to deal with probably much more severe phenomena than those present in the classes of solutions considered so far.