New efficient algorithm for the isometric embedding of 2surface metrics in 3 dimensional Euclidean space
Abstract
We present a new numerical method for the isometric embedding of 2geometries specified by their 2metrics in three dimensional Euclidean space. Our approach is to directly solve the fundamental embedding equation supplemented by six conditions that fix translations and rotations of the embedded surface. This set of equations is discretized by means of a pseudospectral collocation point method. The resulting nonlinear system of equations are then solved by a NewtonRaphson scheme. We explain our numerical algorithm in detail. By studying several examples we show that our method converges provided we start the NewtonRaphson scheme from a suitable initial guess. Our novel method is very efficient for smooth 2metrics.
pacs:
02.70.Hm, 02.60.Lj, 02.40.k, 04.70.Bw1 Introduction
The famous Weyl problem [1] is the question whether each positive definite 2metric with positive intrinsic Gauss curvature given on the unit sphere can be realized as the metric on some 2dimensional surface embedded in three dimensional Euclidean space . This realization is referred to as isometric embedding. Important answers to this question were given first by Lewy [2] and later by Alexandrov and Pogorelov [3, 4] and independently by Nirenberg [5] in his groundbreaking work. The latest answer due to Heinz [6] is in the affirmative if the metric is three times differentiable. Note, however, that these answers are complex existence proofs which cannot easily be used to compute the embedded surface. In order to explicitly construct this surface we have to find a set of equations that can be solved at least numerically for any given metric.
The explicit solution of the problem is of interest for example if we want to visualize black hole horizons. For example current numerical simulations of black holes routinely compute apparent horizons [7, 8]. These horizons are usually found in particular coordinates that are convenient for the stability of the simulations, but horizons plotted in these coordinates can be misleading as deformations can simply be coordinate effects. However, since the intrinsic horizon metric is known in such simulations its isometric embedding in will show the true horizon shape.
Another reason for renewed interest in the isometric embedding problem is the quest for a quasilocal definition of the energy that is enclosed inside a specified 2dimensional surface. A recent definition for this energy by Wang and Yau [9, 10] seems to be particularly promising. It is given in Definition 5.2 of [10] and involves a minimization. The minimum is found by solving Eq. (6.6) of [10]. The latter equation is elliptic and contains the extrinsic curvature of a 2surface that has to be embedded into . Thus in order to find the energy one has to solve an isometric embedding problem as an intermediate step in the calculation. Since elliptic problems are well studied, the hardest part is to find an efficient algorithm that performs this isometric embedding so that the elliptic equation can be solved and used to find the minimum needed for the definition of the energy.
The starting point for the embedding problem is simple and given by
(1) 
where is the vector in that describes the embedded surface as a function of the coordinates (with ) of the given metric on the unit sphere. Since the metric is symmetric Eq. (1) really consists of three equations for the unknown functions . From geometric considerations it is obvious that the embedded surface is unique only up to arbitrary rigid rotations, translations and reflections. In addition, Eq. (1) is not of a well known type such as e.g. hyperbolic, parabolic or elliptic equations. Thus standard methods do not necessarily apply. According to Khuri, little information has been obtained by studying it directly [11]. For this reason Eq. (1) is usually reformulated.
One common approach to reformulating Eq. (1) is to multiply it by with the result
(2) 
Since the left hand side is obviously a flat metric, the metric on the right must also be flat. Setting its intrinsic scalar curvature to zero leads to an equation for given by
(3) 
where the Gauss curvature is obtained from the Ricci Scalar of . This equation is called the Darboux equation and contains at most second derivatives of (because all third derivative terms cancel). The Darboux equation is of MongeAmpere type. It is elliptic if is positive and hyperbolic if is negative. Once is known and can be directly integrated. This approach has been used [12] to study Misner initial data for the collision of black holes.
It is also possible to derive a Darboux equation for
(4) 
It reads [13]
(5) 
This formulation has the advantage that we can obtain the extrinsic curvature (see Eq. (26)) of the embedded surface directly from without having to first compute . The disadvantage of the Darboux equation for is that it is not known what extra conditions we have to impose on to ensure a unique solution. For example a translation of would change the value of as well. Thus for any solution of Eq. (5)) there are infinitely many nearby solutions that can be generated by infinitesimal translations. Without any extra conditions on this approach is not particularly suitable for numerical calculations, because most iterative numerical methods will fail to converge under such conditions.
Another approach [14] uses three dimensional wireframes whose edge lengths are chosen such that they correspond to the lengths between points computed with the metric . In this approach a system of equations for the points in the wire frame is solved directly. The approach has been criticized in [15] for allowing multiple solutions that are not all smooth. For example if the embedded surface is supposed to be a sphere one could obtain a wire frame that corresponds to a sphere whose top third has been inverted, since that would not change the distances of any neighboring points.
Yet another approach [15] is based on expanding the embedded surface in spherical harmonics and minimizing the differences between its metric and . As with any minimization scheme there is the danger that the algorithm gets stuck in a local minimum that does not correspond to the sought after global minimum where both metrics agree. For this reason this approach is computationally intensive.
A more recent approach for numerically solving the embedding problem has been given by Jasiulek and Korzyński [16]. It uses Ricci flow to find conformal relations between the original metric, the round sphere metric and all intermediate metrics. The round sphere metric is embedded in and used as the starting point for an embedding flow back to the original metric. To step from one surface to another, Eq. (1) is linearized and solved for the change in the vector . This method requires an inversion of the extrinsic curvature tensor at each step therefore limiting its application to strictly positive scalar curvature surfaces.
Another recent algorithm called AIM [17] starts from a triangulation of . It finds an approximation to the unique Alexandrov polyhedron for a given triangulated surface with pointwise convex polyhedral metric by using combinatorial Ricci flow followed by an adiabatic pullback with annealing. The main advantage of AIM is that it does not require an initial guess for the embedded surface. The resulting polyhedra do not have inverted regions unlike the wireframes in [14]. While AIM is computationally demanding, it could provide a good initial guess for other algorithms that start from an initial guess for the embedded surface.
2 Numerical method
In this section we present a new numerical method to find the isometric embedding of a 2metric on the unit sphere into . Our approach is to use pseudospectral methods to directly solve Eq. (1). We describe an efficient implementation of this method with the well tested SGRID code [18, 19, 20, 21].
2.1 The pseudospectral collocation method
In one spatial dimension, spectral methods are based on expansions
(6) 
of every field in terms of suitable basis functions with coefficients . Once the coefficients are known it is easy to compute derivatives of from
(7) 
since the derivatives of the basis functions are known analytically.
However, instead of directly storing and manipulating the coefficients up to some desired order in , we make use of the fact that (for the basis functions of interest) we can derive from the values of at certain collocation points. If is known to have the values at the collocation points for it is possible to invert the equations
(8) 
and to exactly solve for the coefficients in terms of the . The location of the different collocation points depends on the basis functions used. For example, for Fourier expansions the collocation points have to be equally spaced in the interval considered. This approach of storing the field’s values at the collocation points is called a pseudospectral collocation method. It has the advantage that nonlinear terms such as can be computed from simple multiplications.
The generalization to two dimensions, is straight forward and can be summarized by
(9) 
I.e. we are using basis functions which are products of functions that depend only on one coordinate. Note that the coordinates and need not be Cartesian coordinates.
In this work we have chosen and to be the standard spherical coordinates and , such that a point on the unit sphere is given by . For both angular directions we use a Fourier basis so that
(10)  
(11) 
where
(12) 
The collocation points are chosen to be
(13)  
(14) 
where
(15)  
(16) 
This choice ensures that there are no collocation points at the coordinate singularities located at and . Notice that, since both and are between 0 and , we have a double covering of the entire domain. This double covering is necessary to ensure the periodicity required for Fourier expansions in both angles.
Any function can then be expressed in this basis as
(17) 
However, in our code we never really uses this expansion to compute all the coefficients . Rather, we only ever expand in one direction and instead use
(18)  
(19) 
to compute the coefficients or along a line in the  or direction. This suffices to compute partial derivatives with respect to our coordinates and .
2.2 Solving the isometric embedding equation
In order to solve the isometric embedding Eq. (1) we represent each component of the vector by its values at the collocation points of Eqs. (13) and (14). This results in unknowns that have to be determined. Since Eq. (1) is symmetric in and we have 3 equations per collocation point which results in equations in total. However, as already mentioned in the introduction, Eq. (1) has a unique solution only up to arbitrary rigid rotations, translations and reflections. This means that not all these equations are independent. In order to fix the position and orientation of the embedded surface we replace six of the equations by the following six conditions. The first three conditions replace the equation for in Eq. (1) at the coordinates with indices , and by
(20) 
Note that here means: round to the closest integer below or at . The conditions (20) ensure that three points with will have . Once is fixed in this way, the embedded surface can still be translated in the  and directions and also rotated about the axis. The two translations are fixed by demanding
(21)  
(22) 
These two conditions replace the equation for at the point and the equation for at the point . They ensure that the embedded surface is centered near . The remaining rotation freedom about the axis is fixed by replacing the equation for at the point with index by
(23) 
Note that has to be computed by spectral interpolation because there is no collocation point at . This last condition fixes the orientation of the embedded surface such that at the point .
The particular components and collocation points at which we replace Eq. (1) by any six conditions that fix translations and rotations do not matter in principle. Our particular choices above are mainly motivated by trying to obtain simple expressions for the conditions. Choosing points near the poles and the equator of the coordinate system is one way of obtaining such simple equations, but other choices are possible.
If we take all these conditions into account we obtain nonlinear equations of the form
(24) 
for the unknown at all collocation points which make up the solution vector . We solve this system of equations by a NewtonRaphson scheme. In order to solve the linearized equations
(25) 
in each NewtonRaphson step, we note that contains spectral derivatives of in different directions, so that the matrix is sparse in the sense that it contains about 95% zeros. We use the sparse matrix solver UMFPACK [22, 23, 24, 25, 26] to numerically solve the linearized Eq. (25). This scheme requires an initial guess, for which we simply use a spherical surface.
We have observed that the NewtonRaphson scheme only reliably converges if both and are chosen to be odd. The reason for this is likely related to the fact that only for an odd number of collocation points we have the same number of sine and cosine functions in the expansions given by Eqs. (10), (11), and (12). For an even number of collocation points, there is only a cosine function at the highest wavenumber. So the derivative of this term (which would be a sine function) cannot be represented with our expansion, which in turn means that the highest coefficient never enters our system of equations and thus cannot be determined.
3 Results
In this section we will use our method to numerically find the embedding of several 2metrics. In order to compare with analytically known results it is useful to compute the extrinsic curvature [13]
(26) 
of the surface embedded in . We now demonstrate our method for a few simple examples.
3.1 Ellipsoid
We first consider the following 2metric:
(27) 
This is the metric of an ellipsoid in with an aspect ratio of . Thus the embedding is known analytically and given by
(28) 
and
(29) 
We have used our numerical method with the metric given in Eq. (3.1) for . We start our algorithm from an initial spherical surface with a radius of . For a resolution of a solution is found after 5 NewtonRaphson steps. The radius of the initial surface is not important. Our method also works with other initial guesses, such as larger radii or offcentered spheres.
The embedded surface for is shown in Fig. 1. We see that our method works as expected and recovers the ellipsoid from the 2metric.
The resulting extrinsic curvature is shown in Fig. 2 together with the intrinsic Gauss curvature.
The norm of its error
3.2 Deformed ellipsoid
The next test is to find an embedding for the metric
(30) 
This metric belongs to a deformed ellipsoid. Thus the analytic solution of the embedding problem is known to be
(31) 
Using the metric given in Eq. (3.2) for , we have started our algorithm from an initial spherical surface with a radius of . For a resolution of a solution is again found after 5 NewtonRaphson steps and the error is again only of order . Again the initial radius is not important and is just the exact value we have used to produce the data for the figures in this subsection. Other initial radii also work, but may need more NewtonRaphson steps. For example for we need 11 NewtonRaphson steps. The embedded surface is shown in Fig. 3.
The resulting curvatures are shown in Fig. 4.
Despite the more irregular shape of the Gauss curvature the algorithm has no problem finding the correct isometric embedding.
3.3 Horizon shape of a Kerr black hole
As already mentioned in the introduction, the isometric embedding of the metric on a black hole horizon can be useful for visualizing the intrinsic horizon geometry. Here we consider the horizon of an axisymmetric Kerr black hole. The metric on the horizon in this case is well known. For a black hole of mass and spin parameter it is given by
(32) 
where
(33) 
is the horizon radius and we are using units of . Note that a horizon exists only for .
Figure 5 shows the Gauss curvature of the horizon for the case of . As we can see it becomes negative near the poles of the black hole. This horizon cannot be embedded in at all. In fact it is well known that only horizons with can be embedded in [15, 27]. Such horizons have everywhere. Thus our algorithm of course can only work if .
For example, for which is less than 1% from the limit our algorithm still works. In this case the embedded surface (shown in Fig. 6) becomes flattened at the poles with nearly zero.
The corresponding curvatures are plotted in Fig. 7.
3.4 Dented sphere
As we have seen from the example of the Kerr black hole horizon, an embedding for into does not always exist. On the other hand there are plenty of closed surfaces in that have negative curvature in some places. Imagine for example a sphere with a dent (with continuous extrinsic curvature). Even though the Gauss curvature is expected to be negative around the dent, the metric on a dented sphere is certainly embeddable in .
In this subsection we consider a particular example of a deformed sphere that while still axisymmetric has a dent at each pole. Our particular metric has the form
(34) 
Its embedding is given by
(35) 
When we use our method with the metric given in Eq. (3.4) and start our algorithm from an initial spherical surface with radius the NewtonRaphson scheme does not converge. This problem occurs for all spherical initial guesses for all radii. We have not fully investigated if it can be circumvented by a more sophisticated algorithm. However, we have tried a modified NewtonRaphson scheme that uses backtracking, i.e. it searches for the optimum step along the direction found by the linear solver but may take a smaller step if that reduces the residual error. This algorithm also fails. It eventually gets stuck because the optimum step becomes a step of length zero.
On the other hand if we start with
(36) 
as initial guess, or something similar, the scheme converges and finds the correct solution. The embedded surface found by the scheme in this case is shown in Fig. 8. We see that this surface has large dents in the polar regions.
In Fig. 9 we show the resulting curvatures. It is obvious that there are large regions where around both poles.
The reason why an initial guess similar to Eq. (3.4) works, while a spherical initial guess fails, may be related to the fact that the guess in Eq. (3.4) already has two dents just like the true embedded surface.
In order to verify our solution we have also conducted a convergence test for the dented sphere case. Notice, however, that the particular 2metric in Eq. (3.4) can be expanded in terms of a finite number of sine and cosine functions. Since our spectral method basically consists of expanding everything in sine and cosine functions, our truncated expansion becomes exact as soon as we expand to high enough order, i.e. if we use enough collocation points. As is evident from Eq. (3.4) the sine and cosine functions with the highest frequency will have as an argument. From Eq. (12) we can thus see that we need at least to include all these functions in our spectral expansions. So a convergence plot is expected to show errors that fall with increasing number of collocation points only for . For the truncation error will be zero and the residual errors will be dominated by numerical round off errors due to the use of floating point arithmetic, i.e. the error should level off at ca. .
In Fig. 10 we show the errors in the extrinsic curvature for different numbers of collocation points . We see the expected approximately exponential convergence for , while as anticipated we reach numerical round off for .
Similar arguments also apply to all the other examples discussed in the previous sections. As soon as is large enough the spectral expansions become exact for all our examples. Notice, however, that this only happens due to the simplicity of our examples. For a generic smooth 2metric the spectral expansions would not be exact for any finite .
4 Discussion
The purpose of this paper is to introduce a new numerical method for the computation of the isometric embedding of 2metrics in . Our method directly solves the embedding Eq. (1) together with 6 conditions (given in Eqs. (20), (21) and (23)) that fix the center and orientation of the embedded surface. We use a pseudospectral collocation point method to represent the derivatives in Eq. (1). This results in three nonlinear equations per collocation point for the three components of the embedding vector. This set of equations is then solved with a NewtonRaphson scheme.
As we have seen it is important to work with an odd number of collocation points in both coordinate directions if one uses Fourier expansions. As long as the given 2metric has positive Gauss curvature everywhere our method seems to always converge if we initialize the NewtonRaphson scheme with a spherical surface as initial guess. For cases where the Gauss curvature is negative in some regions, but an isometric embedding is known to exist, our method can fail if we use a spherical surface as initial guess for the NewtonRaphson scheme. Nevertheless the scheme can be made to converge by choosing a better initial guess.
Our algorithm is very efficient. On a laptop computer for it only takes about half a second to find the isometric embedding of the deformed metric of Eq. (3.2) shown in Fig. 3. For comparable accuracies this is orders of magnitude faster than the AIM algorithm [17]. However, since most other previous methods have used different examples and since the papers describing them do not give exact timing information we cannot directly compare. Nevertheless, our method has the advantage of not allowing multiple solutions, so that it does not suffer from the drawbacks that can cause problems for wireframe approach [14]. We also expect it to be faster than the minimization scheme discussed in [15].
Footnotes
 The norm of the error is not normalized. For a component like we compute it from where we sum over all grid points.
References
 Hermann Weyl. Über die Bestimmung einer geschlossenen konvexen Fläche durch ihr Linienelement. Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, 61:148–178, 1916.
 Hans Lewy. On the existence of a closed convex surface realizing a given Riemannian metric. Proceedings of the National Academy of Sciences of the United States of America, 24:104–106, 1938.
 Aleksandr Danilovic Aleksandrov. Existence of a polyhedron and of a convex surface with a given metric. Doklady Akademii Nauk SSSR, 50:103–106, 1941.
 Aleksei Vasil’evich Pogorelov. On the proof of Weyl’s theorem on the existence of a closed analytic convex surface realizing an analytic metric with positive curvature given on the sphere. Uspekhi Matematicheskikh Nauk, 4:183–186, 1949.
 Louis Nirenberg. The Weyl and Minkowski problems in differential geometry in the large. Communications on Pure and Applied Mathematics, 6:337–394, 1953.
 Erhard Heinz. On Weyl’s embedding problem. Journal of Mathematics and Mechanics, 11:421–454, 1962.
 Luis Lehner. Numerical relativity: a review. Classical and Quantum Gravity, 18(17):R25, 2001.
 Jonathan Thornburg. Event and apparent horizon finders for 3+1 numerical relativity. Living Reviews in Relativity, 10(3), 2007.
 MuTao Wang and ShingTung Yau. Quasilocal mass in general relativity. Phys. Rev. Lett., 102:021101, Jan 2009.
 MuTao Wang and ShingTung Yau. Isometric embeddings into the minkowski space and new quasilocal mass. Communications in Mathematical Physics, 288(3):919–942, 2009.
 Marcus A. Khuri. On the local solvability of darboux’s equation. Discrete and Continuous Dynamical Systems Supplement, pages 451–456, 2009.
 Joseph D. Romano and Richard H. Price. Embedding initial data for black hole collisions. Class. Quantum Grav., 12:875–893, 1995.
 Qing Han and JiaXing Hong. Isometric Embedding of Riemannian Manifolds in Euclidean Spaces, volume 130. American Mathematical Society, 2006.
 HansPeter Nollert and Heinz Herold. Visualization in curved spacetimes II. visualisation of surfaces via embedding. In F. W. Hehl, R. A. Puntigam, and H. Ruder, editors, Relativity and Scientific Computing, page 330, Berlin, 1998. Springer Verlag.
 Mihai Bondarescu, Miguel Alcubierre, and Edward Seidel. Isometric embeddings of black hole horizons in threedimensional flat space. Class. Quantum Grav., 19(2):375–392, 21 January 2002.
 Michael Jasiulek and Miko Korzynski. Isometric embeddings of 2spheres by embedding flow for applications in numerical relativity. Class.Quant.Grav., 29:155010, 2012.
 Shannon Ray, Warner A. Miller, Paul M. Alsing, and ShingTung Yau. Adiabatic Isometric Map to Embed Polyhedral Metrics in Euclidean 3Space. 2014. in preparation.
 Wolfgang Tichy. Black hole evolution with the BSSN system by pseudospectral methods. Phys. Rev., D74:084005, 2006.
 Wolfgang Tichy. A new numerical method to construct binary neutron star initial data. Class. Quant. Grav., 26:175018, 2009.
 Wolfgang Tichy. Long term black hole evolution with the BSSN system by pseudospectral methods. Phys. Rev., D80:104034, 2009.
 Wolfgang Tichy. Constructing quasiequilibrium initial data for binary neutron stars with arbitrary spins. Phys.Rev., D86:064024, 2012.
 Timothy A. Davis and Iain S. Duff. An unsymmetricpattern multifrontal method for sparse LU factorization. SIAM J. Matrix Anal. Applic., 18(1):140–158, 1997.
 Timothy A. Davis and Iain S. Duff. A combined unifrontal/multifrontal method for unsymmetric sparse matrices. ACM Trans. Math. Softw., 25(1):1–20, 1999.
 Timothy A. Davis. Algorithm 832: Umfpack v4.3—an unsymmetricpattern multifrontal method. ACM Trans. Math. Softw., 30(2):196–199, 2004.
 Timothy A. Davis. A column preordering strategy for the unsymmetricpattern multifrontal method. ACM Trans. Math. Softw., 30(2):165–195, 2004.

Timothy A. Davis.
UMFPACK a sparse linear systems solver using the Unsymmetric
MultiFrontal method:
http://www.cise.ufl.edu/research/sparse/umfpack/.  Larry Smarr. Surface geometry of charged rotating black holes. Phys. Rev. D, 7:289, 1973.