text:
# Dynamical control of the constraints growth in free evolutions of Einstein’s equations.

## I Motivation and overview

The Einstein equations are typically solved in what is called free or unconstrained evolutions. The equations are split into a set of evolution equations and another one of constraints. The standard procedure is to use initial data that satisfy the constraints and later solve only for the set of evolution equations. The logic is that by virtue of Bianchi identities the constraints should be satisfied at all times if they do so initially and the evolution equations are solved. This is actually true only in the absence of boundaries, or in the future domain of dependence of the initial data; in the presence of boundaries appropriate boundary conditions have to given in such a way that the associated initial-boundary value problem is well posed, a highly non-trivial issue on its own (see boundaries  for related work).

However, problems arise in free evolutions even neglecting the presence of boundaries. Although the constraints should be exactly satisfied at the continuum level, in typical, fixed resolution, numerical simulations of strong fields the constraints quickly grow and eventually the code crashes. For a numerically stable scheme (throughout this paper, the term numerical stability is used as equivalent to convergence, in the sense of Lax’s theorem 111That is, the existence of a precise bound for the numerical solution in terms of the initial data, with the bound being independent of resolution.) this growth should go away with resolution, but very high resolutions are usually needed if one wants to keep these errors under control by just adding more points to the simulation. Therefore this procedure is not practical and it might not be even feasible.

Sometimes it happens there are no growing modes at the continuum, but the discretization one is using, even if numerically stable, can introduce errors that grow quickly in time. Examples of this potential source of instabilities are shown in Ref. exc ; it is there also explained how to rearrange the semidiscrete equations in order to prevent this when there is a “conserved” quantity at the continuum (for example, in the case of wave propagation on a stationary spacetime the physical energy would play this role).

The picture that has emerged in other cases is that rapidly growing discrete constraint violating modes are numerical excitations of unstable continuum modes. In other words, solutions of Einstein’s evolution equations that are initially slightly off the constraint surface, but then deviate from it very quickly. The standard approach is then to seek formulations that are as stable as possible under continuum constraint violations. A trivial observation here is that when rapidly growing discrete constraint violations appear it is because the given numerical method allows for them.

Therefore, I argue that it is desirable to control the discretized constraints, taking into account possible growth introduced both by the formulation of the equations, and by the numerical method. The purpose of this paper is to present a method for doing so.

As one is interested in solutions that satisfy the constraints, the Einstein evolution equations can be written in infinitely different ways by changing the system’s behavior off of the constraint surface. For definiteness I will now concentrate on first order (both in space and time) formulations, though the ideas of this paper clearly do not depend on this and can be applied to second order formulations as well. Beginning with any formulation, e.g. a symmetric hyperbolic one, one can add to the right hand side (RHS) terms proportional to the constraints. That is

$$\dot{u}=\underset{j}{}A^j(u,t,\stackrel{}{x})_ju+B(u,t,\stackrel{}{x})+\mu C$$
(1)
where $`B`$ and $`u`$ are a vector valued functions (containing the metric and related variables),

$$C=\left(C_1\mathrm{}C_n\right),$$
(2)
with $`C_i=C_i(u,_ju)`$, is a vector valued constraint (containing the Hamiltonian, momentum, and perhaps other non-physical constraints that appear as integrability conditions when the system is written in first order form), and $`A,\mu `$ are matrix-valued functions. The constraints (2) in Einstein’s equations are quasilinear. That is, nonlinear in the field variables $`u`$, but linear in the spatial derivatives, $`_iu`$. Adding the constraints in a linear way, as in equation (1), yields a quasilinear system of hyperbolic partial differential equations, for which many results from the mathematical literature apply 222During preparation of this manuscript, related by work by David Fiske appeared fiske , in which the motivation is to make the constraints decay exponentially to zero for any spacetime (i.e. without any a priori knowledge). Presumably the resulting equations are not quasilinear (or linearly degenerate in the vacuum case)..

The freedom in choosing $`\mu `$ was originally used to cast Einstein’s equations in symmetric hyperbolic form, later to obtain “physical” characteristic speeds (see hyp  for reviews), and more recently to improve the stability properties of formulations around fixed backgrounds (the number of papers in the area is too large to be reviewed here; see stability  and references therein). Typically $`\mu `$ is chosen as a constant matrix; however, it does not appear to be the most effective way of controlling deviations off the constraint surface in a generic case for two reasons: i) some information about the solution (typically, the assumption that it will be close -or actually identical - to some known background) is usually needed in order to choose an optimal $`\mu `$. However, this kind of information will generically not be available. ii) In principle, constants do not seem the best way to control varying fields.

Choosing constant values for $`\mu `$ seems to be favored for historical reasons, and it is actually not required. Symmetric hyperbolicity and physical characteristic speeds, for example, can be achieved even if using functions of spacetime, provided they are a priori given. See st  for one example.

## II Dynamical minimization of the constraints growth: the main idea

Consider the constraints $`C=C(u,D_iu)`$, where the operator $`D`$ may be either a continuum or discrete derivative. Later $`D`$ will represent discrete derivatives, since the motivation is to control not only continuum constraint violations, but also violations that arise in the discrete equations. For any slicing of spacetime $`S_t\times `$, define a norm $`N_c=N_c(t)`$ for the constraints. Note I will always refer to this norm, not the norm associated with the main field variables (will comment on this in the last section). For simplicity I choose the $`L_2`$ norm,

$$N_c(t)=\frac{1}{2}_{S_t}\underset{i}{}C_i^2$$

where integration is on the spatial hypersurface $`S_t`$ 333If the derivative is discrete this is to be replaced by a discrete sum, of course..

One wants this norm either not to grow during an evolution, or least to grow as slowly as possible. As $`C`$ is a function of the field variables, any solution to the evolution equations automatically determines the rate growth for $`N_c`$, independent of whether one is considering the fully discrete, semidiscrete or continuum systems. I will concentrate on the last two cases (i.e., with continuous time). In principle the dependence of $`N_c`$ on the evolution equations is complicated. However, given that $`\mu `$ is added linearly to the right hand side of the evolution equations, as in Eq.1, the time derivative of $`N_c`$ can be written as

$`\dot{N}_c`$ $`=`$ $`{\displaystyle _{S_t}}{\displaystyle \underset{i}{}}C_i\dot{C_i}`$
$`=`$ $`{\displaystyle _{S_t}}{\displaystyle \underset{i,j}{}}C_i\left[{\displaystyle \frac{C_i}{u_j}}\dot{u_j}+{\displaystyle \underset{k}{}}{\displaystyle \frac{C_i}{(D_ku_j)}}D_k\dot{u_j}\right]`$
Using the evolution equations Eq.(1) and allowing $`\mu `$ to depend only on time (not on space), the time derivative of the norm can then be written as

$$\dot{N}_c=I^{(1)}+\text{trace}(\mu \times I^{(2)})$$
(3)
where

$`I^{(1)}`$ $`=`$ $`{\displaystyle }{\displaystyle \underset{i,j}{}}C_i[{\displaystyle \frac{C_i}{u_j}}+{\displaystyle \underset{k}{}}{\displaystyle \frac{C_i}{D_ku_j}}D_k]\times `$ (4)
$`\left[{\displaystyle \underset{l}{}}(A^lD_lu_j)+B_j\right]`$
$`I_{jl}^{(2)}`$ $`=`$ $`{\displaystyle \underset{i}{}C_i\left[\frac{C_i}{u_j}+\underset{k}{}\frac{C_i}{D_ku_j}D_k\right]C_l}`$ (5)
The quantities $`\{N_c,I^{(1)},I^{(2)},\mu \}`$ depend only on time, in the case of $`\mu `$ by assumption and in the others because spatial integration is carried out. The splitting (3) is a crucial step in the method presented here, since one gains analytical control on the norm growth as a function of $`\mu `$. More explicitly: $`\dot{N_c}`$ is a linear function of $`\mu `$ (because the constraints are added linearly to the evolution equations).

The idea now is to choose $`\mu `$ such that $`\dot{N_c}`$ has some desired behavior within the freedom allowed by a given formulation of the evolution equations. As a simple example, consider a single scalar equation. One could choose $`\dot{N}_c=0`$ by defining $`\mu =I^{(1)}/I^{(2)}`$. Or one could choose the norm to decay at a constant or exponential rate by defining $`\mu =(I^{(1)}+d^2)/I^{(2)}`$ or $`\mu =(I^{(1)}+d^2N_c)/I^{(2)}`$, respectively, with $`d`$ some constant. On the other hand, if the different components of the $`\mu `$ matrix are restricted to some sets (this is a restriction that appears quite often when requiring symmetric hyperbolicity, physical characteristic speeds, or both) one could define $`\mu (t)`$ as the one that gives the “smallest” growth of $`\dot{N}_c(t)`$ in those sets, etc.

This could be done at the continuum. However, the problem is that $`\mu `$ would then depend on the constraints, which in turn depend on the derivatives of the field variables. As a consequence, the resulting evolution equations would not be quasilinear, and it is far from clear that they would define a well posed initial, or initial-boundary, value problem. Indeed, in Section IV I present numerical evidence that strongly suggests that the resulting system would be ill posed.

Therefore, I propose to define $`\mu `$ through a single resolution run and interpolate the obtained discrete function in order to have it defined at all times. In an actual computation one would choose a fixed resolution and dynamically compute $`I^{(1)},I^{(2)}`$ and thus a $`\mu `$ that gives the desired norm behavior, making sure that the resulting $`\mu (t)`$ is within the range allowed by symmetric hyperbolicity (symmetric hyperbolicity is not a requirement of the method, but ensures well posedness of the initial-boundary value problem, provided appropriate boundary conditions are given). In order to ensure numerical stability this $`\mu `$ matrix valued function has to be kept fixed at other resolutions. $`\mu `$ will depend on the original resolution used to define it, but this is not a problem. The important thing is to keep the constraints under control for a given resolution. They will also remain under control for better resolutions if one is in the convergence regime.

$`\mu `$ will also depend on the given problem of interest, and on the numerical method used to solve the equations. That is, for any chosen desired behaviour for $`\dot{N}_c`$, there will be one $`\mu `$, and one well posed formulation of the problem for each physical situation one wants to solve for with a given numerical method. This is not a practical problem and, in fact, it seems difficult to find a formulation that has optimum stability properties for all possible solutions and all possible numerical methods.

## III Controlling the constraints norm growth: more details

If possible, one might want the constraints norm growth to be identically zero at all times $`\dot{N}_c=0`$. A potential problem with this approach is that the resulting $`\mu `$ might have large values, rapid variations in time, or both. The resulting set of PDE’s would then be stiff, having equivalent short timescales, in which case a very small Courant factor would be required. In the next section I show an example where this happens and argue that better results may be obtained by not enforcing $`\dot{N}_c=0`$. Even if one has enough freedom to enforce $`\dot{N}_c=0`$, it appears that allowing $`N_c`$ for some fluctuations around the initial value gives evolutions with better stability. I now elaborate on one way of doing so, assuming one has enough freedom in $`\mu `$. This approach is highly non unique, and many variations of the fundamental idea presented in this paper seem possible, and worth further exploration. If one chooses

$$\dot{N_c}=aN_c$$
(6)
with $`a>0`$, any violation in the constraints will decay exponentially

$$N_c(t+\mathrm{\Delta }t)=N_c(t)e^{a\mathrm{\Delta }t}$$
(7)
Choosing a constant $`a`$ makes the constraints decay for all times, I will discuss this in the last section. Here I take a different approach. I choose a a tolerance value for the norm, $`N_c=T`$, and solve for $`\mu `$ such that the constraints decay to this tolerance value after a given relaxation time. More precisely, I choose $`a`$ such that after time $`n_a\mathrm{\Delta }t`$ the constraints will have the value $`T`$. Replacing $`N_c(t+\mathrm{\Delta }t)`$ by $`T`$ in equation (7) and solving for $`a`$ gives

$$a(t)=\frac{1}{n_a\mathrm{\Delta }t}\mathrm{ln}\left(\frac{T}{N_c(t)}\right).$$
(8)
I now solve $`\dot{N_c}=aN_c=I^{(1)}+\text{trace}(\mu I^{(2)})`$ for $`\mu `$. One solution (non unique, since the equation is scalar and $`\mu `$ is a matrix) is, assuming $`I^{(2)}`$ is invertible, and ommiting all the indexes,

$$\mu =\frac{\left[I^{(2)}\right]^1\left(aN_c+I^{(1)}\right)}{dim\left[I^{(2)}\right]},$$
(9)
with $`a(t)`$ given by eq.(8). Thus, for any value of $`N_c`$ at time $`t`$, at time $`t+n_a\mathrm{\Delta }t`$ the value will be $`T`$. Later in this paper $`\mathrm{\Delta }t`$ will be the discrete time step and, in particular, if $`n_a`$=1, then $`N_c=T`$ at each time step.

It must be emphasized that these results hold in the semidiscrete case. I.e., for any \- not necessarily high - spatial resolution. But since time has been assumed continuous, the Courant factor must be such that the fully discrete simulation faithfully represents the above semidiscrete calculations (another option would be to perform a fully discrete analysis). In the next section I present some numerical experiments to study, among other issues, the extent to which the fully discrete system represents the semidiscrete analysis. As we will see, at least in the examples here considered, standard values of the Courant factor already yield good results.

## IV Some numerical examples

In these numerical experiments I use the method of lines and finite differencing. The difference operator chosen is the simplest one that satisfies summation by parts (second order in the interior, and first order in the boundaries), numerical dissipation is introduced taken into account modifications at boundaries, and third order Runge Kutta is used as time integrator (see exc ).

### IV.1 The equations

I use the standard ADM equations in spherical symmetry with the exact lapse-area locking choices for the lapse and shift area  to illustrate and study the method presented here. This choice of lapse and shift amounts to specifying in an arbitrary, but a priori, way the lapse as a function of spacetime, and the shift as $`\beta =\alpha rK_b`$. This choice of shift implies that the radial area does not change in time, $`\dot{g}_{\theta \theta }=0`$, and one can choose $`r=g_{\theta \theta }^{1/2}`$ as a coordinate. The 3-metric and extrinsic curvature in coordinates $`r,\theta ,\varphi `$ then take the form $`g_{ij}=diag(a^2,r^2,r^2)`$, $`K_j^i=diag(K_a,K_b,K_b)`$.

The evolution equations, with the product of the Hamiltonian constraint and $`\mu =\mu (t)`$ added to the time derivative of the $`K_a`$ equation, are

$`\dot{a}`$ $`=`$ $`\alpha r(aK_b)^{}+a\left[\alpha (K_bK_a)+rK_b\alpha ^{}\right]`$ (10)
$`\dot{K_a}`$ $`=`$ $`(\alpha ^{}r+2\alpha )a^3r^1a^{}a^2\alpha ^{^{\prime \prime }}+`$ (11)
$`\alpha \left[rK_bK_a^{}+K_a(K_a+2K_b)\right]+\mu H`$
$`\dot{K_b}`$ $`=`$ $`\alpha \left(a^3ra^{}+rK_bK_b^{}\right)a^2r^1\alpha ^{}+`$ (12)
$`\alpha \left[K_b(K_a+2K_b)+r^2(1a^2)\right]`$
These equations constitute a strongly hyperbolic system, and the characteristic speeds are “physical”(along the light cone or normal to the spatial hypersurfaces): $`\beta ,\beta \pm \alpha /a`$, for any function $`\mu (t)`$. That is, the equations are strongly hyperbolic even for $`\mu =0`$, i.e., the standard ADM equations in spherical symmetry with this choice of gauge are strongly hyperbolic.

The Hamiltonian and momentum constraints are, respectively,

$`H`$ $`=`$ $`{\displaystyle \frac{4}{ra^3}}a^{}+{\displaystyle \frac{2}{r^2}}(1a^2)+2K_b(K_b+2K_a)`$ (13)
$`M`$ $`=`$ $`K_b^{}+{\displaystyle \frac{K_bK_a}{r}}`$ (14)
I perform these experiments with the Painlevé-Gullstrand (PG) slicing of the Schwarzschild black hole spacetime,

$$a=1,K_a=\frac{\beta }{2r},K_b=\frac{\beta }{r}$$
(15)
$$\alpha =1,\beta =\alpha rK_b$$
(16)
This is a stationary, exact solution, of the evolution and constraint equations, (10,11,12), and (13,14), respectively.

As boundary conditions I set the time derivative of the incoming characteristic modes to zero. This will define a well posed initial-boundary value problem but in general will violate the constraints. One possibility would be to derive constraint-preserving boundary conditions for this problem. However, this is not the purpose of this paper but, rather, to devise a mechanism to control the constraint growth having fixed, among other things, the boundary conditions. This is not meant to be a replacement for constraint-preserving boundary conditions but, instead, complementary to them.

### IV.2 Numerically stable but with fast growing errors simulations

Evolutions of the PG spacetime with equations (10,11,12) can give rise to errors that grow fast in time, even with a stable numerical method.

Figure 1 shows $`N_c`$ vs. time for evolutions of PG initial data, at different resolutions. Typical numerical parameters are chosen: Courant factor $`\lambda =0.25`$, inner and outer boundaries at $`r_i=M,r_o=40M`$, respectively, and dissipation parameter $`\sigma =0.1`$. The constraints converge to zero with resolution, but at fixed resolution they grow very fast in time. This is not a peculiarity of this formulation and this numerical method but, indeed, is typical for free evolutions of Einstein’s equations in the strong field regime.

One could attempt to modify the numerics in order to minimize the growth of the constraints, but I do not pursue this approach here.

An interesting feature of the instabilities of the simulations shown in Figure 1 is that the norm growth is initially negative, but it becomes positive after a single time step, triggering the instability, see (17) (there, as in the rest of the paper, $`\dot{N}_c`$ is computed through the semidiscrete expression (3)).

| $`\mathrm{\Delta }r`$ | $`\mathrm{\Delta }t`$ | $`\dot{N}_c(0)`$ | $`\dot{N}_c(\mathrm{\Delta }t)`$ |
| --- | --- | --- | --- |
| $`M/5`$ | $`5.0\times 10^2`$ | $`5.01\times 10^3`$ | $`1.47\times 10^2`$ |
| $`M/10`$ | $`2.5\times 10^2`$ | $`1.81\times 10^3`$ | $`3.85\times 10^3`$ |
| $`M/20`$ | $`1.25\times 10^2`$ | $`5.39\times 10^4`$ | $`9.54\times 10^4`$ |
| $`M/40`$ | $`6.25\times 10^3`$ | $`1.46\times 10^4`$ | $`2.34\times 10^4`$ |
| $`M/80`$ | $`3.12\times 10^3`$ | $`3.82\times 10^5`$ | $`5.76\times 10^5`$ |

(17)
A static choice of $`\mu `$ cannot account for these fluctuations. In particular, choosing $`\mu `$ using only the growth rate around the background solution, which in this case would be the initial data, and not the errors introduced during numerical integration, in general will not be able to correct these kind of fluctuations in $`\dot{N}_c`$.

### IV.3 Strictly enforcing semidiscrete norm preservation

I now discuss some numerical experiments using a dynamically calculated $`\mu (t)`$ such that $`\dot{N_c}=0`$ in the semidiscrete case. For a chosen resolution, $`\mu `$ is defined at each time step by

$$\mu (n\mathrm{\Delta }t)=\{\begin{array}{c}I^{(1)}(n\mathrm{\Delta }t)/I^{(2)}(n\mathrm{\Delta }t),\text{ if }I^{(2)}(n\mathrm{\Delta }t)0\hfill \\ \mu \left((n1)\mathrm{\Delta }t\right),\text{ if }I^{(2)}(n\mathrm{\Delta }t)=0\text{ and }n1\hfill \\ 0,\text{ if }I^{(2)}(n\mathrm{\Delta }t)=0\text{ and }n=0\hfill \end{array}$$
(18)
and $`\mu (n\mathrm{\Delta }t)`$ is in turn used to compute the RHS needed to advance $`u\left(n\mathrm{\Delta }t\right)`$.

Figure 2 shows the resulting norm as a function of time, for evolutions enforcing $`\dot{N}_c=0`$ at different resolutions. The code blows up rather fast, and this happens at earlier times when resolution is increased. This last feature should not be seen as an indication of ill-posedness of the method, since Figure 2 should not be seen as a convergence test. For such a test one should dynamically define $`\mu (t)`$ with a given resolution and then keep that $`\mu `$ fixed for all other runs. Figure 3 shows the result of doing so, defining $`\mu (t)`$ through a run with $`\mathrm{\Delta }r=M/5`$ and then keeping it fixed for a convergence test 444Since $`\mu `$ is obtained for the coarsest resolution used in the convergence test, some interpolation procedure is needed in order to have this function defined at intermediate time steps. The simple procedure I have followed here is to define the needed intermediate values of $`\mu `$ between time steps $`n`$ and $`n+1`$ just as $`\mu (n)`$, though better (e.g. trigonometric) interpolation could be used. Another option would be to define $`\mu `$ through the finest resolution that is going to be used for the convergence test and later coarsen the obtained discrete $`\mu `$ function.. As expected, and in contrast to Figure 2, the norm does decrease as the resolution is increased. Figure 2 strongly suggests that if one insisted in dynamically defining $`\mu `$ at each resolution, instead of fixing it at a single resolution, one would end up with an ill-posed problem. On the other hand, fixing $`\mu `$ leads to a strongly hyperbolic problem and no convergence problems appear.

Considering these runs with $`\dot{N}_c=0`$, it appears that enforcing semidiscrete norm preservation is “too rigid” (recall the discussion at the beginning of Section III). Figure 4 shows the $`\mu `$ functions obtained through the runs of Figure 2. Notice that when increasing resolution the $`\mu `$ that is needed in order to preserve the norm increases quite fast in absolute value. For example, with $`\mathrm{\Delta }r=M/80`$ the initial values of $`\mu `$ are of order $`10^2`$. Since $`\mu `$ appears in the principal part of the equations, having a large value implies the need of a small Courant factor to follow the simulation; while Figures 2 and 4, on the other hand, where obtained using the same Courant factor when changing resolution. A Fourier decomposition of the numerical solution, or an explicit von-Neumann analysis of a linear, constant coefficient problem could give further insight into this question.

A second possibility would be to perform a fully discrete analysis (as opposed to the semidiscrete one used in this paper), and to enforce strict conservation of the fully discrete norm, but these issues are not pursued here.

### IV.4 Dynamically achieving a tolerance value

I now set a tolerance value $`T`$, and choose $`\mu `$ such that $`N_c=T`$ after a given number of timesteps. In doing so I use the semidiscrete expressions (8,9). If $`I^{(2)}`$ is zero, $`\mu `$ will be copied from its previous value, or set to zero if this happens in the first iteration.

Figure 5 shows results with a dynamically defined $`\mu `$, using the same resolution, domain, Courant and dissipation factors as before. $`n_a`$ is chosen to be one, corresponding to the norm returning to the specified tolerance value after every time step.

The discretized constraints are, of course, not zero initially, even when the initial data is analytic, because computing them involves finite differencing the initial data. The initial value for the norm for this resolution is $`N_c=4.8\times 10^4`$. Therefore I choose a tolerance value $`T=5\times 10^4`$, to keep the constraints close to its initial value.

Figure 6 shows the resulting norms for a convergence test performed with the same resolutions used in Figure 1, and the $`\mu `$ shown in Figure 5. The code runs for $`10^4M`$ without any sign of instabilities, even with the coarse resolution that was used to define $`\mu `$.

It could happen that one controls the constraints at the price of introducing other errors. This does not happen, at least in the numerical experiments here considered. In these experiments, not only does the norm associated with the constraints remains close to its initial, truncation value, but the same happens with the norm associated with the main field variables, $`N=(a^2+K_a^2+K_b^2)`$. Since the exact solution is stationary, $`N`$ being close to its initial value means that the method does not introduce new errors. Figure 7 shows $`N_c(t)`$ and $`N(t)`$ for the run used to define $`\mu `$.

#### IV.4.1 Dependence of $`\mu `$ on the tolerance value

In the previous examples I chose the tolerance value, $`T`$, to be roughly the initial, truncation error for $`N_c`$. Although this seems the natural thing to do, I will now discuss sensitivity of the stability method with respect to the chosen value of $`T`$.

Figure 8 shows results when choosing different tolerance values to define $`\mu `$, but keeping all other numerical parameters (in particular, the resolution, i.e. Figure 8 is not a convergence test) fixed. When $`T`$ is much larger than the initial truncation value of $`N_c`$, large errors allow triggering of nonlinear instabilities. On the other hand, values of $`T`$ considerably smaller than the initial truncation one forces $`\mu `$ to change very fast, and to large values. In either case the code crashes. The $`\mu `$ functions for two such cases are shown in Figure 9.

Figure 10 shows the details of Figure 8 near $`t=0`$. Notice that not only all the runs begin with the same (truncation) value for $`N_c`$ (this is expected, since the same resolution is used), but $`N_c`$ is also the same in the very first timestep. The reason for this is that $`\mu =0`$ initially for all the runs, because $`I^{(2)}=0`$ in the initial data. This is a coincidence of the initial data being evolved in this example: the only discretization error in the computation of the Hamiltonian constraint in this model comes from the finite differencing of $`a`$, which in this example is initially $`1`$ for all grid points.

### IV.5 Intentional constraint violations

As a final numerical example, now I consider initial data that violates the constraints at the continuum. The equations still use the exact lapse and area-locking shift given by Eq.(16), but now the initial data are given by

$`a`$ $`=`$ $`\left(1+s_1e^{(rr_0)^2}\right)`$ (19)
$`K_a`$ $`=`$ $`{\displaystyle \frac{\beta }{2r}}\left(1+s_2e^{(rr_0)^2}\right)`$ (20)
$`K_b`$ $`=`$ $`{\displaystyle \frac{\beta }{r}}=\left(1+s_3e^{(rr_0)^2}\right)`$ (21)
When $`s_i=0`$, this data reduces to that of the PG black hole. The constants $`s_i`$ are used to introduce an “asymmetry” in the perturbations of the PG initial data. The perturbations are rather large, with $`s_1`$ or order unity. In this example, $`s_1=1,s_2=0.5,s_3=0.7`$. When $`\mu =0`$ the code crashes before $`t=4M`$ at a resolution $`\mathrm{\Delta }r=M/5`$. At this resolution $`N_c`$ for the initial data is $`N_c(0)=2.66\times 10^3`$ (one order of magnitude larger than the associated truncation error for the case $`s_i=0`$, which for this resolution is $`N_c=4.8\times 10^4`$ ). Thus, I choose a tolerance value of $`T=5\times 10^3`$, $`n_1=10`$ in order to avoid fast variations, and rerun, dynamically obtaining $`\mu (t)`$. Results for the unmodified equations ($`\mu =0`$), and with a dynamical $`\mu `$ are shown in Figures 11, 12. The difference in stability is quite striking, especially considering that the perturbation of the PG initial data is so large.

## V Remarks

There is a lot of interest in finding formulations of Einstein’s equations with good stability properties off of the constraint surface. So far, most, if not all, of the work has been concentrated in stability properties around known solutions. The purpose of this paper has been to introduce new ideas that hopefully will be helpful in generic evolutions, where one has no a priori knowledge of the solution. This is one fundamental difference with previous approaches.

Another difference is that previous work has concentrated on stability properties at the continuum. Although this is an important issue, controlling constraint violating modes introduced by the numerical method is also crucial. Thus, this proposal aims to control discrete constraint violations introduced both by the underlying formulation of equations at the continuum, and whatever numerical method one has chosen.

Sometimes it seems that one would ideally want the constraints to decay exponentially to zero lambda ,fiske . Although this might be attractive at the continuum, in what concerns the discrete constraints it might be better to maintain them around the initial, truncation error. For example, beginning with initial data corresponding to a stationary spacetime, it seems difficult to imagine that by making the constraints decrease in time through numerical integration (keeping resolution fixed), one would end up with a numerical solution that has smaller errors than the initial data, which was computed by just numerically evaluating the exact solution at given gridpoints. Thus, I have concentrated here on keeping the discrete constraints at some tolerance value, instead of forcing them to decay to zero. However, there may be other scenarios where it is preferable to make them decay to values smaller than the initial one, and the method here presented allows for this as well. Along these lines, an extension of the work presented here is to use the evolution equations applied to initial data off the constraint surface to produce solutions of the constraints as an alternative to relaxing the constraints themselves 555L. Lehner, private communication.

In some other approaches ls , a norm for the main variables, say $`N(t)=_{S_t}_iu_i^2`$, is minimized with respect to constraint violating perturbations (around a given, known background). I could have here chosen to dynamically minimize this quantity. One potential problem with this is that it is not clear how much minimization makes the solution closer to the “exact” one. For example, in the spherically symmetric example here considered the time derivative of the norm for the main variables has the form

$$\dot{N}=\stackrel{~}{I^{(1)}}+\mu \stackrel{~}{I^{(2)}}$$

Since $`\mu `$ is completely free in this model, one is able to make $`\dot{N}`$, for some chosen but arbitrary resolution, as negative as one wants. What happens then is that all the field variables decay to zero (numerical experiments that I have done do confirm this). Clearly this is not what one wants, since the “exact” solution is not identically zero. On the other hand, one could also attempt to keep $`N`$ close to its initial, truncation value. However, this would be a good idea only if the solution is stationary, since otherwise there is no reason for the exact $`N(t)`$ to be close to $`N(0)`$ for all times. Thus I have here considered a norm associated with the constraints, since they are analytically zero, regardless of the solution.

## VI Acknowledgments

I would like to thank Oscar Bruno, Gioel Calabrese, Luis Lehner, David Neilsen, Jorge Pullin, Oscar Reula, Olivier Sarbach, Ed Seidel and Jonathan Thornburg for very helpful discussions, comments and suggestions. This work was supported in part by NSF grant PHY9800973, the Horace Hearne Jr. Institute for Theoretical Physics, and Fundación Antorchas.


format: text

text:
# Global structure of the Zipoy-Voorhees-Weyl spacetime and the 𝛿=2 Tomimatsu-Sato spacetime

## I Introduction

Almost thirty years have passed since Tomimatsu and SatoTomimatsu.A&Sato1972  found a series of exact solutions, the so-called TS solutions, to the vacuum Ernst equation, which can be regarded as representing quadrupole deformations of the Kerr solution. However, the global structure of these solutions has not been studied well, or at least, its details are not known widely, though some basic features such as the existence of naked ring singularities and causality violating regions are well knownTomimatsu.A&Sato1973 ; Tomimatsu.A&Sato1973a ; Gibbons.G&Russell-Clark1973 . Furthermore, some conflicting or wrong statements on their global and singularity structures are sometimes found in the literature for the case the deformation parameter $`\delta `$ is even. For example, Gibbons and Russell-ClarkGibbons.G&Russell-Clark1973  claimed that the TS (Tomimatsu-Sato) solution with $`\delta =2`$, which is referred to as TS2 in the present paper, does not have a horizon, by showing that the spacetime is inextendible across the segment I: $`\rho =0,\sigma <z<\sigma `$ in the Weyl coordinates, which corresponds to a horizon for odd $`\delta `$Tomimatsu.A&Sato1973a . Although they noticed that TS2 has quasi-regular directional singularities at P±: $`(\rho ,z)=(0,\pm \sigma )`$ in the Weyl coordinates, and that the spacetime has two-dimensional extensions along the $`z`$-axis across these points, they failed to recognize that these points are horizons. Subsequently, ErnstErnst.F1976  pointed out that these points are hypersurfaces, and investigated the behavior of geodesics crossing them by introducing a polar coordinate system around these points. However, he gave the incorrect conclusion that the corresponding surfaces are time-like and a particle which crossed one of the hypersurfaces can come out again through the same hypersurface. Later, Papadopoulos et alPapadopoulos.D&Stewart&Witten1981  recognized that the two points P± of the ZVW spacetime with $`\delta 2`$, which is a Weyl solution described by the Zipoy-Voorhees metricZipoy.D1966 ; Voorhees.B1970 ; Stephani.H&&2003B  and a static limit of the TS solution for integer values of $`\delta `$, are also hypersurfaces and tried to construct a coordinate system providing a smooth extension of the spacetime across these hypersurfaces. They found a two-dimensional extension along the $`z`$-axis, but were not able to find a coordinate system giving a four-dimensional extension.

One main reason for the lack of detailed studies may be the fact that the TS solutions and the Weyl solutions have pathological features including the existence of naked singularities. If the cosmic censorship hypothesisPenrose.R1969  were correct, such solutions would not describe the real world. However, the cosmic censorship is still a hypothesis, and there are now lots of counter examplesPapapetrou.A1985B ; Hollier.G1986 ; Singh.T1996A ; Gundlach.C2000A ; Harada.T&Iguchi&Nakao2002 . Of course, it is probable that the cosmic censorship hypothesis turns out to be correct in practice, since most of the analytic counter examples found so far are spherically symmetricWald.R1997A . Hence, in order to clarify whether the cosmic censorship is correct or not, it is inevitable to study non-spherically symmetric systemsLake.K1992 . For that purpose, knowledge on possible structures of spacetimes with naked singularities will give useful information. Such knowledge will be also important if the cosmic censorship turns out to be incorrect. From this point of view, in the present paper, we investigate in detail the global structure of the ZVW and TS spacetimes, with focus on singularities and horizons. Although we do not restrict the deformation parameter $`\delta `$ for the ZVW solution, we only consider the $`\delta =2`$ case for the TS solution, because it is hard to analyse the TS solutions with $`\delta 3`$, even by using symbolic computations by computers.

We will reveal lots of new features of the spacetimes described by these solutions concerning the singularity and causal structure. In particular, we will prove that the ZVW spacetimes with $`\delta =2,3`$ and $`\delta 4`$ have a degenerate Killing horizon at P±, as anticipated by Papadopoulos et al, and show that the horizon has a ring singularity at the equator. We further show that at least for $`\delta =2`$, the Tomimatsu-Sato spacetime has degenerate Killing horizons at P±, each of which is a sphere with conic singularity. We will also point out that the ZVW solution with $`\delta >0`$ ($`\delta 1`$) provides a rare example of non-trivial naked singularity with positive gravitational mass.

The paper is organized as follows. In section II, we first study in detail the structure of the singularities of the ZVW spacetimes, and then show that for $`\delta =2,3`$ and $`\delta 4`$, the directional quasi-regular singularities at P± are regular parts of a degenerate Killing horizon, by giving explicit regular extensions of the spacetime across these points. In section III, we extend the analysis to TS2. After a brief review of the basic features of this spacetime, we prove that the directional quasi-regular singularities at P± are degenerate Killing horizons. We also clarify the structure of singularity at the segment I on the $`z`$-axis. Section IV is devoted to summary and discussion. All symbolic computations in the present paper were done by Maple V.

## II Zipoy-Voorhees-Weyl Spacetime

The ZVW spacetimes are a two-parameter family of static, axisymmetric and asymptotically flat vacuum solutions to the Einstein equations described by the Zipoy-Voorhees metricZipoy.D1966 ; Voorhees.B1970 ; Stephani.H&&2003B , which include the flat and Schwarzschild solutions. In the canonical form for the static axisymmetric metric

$$ds^2=fdt^2+f^1[e^{2\gamma }(d\rho ^2+dz^2)+\rho ^2d\varphi ^2],$$
(1)
they are expressed asVoorhees.B1970

$$f=\left(\frac{x1}{x+1}\right)^\delta ,e^{2\gamma }=\left(\frac{x^21}{x^2y^2}\right)^{\delta ^2}.$$
(2)
Here, $`x`$ and $`y`$ are the prolate spheroidal coordinates related to the Weyl coordinates $`(\rho ,z)`$ by

$$\rho =\sigma \sqrt{(x^21)(1y^2)},z=\sigma xy.$$
(3)
In this prolate spheroidal coordinate system $`(t,\varphi ,x,y)`$, the ZVWmetric is written as

$$ds^2=fdt^2+R^2d\varphi ^2+\mathrm{\Sigma }^2\left(\frac{dx^2}{x^21}+\frac{dy^2}{1y^2}\right),$$
(4)
where

$`R^2=\sigma ^2\left({\displaystyle \frac{x+1}{x1}}\right)^{\delta 1}(x+1)^2(1y^2),`$ (5)
$`\mathrm{\Sigma }^2=\sigma ^2{\displaystyle \frac{(x+1)^{\delta ^2+\delta }}{(x1)^{\delta ^2+\delta }}}(x^2y^2)^{1\delta ^2}.`$ (6)
Here, note that the metric is invariant under the transformation $`\delta \delta ,xx`$. Further, the metric has a curvature singularity at some value of $`x`$, except for the case $`\delta =0`$, as we see below. Hence, we can assume that the asymptotically flat region corresponds to the range $`x>1`$ and $`1y1`$. In this asymptotic flat region, for the spherical coordinates $`(r,\theta ,\varphi )`$ defined by $`\sigma x=r\delta \sigma +\mathrm{O}(1/r)`$ and $`y=\mathrm{cos}\theta `$, the metric components in (4) are expanded as

$`f=1{\displaystyle \frac{2\sigma \delta }{r}}+\mathrm{O}\left({\displaystyle \frac{1}{r^3}}\right),`$ (7a)
$`R^2/\mathrm{sin}^2\theta =r^2+2\sigma \delta r+(4\delta ^21)\sigma ^2+\mathrm{O}\left({\displaystyle \frac{1}{r}}\right),`$ (7b)
$`\mathrm{\Sigma }^2=r^2+2\sigma \delta r+[4\delta ^21+(1\delta ^2)\mathrm{sin}^2\theta ]\sigma ^2+\mathrm{O}\left({\displaystyle \frac{1}{r}}\right).`$ (7c)
From these expressions, we see that the gravitational mass of the system is represented as $`M=\delta \sigma `$, and that the parameter $`\delta `$ taking a real number describes a quadrupole deformation of the spacetimeZipoy.D1966 ; Voorhees.B1970 : each $`f=`$constant surface has a prolate shape for $`|\delta |<1`$ and an oblate shape for $`|\delta |>1`$. In particular, $`\delta =0`$ corresponds to the flat metric, and $`\delta =\pm 1`$ to the Schwarzschild metric with positive or negative mass.

### II.1 Singularity structure

From (4), it is clear that the metric can have curvature singularity only at $`x=\pm 1`$ and $`x=\pm y`$. In this section, we first confirm that there is actually singularity at $`x=\pm 1`$ for $`\delta 0,1`$ and investigate the structure of the singularity there. In the present paper, we focus on singularities accessible from the asymptotically flat region $`x>1,1y1`$, and do not consider possible singularities at $`x=\pm y`$, which are separated from the asymptotically flat region by naked singularities or horizons.

For a static and axisymmetric metric, $`\mathrm{\Psi }_1=\mathrm{\Psi }_3=0`$ and all scalar polynomials of the Weyl curvature can be written as polynomials of $`\mathrm{\Psi }_2`$ and $`\mathrm{\Psi }_0\mathrm{\Psi }_4`$, where $`\mathrm{\Psi }_0\mathrm{\Psi }_4`$ are the Newman-Penrose components of the Weyl curvature with respect to the complex null tetrad $`(k,l,m,\overline{m})`$ defined by

$$k_\mu dx^\mu =\frac{\sqrt{f}dt+Rd\varphi }{\sqrt{2}},l_\mu dx^\mu =\frac{\sqrt{f}dtRd\varphi }{\sqrt{2}},m_\mu dx^\mu =\frac{e^\gamma (d\rho idz)}{\sqrt{2f}}.$$
(8)
For the metric (4), these curvature invariants are given by

$`\mathrm{\Psi }_2={\displaystyle \frac{\delta (x\delta )(x^2y^2)^{\delta ^21}}{2(x1)^{\delta ^2\delta +1}(x+1)^{\delta ^2+\delta +1}}},`$ (9b)
$`\mathrm{\Psi }_0\mathrm{\Psi }_4=`$ $`{\displaystyle \frac{\delta ^2(x^2y^2)^{2\delta ^23}}{4[(x1)^{\delta ^2\delta +1}(x+1)^{\delta ^2+\delta +1}]^2}}\times [9(x1)^2(x^2y^2)`$
$`+(\delta 1)\{[12(\delta +1)x^26(2\delta ^2+2\delta +3)x+(\delta +1)(4\delta ^2+5)](1y^2)`$
$`9(2x\delta 1)(x^21)\}].`$
From these expressions, we immediately see that the curvature invariants diverge on the open segment I: $`x=1,1<y<1`$, if $`\delta 0,1`$. This segment corresponds to the segment $`\rho =0,\sigma <z<\sigma `$ in the Weyl coordinates, which is depicted by a thick line in Fig.2. It is also easy to see that one can reach any point on this segment by a curve with a finite proper length in the $`xy`$ plane. As is expected from this, this segment is a naked singularity. In fact, we can show that in the case $`L=0`$ and $`ϵ=0`$ (null geodesics), the geodesic equation (62d) has a regular solution expressed by a power series of $`x1`$ as

$$y=y_0+\frac{\delta +1}{\delta 1}y_0(x1)+\mathrm{},$$
(10)
for any $`y_0`$ in the open interval $`(1,1)`$, provided that $`\delta 1`$. From (62c), we see that the null geodesic represented by this solution reaches the point $`(x,y)=(1,y_0)`$ in a finite affine parameter. This implies that each point in the open segment $`\rho =0,\sigma <z<\sigma `$ is a naked curvature singularity in the standard sense. By looking for a solution of the type $`y=y_0+a(x1)^p+\mathrm{}`$ ($`p>0`$) to (62d), we can also show that for $`\delta <0`$ or $`\delta >2`$ there exists no other solution that reaches the point $`(1,y_0)`$ for $`y_0^2<1`$, while for $`0<\delta <1`$ or $`1<\delta 2`$ there exist one parameter family of such solutions with $`p=(1+2\delta \delta ^2)/2`$. Hence, for the latter case, one can reach every point on the open segment I from a point with $`x>1`$ by a null geodesic. Fig. 3 gives two examples of such a family of null geodesics for the case $`\delta =2`$.

The curvature invariants also exhibit singular behavior at the two points P±: $`(\rho ,z)=(0,\pm \sigma )`$. To see this, we introduce the new coordinate system $`(X,Y)`$ defined by

$$X^2=\frac{1y^2}{x^21},Y=\frac{y}{x}.$$
(11)
The coordinates $`x`$ and $`y`$ are expressed in terms of $`X`$ and $`Y`$ as

$$x=\sqrt{\frac{X^2+1}{X^2+Y^2}},y=Y\sqrt{\frac{X^2+1}{X^2+Y^2}}.$$
(12)
The region $`x1,|y|1`$ is mapped to the region $`X0,|Y|1`$, and $`X=`$const curves and $`Y=`$const curves are represented by circles with centers on the $`\rho `$-axis and on the $`z`$-axis, respectively, as shown in Fig.2.

Inserting these expressions into (9), we find that the curvature invariants behave in the limit $`Y^21`$ with fixed $`X`$ as

$`\mathrm{\Psi }_2{\displaystyle \frac{\delta (1\delta )}{2^{2\delta +1}}}\left(1+X^2\right)^{\delta ^2\delta +1}(1Y^2)^{\delta 2},`$ (13a)
$`\mathrm{\Psi }_0\mathrm{\Psi }_4{\displaystyle \frac{\delta ^2(\delta 1)^2}{2^{4\delta +4}}}\left(1+X^2\right)^{2\delta ^22\delta +1}\left(9+(2\delta 1)^2X^2\right)\left(1Y^2\right)^{2(\delta 2)}.`$ (13b)
From this, we see that the two points P±(i.e., $`Y=\pm 1`$) are curvature singularities for $`\delta <2`$ ($`\delta 0,1`$). In this case, one can reach each of these points by a spatial curve as well as by a null geodesic along the $`z`$-axis in the $`\rho z`$ plane with a finite affine length $`\mathrm{}`$, because it follows from (62c) that the affine parameter is proportional to the $`x`$-coordinate for such geodesics. Hence, they are naked.

In contrast, for $`\delta >2`$, the curvatures vanish as one approaches P±. Hence, they are not curvature singularity. Furthermore, there exists no spatial curve that reaches these points within a finite proper length, as was pointed out by Papadopoulos et alPapadopoulos.D&Stewart&Witten1981 . However, since null geodesics along the $`z`$-axis reach these points in finite affine parameters, they are quasi-regular singularity for the region $`x>1`$Papadopoulos.D&Stewart&Witten1981 . In the next subsection, we will show that these points are in reality Killing horizons at least for $`\delta =3`$ and $`\delta 4`$.

The situation for $`\delta =2`$ is similar but more subtle. The curvature invariants approach non-vanishing values in the limit $`Y\pm 1`$, and these limits depend on $`X`$, i.e., on the direction one approaches the points P±. They are again at spatial infinity, but can be reached by null geodesics with finite affine parameters. Hence, they are apparently directional quasi-regular singularitiesPapadopoulos.D&Stewart&Witten1981 . We will also show that these points are Killing horizons.

Next, let us examine geometrical shapes of these singularities. One simple quantity characterizing shapes is the circumferential radius $`R`$ of circles generated by the Killing vector $`\eta =_\varphi `$, given by (5). The asymptotic behavior of this quantity near the singularity is summarized in Table 1. This asymptotic behavior suggests that the singularity is point-like or string like for $`\delta <1`$ ($`\delta 0`$) and is ring-like for $`\delta >1`$.

To confirm this, we estimate the length of the singularity in the $`z`$-direction. There are various candidates for the definition of such a length. One natural candidate is the length of $`x=`$const curves,

$`L_x`$ $`:={\displaystyle _1^1}{\displaystyle \frac{\mathrm{\Sigma }dy}{\sqrt{1y^2}}}`$ (14)
$`=\sigma (x+1)^{\delta (\delta +1)/2}(x1)^{\delta (\delta 1)/2}{\displaystyle _0^{\pi /2}}𝑑\theta (x^21+\mathrm{sin}^2\theta )^{(1\delta ^2)/2}.`$
The second candidate is the length of $`X=`$const curves, which is expressed in the limit $`X+\mathrm{}`$ as

$$L_X2^\delta \sigma X^{\delta (1\delta )}_1^1\frac{dY}{(1Y^2)^{\delta /2}}.$$
(15)
Here, the integral on the right-hand side diverges for $`\delta 2`$, hence $`L_X`$ can have a finite limit only for $`\delta <2`$. The asymptotic behavior of these quantities in the limits $`x1`$ and $`X+\mathrm{}`$ are summarized in Table 1. Comparison of this behavior of $`L_x`$ and $`L_X`$ with that of $`R`$ shows that the singularity is point-like for $`\delta <0`$, string-like for $`0<\delta <1`$, and ring-like for $`1<\delta <2`$. In contrast, these quantities do not seem to describe the length of singularities for $`\delta 2`$, because the two points P± are horizons and not singularity, as we will show later. For these cases, we have to measure the length of the singularity by the limit length of curves which are tangent to the segment I at $`(\rho ,z)=(0,\pm \sigma )`$. Therefore, we define the length of I by the $`a\mathrm{}`$ limit of the length $`L_n`$ of curves with $`X(1Y^2)^n=a`$(constant) for $`n>0`$. In this limit, $`L_n`$ is expressed as

$$L_n2^\delta \sigma a^{\delta (1\delta )}_1^1𝑑Y(1Y^2)^{n\delta (\delta 1)\delta /2}.$$
(16)
The integral on the right-hand side of this equation converges for $`2\delta 2n\delta (1\delta )>0`$. For any $`\delta `$, there exists $`n>0`$ that satisfies this inequality, and for such $`n`$, the asymptotic behavior of $`L_n`$ is simply determined by the sign of $`\delta (1\delta )`$, as summarized in Table 1. This behavior implies that the singularity is ring-like for $`\delta 2`$.

To summarize, the segment singularity I is geometrically point-like for $`\delta <0`$, string-like for $`0<\delta <1`$ and ring-like for $`\delta >1`$. This feature is consistent with the behavior (7) of the metric at spatial infinity for $`\delta >0`$, although such a good correspondence does not exist for $`\delta <0`$.

Here, note that this characterization of the singularity is based on the purely local geometry and may not coincide with its physical appearance. For example, if an observer in the region $`x>1`$ looks at the singularity by light rays, he will observe a cylinder-like object for the case $`1<\delta 2`$, because light rays starting from different points on I can reach the same point in $`x>1`$. For the other values of $`\delta `$, the geometrical shape and the physical appearance coincide.

Finally, we examine the Komar mass of the singularity. The Komar mass is in general defined in terms of the time-like Killing vector $`\xi =_t`$ asKomar.A1959 ; Komar.A1962 ; Heusler.M1996B

$$M=\frac{1}{16\pi }_\mathrm{\Sigma }ϵ_{\mu \nu \lambda \sigma }^\mu \xi ^\nu dx^\lambda dx^\sigma .$$
(17)
Here, the orientation of the spatial two-surface $`\mathrm{\Sigma }`$ should be chosen so that for a future-directed time-like vector $`u`$ Can outward normal $`n`$ to $`\mathrm{\Sigma }`$ and a basis $`(v,w)`$ giving the orientation of $`\mathrm{\Sigma }`$, $`(u,n,v,w)`$ gives the positive orientation of the spacetime.

For the ZVW metric, in the case in which $`\mathrm{\Sigma }`$ is an axisymmetric surface determined by a curve $`\gamma `$ in the $`xy`$ plane, this expression reduces to

$$M=\frac{\delta \sigma }{2}_\gamma 𝑑y,$$
(18)
where the orientation of $`\gamma `$ should be taken so that for the tangent vector $`V`$ of $`\gamma `$ and a normal $`N`$ to $`\gamma `$ directed outward from the region enclosed by $`\gamma `$, $`(N,V)`$ has the same orientation as that of $`(_x,_y)`$. It is clear from this expression that each portion of the segment singularity has a positive mass and the total mass carried by the singularity coincides with the total mass of the system, $`M=\delta \sigma `$. This result is quite interesting because most of the known examples of naked singularities have negative or zero massSingh.T1996A ; Gergely.L2002 . For example, central shell-focusing singularities formed by spherically symmetric gravitational collapse always have non-positive local gravitational mass if they are locally nakedLake.K1992 . One exception is a shell-crossing singularity, which can have a positive mass. However, this type of singularity is weak in the sense that the spacetime has a $`C^0`$ extension across the singularityNewman.R1986 . Another counter example is the time-like singularity of the Janis-Newman-Winicour solutionJanis.A&Newman&Winicour1968 ; Virbhadra.K1997  for the spherically symmetric Einstein-Massless-Scalar (EMS) system. Virbhadra showed that this singularity has a positive mass if one defines the mass by the spatial integration of the Einstein energy-momentum complexVirbhadra.K1999 . This mass is identical to the Komar mass, whose value is independent of the location of the integration surface $`\mathrm{\Sigma }`$ for a static solution of the spherically symmetric EMS system. However, it does not coincide with the local gravitational mass used in the mass inflation argumentPoisson.E&Israel1990 . We can show that the latter becomes negatively infinite for the singularity of the Janis-Newman-Winicour solution. This discrepancy seems to come from the difference in the matter contribution to the mass. In contrast, because the ZVW solution discussed in the present paper is a vacuum solution, this kind of ambiguity associated with the matter contribution does not exist in the definition of the mass for the singularity of the ZVW spacetime. Thus, it seems to provide a clearer non-trivial example of naked singularities with positive gravitational mass.

### II.2 Extension across P±

In this subsection, we will demonstrate that the two points P± are Killing horizons for $`\delta =2,3`$ and $`\delta 4`$, as anticipated by Papadopoulos et alPapadopoulos.D&Stewart&Witten1981 , by constructing smooth extensions of the spacetime across these points.

The coordinate system around the points P± ($`(X,Y)=(0,\pm 1)`$) that we adopt to extend the spacetime is

$`S=2^\delta \sigma {\displaystyle \frac{X}{\sqrt{1+X^2}}}\left({\displaystyle \frac{1Y^2}{Y^2}}\right)^{(2\delta )/2},`$ (19a)
$`W={\displaystyle \frac{1}{4}}(1Y^2)(1+X^2)^{\delta 2}.`$ (19b)
For $`\delta =2`$, $`S`$ and $`W`$ depend only on $`X`$ and $`Y`$, respectively, and these two coordinate systems are essentially equivalent. However, for $`\delta >2`$, each $`S=`$const curve is tangential to the $`z`$-axis in the $`\rho z`$ plane, as is shown in Fig. 4. Hence, the $`(S,W)`$ coordinates expand cusp regions at P± in the $`(X,Y)`$ coordinates. Such a singular coordinate system is required because $`R`$ diverges as $`Y`$ tends $`\pm 1`$ along $`X=`$const curves for $`\delta >2`$, as is shown in Table 1.

In this new coordinate system, the ZVW metric is expressed as

$$ds^2=aW^\delta \left(dt^2h^2\frac{dW^2}{W^{2\delta }}\right)+bS^2d\varphi ^2+cdS^2,$$
(20)
with

$`a=\left({\displaystyle \frac{2}{x+1}}\right)^{2\delta }(X^2+Y^2)^\delta (1+X^2)^{\delta (2\delta )},`$ (21a)
$`b=\left({\displaystyle \frac{x+1}{2}}\right)^{2\delta }Y^{2(2\delta )}(1+X^2)(X^2+Y^2)^{\delta 2},`$ (21b)
$`c=\left({\displaystyle \frac{x+1}{2}}\right)^{2\delta }{\displaystyle \frac{(1+X^2)^{3\delta ^2}(X^2+Y^2)^{\delta 2}}{Y^{2(\delta 3)}[Y^2+(\delta 2)^2X^2]}},`$ (21c)
$`h^2=4\sigma ^2\left({\displaystyle \frac{x+1}{2}}\right)^{4\delta }{\displaystyle \frac{(1+X^2)^{\delta ^26\delta +6}(X^2+Y^2)^{2(\delta 1)}}{Y^2+(\delta 2)^2X^2}}.`$ (21d)
Here, the dependence of $`x`$, $`X`$ and $`Y`$ on $`S`$ and $`W`$ is implicitly determined by the relations

$`x=\sqrt{{\displaystyle \frac{X^2+1}{X^2+Y^2}}},`$ (22a)
$`Y^2=1{\displaystyle \frac{4W}{(1+X^2)^{\delta 2}}},`$ (22b)
$`X^2={\displaystyle \frac{S^2}{16\sigma ^2}}{\displaystyle \frac{(1+X^2)W^{\delta 2}}{[(1+X^2)^{\delta 2}4W]^{\delta 2}}}.`$ (22c)
This coordinate system covers only a neighbourhood of P+ or P- in general, and $`P_\pm `$ corresponds to $`W=0`$. For $`\delta =2`$, $`a,b`$ and $`c`$ have non-vanishing finite limits at $`W=0`$, which are expressed as regular functions of $`X`$, or equivalently of $`S`$. On the other hand, for $`\delta >2`$, in the limit $`W0`$ with $`S`$ fixed, $`X^2`$ behaves as

$$X^2=\frac{S^2}{16\sigma ^2}W^{\delta 2}\left(1+\mathrm{O}\left(W\right)\right).$$
(23)
Hence, $`a,b`$ and $`c`$ approach unity in the same limit.

Here, note that if $`\delta `$ is an integer equal to or greater than 2, $`a,b`$ and $`c`$ have unique regular analytic extensions to the region $`W<0`$. A $`C^2`$ extension is possible even for non-integer $`\delta `$ by replacing $`W^\delta `$ in (22c) by $`|W|^\delta `$, if $`\delta 4`$.

In the extension of the spacetime, a crucial point is the asymptotic behavior of $`h`$ at $`W=0`$. First, we examine it for $`\delta >2`$. Let us define the function $`h_0(W)`$ by

$$h_0(W):=h(S=0,W)=\frac{2\sigma }{(14W)^{3/2}}\left[\frac{1}{2}2W+W^2+\frac{1}{2}(12W)\sqrt{14W}\right]^{\delta /2}.$$
(24)
Then, we can show that $`h^2h_0^2`$ is expressed as

$`h(S,W)^2h_0(W)^2`$ $`=\left[16\delta (\delta 1)W^2+\mathrm{O}\left(W^3\right)\right]\sigma ^2X^2`$ (25)
$`2(\delta 1)(\delta 2)^2(\delta 3)\sigma ^2X^4+\mathrm{O}(WX^4,X^6)`$
$`=\delta (\delta 1)W^\delta S^2{\displaystyle \frac{(\delta 1)(\delta 2)^2(\delta 3)}{128\sigma ^2}}W^{2\delta 4}S^4`$
$`+\mathrm{O}(W^{\delta +1},W^{2\delta 3}).`$
From this, it follows that for $`\delta =3`$ or $`\delta 4`$, $`h^2`$ is expressed in terms of a function $`k(S,W)`$ that is continuous at $`W=0`$ as

$$h^2(S,W)=h_0^2(W)+W^\delta k(S,W).$$
(26)
Next, for $`\delta =2`$, $`h^2`$ can be expanded with respect to $`W`$ as

$$h^2=4\sigma ^2\left[1+4W+\left(20\frac{S^2}{2\sigma ^2}+\frac{S^4}{64\sigma ^4}\right)W^2+\mathrm{O}\left(W^3\right)\right].$$
(27)
Hence, the expression (26) also holds in this case with a function $`k`$ that is smooth at $`W=0`$.

Now, we can construct a regular extension across P± by introducing the retarded/advanced-type coordinate

$$u_\pm :=t\pm ^Wh_0(W)\frac{dW}{W^\delta }.$$
(28)
In term of the coordinate system $`(t,\varphi ,u_\pm ,S)`$, the ZVW metric is written

$$ds^2=aW^\delta du_\pm ^2\pm 2ah_0dWdu_\pm +akdW^2+bS^2d\varphi ^2+cdS^2.$$
(29)
It is clear from the comment above on the behavior of $`a`$, $`b`$, $`c`$ and $`k`$ that this metric can be regularly extended to the region $`W<0`$ for $`\delta =2,3`$ and $`\delta 4`$. Here, it is understood that $`W^\delta `$ is replaced by $`|W|^\delta `$ for non-integer $`\delta `$, as mentioned before. This extension connects the regular solutions in two regions with $`x>1`$ and $`1<x<1`$ as shown in Fig. 5, and is analytic and unique in the case $`\delta `$ is an integer (the extension path in the $`xy`$ plane for a non-integer $`\delta `$ is the same as that in the case $`\delta `$ is an even integer). In contrast, for non-integer $`\delta `$ such that $`\delta >4`$, the extension is at least $`C^{[\delta ]2}`$, where $`[\delta ]`$ represents the maximum integer less than or equal to $`\delta `$, but the $`([\delta ]1)`$-th derivatives of the metric are divergent at $`W=0`$ in general.

Thus, we have shown that the ZVW spacetime with $`\delta =2,3`$ or $`\delta 4`$ has Killing horizons at P±: $`(\rho ,z)=(0,\pm \sigma )`$. These two horizons share the ring singularity corresponding to the open segment I as the common boundary. As a whole, they can be regarded as a single spherical horizon with a ring singularity at the equator whose circumference is infinite, as is illustrated in Fig. 6. Since the $`W`$-derivative of $`g_{tt}`$ vanishes at $`W=0`$, this horizon is degenerate in its regular part in the sense that the surface gravity vanishes.

Here, note that the four-dimensional extension constructed here coincides with the two-dimensional extension along the $`z`$-axis constructed by Papadopoulos et alPapadopoulos.D&Stewart&Witten1981 . In particular, the conformal diagram of this two-dimensional sector of the extended spacetime is given by Fig. 3 for even $`\delta `$ and by Fig. 4 for odd $`\delta `$ in their paper. We can also easily draw a conformal diagram of the two-dimensional sector corresponding to the equatorial plane. The result is just the same as the diagram for a negative-mass Schwarzschild black hole, for any $`\delta 1`$. Hence, the singularity at the segment I is time-like.

Finally, we comment on the area of the horizon. From (29), its area $`A`$ for $`\delta >2`$ is infinite:

$$A=\underset{W0}{lim}2\times 2\pi _0^{\mathrm{}}\sqrt{b}S𝑑S=+\mathrm{}.$$
(30)
In contrast, for $`\delta =2`$, it is given by the finite value

$$A=\underset{W0}{lim}2\times 2\pi _0^{4\sigma }\sqrt{bc}S𝑑S=32\pi \sigma ^2.$$
(31)
This horizon area for the $`\delta =2`$ case is half of that for the Schwarzschild black hole with the same mass $`M=2\sigma `$ and equal to twice of that for the Schwarzschild black hole with mass $`M/2`$. To be precise, we can show that the Schwarzschild solution with mass $`M=2\sigma `$ can be represented by the Israel-Kahn solutionIsrael.W&Kahn1964  with two black holes of the same mass whose centers are separated by $`2\sigma `$ in the Weyl coordinates, and when these centers come closer, the total horizon area decreases monotonically. The above area of the $`\delta =2`$ ZVW solution is obtained in the limit the two centers coincide.

## III Tomimatsu-Sato spacetime

In this section, we extend the analysis on the ZVW solution to the Tomimatsu-Sato solution. In the present paper, we only consider TS2, i.e., the $`\delta =2`$ case.

The metric of TS2 is expressed in the prolate spheroidal coordinates asTomimatsu.A&Sato1972

$`ds^2`$ $`=fdt^2+2f\omega dtd\varphi +R^2d\varphi ^2+\mathrm{\Sigma }^2\left({\displaystyle \frac{dx^2}{x^21}}+{\displaystyle \frac{dy^2}{1y^2}}\right),`$ (32)
with

$$f=\frac{A}{B},f\omega =\frac{4\sigma q}{p}\frac{(1y^2)C}{B},\mathrm{\Sigma }^2=\frac{B}{p^4(x^2y^2)^3},R^2=\frac{\sigma ^2}{p^2}\frac{(1y^2)D}{B}.$$
(33)
Here, $`A,B`$, $`C`$ and $`D`$ are polynomials given by

$`A=`$ $`(1y^2)^4g(Z)g(Z),`$ (34a)
$`B=`$ $`(g(x)+q^2y^4)^2+4q^2y^2\left[px^3+1(px+1)y^2\right]^2,`$ (34d)
$`C=`$ $`q^2(px+1)y^4(y^2+3)+\left\{2q^2(px^3+1)+(px+1)g(x)\right\}y^2`$
$`(2px^3px+1)g(x),`$
$`D=`$ $`p^2q^2(x^21)y^8+4q^2(px+1)(p^3x^3+3p^2x^2p^3x+4px3p^2+4)y^6`$
$`2q^2(3p^4x^6+4p^3x^53p^4x^4+8p^3x^3+37p^2x^212p^3x+48px13p^2+24)y^4`$
$`+4q^2(p^4x^8p^4x^6+4p^3x^54p^3x^3+15p^2x^2+24px3p^2+12)y^2`$
$`+g(x)(p^4x^6+6p^3x^5p^4x^4+16p^2x^412p^3x^3+32px^3`$
$`+15p^2x^2+6p^3x15p^2+16),`$
where

$$g(x)=p^2x^4+2px^32px1,Z^2:=\frac{x^2y^2}{1y^2}.$$
(35)
This solution is asymptotically flat, and the parameters $`\sigma ,p`$ and $`q`$ are related to its mass $`M`$ and angular momentum $`J`$ as

$$p^2+q^2=1,M=\frac{2}{p}\sigma ,J=M^2q.$$
(36)
In particular, in the non-rotating case $`q=0`$, this solution reduces to the $`\delta =2`$ ZVW solution, which is referred to as ZVW2 from this point. In the present paper, we only consider the case $`0<p,q<1`$.

### III.1 Basic properties

Although TS2 reduces to ZVW2 in the limit $`q=0`$, the former solution has quite different features compared with the latterTomimatsu.A&Sato1973 ; Tomimatsu.A&Sato1973a ; Gibbons.G&Russell-Clark1973 .

#### III.1.1 Ergo region

$`f=g_{tt}=\xi \xi `$ is not positive definite in the asymptotically flat region $`x>1,1y1`$, and the so-called ergo region appearsGibbons.G&Russell-Clark1973 . From (34a) and (34d), $`f`$ becomes negative in the region $`A<0`$. Since $`g(x)=0`$ has two single real roots $`x_1`$ and $`x_0`$ with $`1<x_0<x_1`$ for $`0<p<1`$, the ergo region in $`x1`$ is given by

$$x_0^2(x_0^21)y^2<x^2<x_1^2(x_1^21)y^2.$$
(37)
Note that since $`Z`$ in (35) is expressed in terms of the $`(X,Y)`$ coordinates in (11) as $`Z^2=1+1/X^2`$, the boundaries of this ergo region, i.e., the infinite redshift surfaces, are represented as $`X=`$const (see Fig. 7).

#### III.1.2 Causality violation

Since $`D`$ is not positive definite in $`x1`$, the norm $`R^2`$ of the rotation Killing vector $`\eta =_\varphi `$ becomes negative in some regionGibbons.G&Russell-Clark1973 . Each $`S^1`$ orbit in this region becomes a closed time-like curve. This region contains the segment I: $`\rho =0,\sigma <z<\sigma `$, where $`\rho `$ and $`z`$ are the Weyl-type coordinates related to $`x`$ and $`y`$ by (3), and contacts with the inner boundary of the ergo region from (34d), as shown in Fig. 7. The appearance of this causality violating region in the domain of outer communication is the most exotic feature of TS2.

#### III.1.3 Singularity

From (32), TS2 may have curvature singularity at points where $`B=0`$, in addition to the points $`x=\pm 1`$ and $`x=\pm y`$ in the Weyl case. This is confirmed by looking at the behavior of the Weyl curvature. As in the case of the ZVW solutions, $`\mathrm{\Psi }_1=\mathrm{\Psi }_3=0`$ with respect to the null complex tetrad (8) with $`dt`$ replaced by $`dt\omega d\varphi `$, and all polynomial curvature invariants are expressed as polynomials of $`\mathrm{\Psi }_2`$ and $`\mathrm{\Psi }_0\mathrm{\Psi }_4`$. By symbolic computations, we can show that they have the forms

$$\mathrm{\Psi }_2=\frac{P_1}{B^3},\mathrm{\Psi }_0\mathrm{\Psi }_4=\frac{P_2}{B^6},$$
(38)
where $`P_1`$ and $`P_2`$ are polynomials of $`x`$ and $`y`$ that have no non-trivial common divisor with $`B`$. Since the explicit expressions for these curvature invariants are quite long, we only give their values at the equatorial plane $`y=0`$:

$`\mathrm{\Psi }_2(x,y=0)={\displaystyle \frac{p^4x^6(23px+px^3)}{g(x)^3}},`$ (39a)
$`(\mathrm{\Psi }_0\mathrm{\Psi }_4)(x,y=0)={\displaystyle \frac{9p^8x^{10}(x^8p^22p^2x^6+4x^5p+5p^2x^44px^3+8px+4)}{g(x)^6}}.`$ (39b)
From these expressions, it follows that the curvature invariants diverge at the zero of $`g(x)`$, $`x=x_0`$, which is just the intersection of the inner boundary of the ergo region with the equatorial plane. It is easy to see from (34d) that $`B`$ vanishes only at $`(x,y)=(x_0,0)`$ for $`x1`$ and $`0<p<1`$, except for the two points P±: $`(x,y)=(1,\pm 1)`$. Hence, in contrast to the ZVW spacetime, there exists no curvature singularity on the open segment I, and instead there appears the famous ring singularityTomimatsu.A&Sato1972 , which has an infinite circumference, because $`x=x_0`$ is a single root of $`D`$ and a double root of $`B`$ on the equatorial plane $`y=0`$ from (34d) and (34d). This ring singularity reduces to the ring-like singularity at the segment I of ZVW2 in the limit $`q0`$.

$`B`$ also vanishes at P±. However, we can show by symbolic computations that the curvature invariants have finite limits at these points. These limits are direction dependent as in the case of the ZVW spacetime. For example, their values in the cases in which one approaches P± from $`|z|>\sigma `$ and from $`|z|<\sigma `$ along the $`z`$-axis are given by

$`X=0,Y=\pm 1:`$ $`\mathrm{\Psi }_2={\displaystyle \frac{p(p\pm iq)}{8(1+p)}},\mathrm{\Psi }_0\mathrm{\Psi }_4={\displaystyle \frac{9p^2(piq)^2}{64(p+1)^2}},`$ (40a)
$`X=\mathrm{},Y=\pm 1:`$ $`\mathrm{\Psi }_2={\displaystyle \frac{p^4(q\pm ip)}{8q^3(1+p)}},\mathrm{\Psi }_0\mathrm{\Psi }_4={\displaystyle \frac{9p^8(piq)^2}{64q^6(p+1)^2}}.`$ (40b)
Thus, P± are quasi-regular singularitiesGibbons.G&Russell-Clark1973 ; Ernst.F1976 . In the next subsection, we will show that these two points are disconnected Killing horizons, by constructing regular extensions of the spacetime across these points.

All metric components of (32) as well as the polynomial curvature invariants are regular and finite on the open segment I, which is expressed as $`x=1,1<y<1`$ in the $`(x,y)`$ coordinates, in contrast to the case of ZVW2. However, the spacetime turns out to be singular on this segmentGibbons.G&Russell-Clark1973 . To see this, let us rewrite the metric (32) as

$$ds^2=f(dt\omega d\varphi )^2+\mathrm{\Sigma }^2\left(d\chi ^2+\mathrm{}^2\mathrm{sinh}^2\chi d\varphi ^2\right)+\mathrm{\Sigma }^2\frac{dy^2}{1y^2},$$
(41)
where $`\chi `$ is related to $`x`$ by $`x=\mathrm{cosh}\chi `$, and $`\mathrm{}`$ is defined by

$$\mathrm{}^2=\frac{\sigma ^2(1y^2)}{f\mathrm{\Sigma }^2}.$$
(42)
In the neighbourhood of the segment I, $`\omega `$ behaves as

$`\omega ={\displaystyle \frac{4\sigma (p+1)}{qp}}+\omega _1(x,y)\mathrm{sinh}^2\chi ;`$
$`\omega _1={\displaystyle \frac{4\sigma (p+1)[(p1)(y^4+1)+2(3p+1)y^2]}{q^3(1y^2)^2}}+\mathrm{O}\left(\chi ^2\right).`$ (43)
Hence, if we introduce the new time coordinate $`\stackrel{~}{t}`$ by

$$\stackrel{~}{t}=t\omega _0\varphi ;\omega _0=\frac{4\sigma (p+1)}{qp},$$
(44)
(41) can be written as

$$ds^2=f(d\stackrel{~}{t}\omega _1\mathrm{sinh}^2\chi d\varphi )^2+\mathrm{\Sigma }^2\left(d\chi ^2+\mathrm{}^2\mathrm{sinh}^2\chi d\varphi ^2\right)+\mathrm{\Sigma }^2\frac{dy^2}{1y^2}.$$
(45)
From this, it follows that the spacetime is regular on the segment I, if $`\mathrm{}`$ is unity on I. However, since $`\mathrm{}`$ behaves around I as

$$\mathrm{}=\frac{p^2}{q^2}+\frac{p^2(p^2+3)}{q^4(1y^2)}\mathrm{sinh}^2\chi +\mathrm{O}\left(\chi ^4\right),$$
(46)
the spacetime is not regular on I unless $`p=q=1/\sqrt{2}`$ and has a conical singularity there. The deficit angle is positive for $`q^2>p^2`$ and negative for $`q^2<p^2`$.

This is not the whole story. There is another subtlety concerning the spacetime structure around IGibbons.G&Russell-Clark1973 . Since $`\varphi `$ is an angle coordinate and defined modulo $`2\pi `$, the new time coordinate $`\stackrel{~}{t}`$ is not a single-valued function, and each $`\stackrel{~}{t}=`$const surface has a helical structure in the original coordinates, as illustrated in Fig. 8. These helical surfaces cannot be extended regularly to the axis $`\chi =0`$. Hence, $`\varphi `$ cannot be regarded as an angle coordinate on these $`\stackrel{~}{t}=`$const surfaces, but rather should be regarded as a parameter of helical orbits of the Killing vector $`_\varphi +\omega _0_t`$. In this sense, the deficit angle interpretation above does not have a global meaning. This exotic structure arises because the two Killing vectors $`_t`$ and $`_\varphi `$ become parallel on I.

#### III.1.4 Komar mass

For the metric (32), the Komar mass (17) contained in a 3-surface that is bounded by a closed 2-surface determined by a curve $`\gamma `$ in the $`xy`$ plane is expressed as

$$M=_\gamma \left(M_xdx+M_ydy\right),$$
(47)
where

$$(1y^2)M_ydx(x^21)M_xdy=\frac{1}{4\sigma }\left[R^2df+f\omega d(f\omega )\right].$$
(48)
We can show that for any curve $`\gamma `$ whose end points are on the regular part of the $`z`$-axis, $`z>\sigma `$ and $`z<\sigma `$, $`M`$ is equal to the total mass $`M=2\sigma /p`$, if $`\gamma `$ does not pass the ring singularity. This implies that the ring singularity has zero gravitational mass. In contrast, from the asymptotic behavior of $`M_x`$ and $`M_y`$ around $`x=1`$,

$`M_x={\displaystyle \frac{4\sigma y[(1p)(y^6y^4)+(p+7)y^2+9p]}{[(1p)y^4+2(p+3)y^2+1p]^2}}+\mathrm{O}\left(x1\right),`$ (49a)
$`M_y={\displaystyle \frac{4\sigma (1y^2)[(1p)(y^4+1)2(3p+1)y^2]}{p[(1p)y^4+2(p+3)y^2+1p]^2}}+\mathrm{O}\left(x1\right),`$ (49b)
the main contribution to the integral on the right-hand side of (47) comes from the segment I for a curve close to I (see Fig. 9). This seems to indicate that the gravitational mass is carried by the singularity on I. However, this interpretation may not be valid, because a 2-surface corresponding to a curve close to the segment I is time-like. Furthermore, $`M_y`$ is not positive definite, as shown in Fig. 9.

### III.2 Extension across P±

In order to construct an extension of TS2, we utilize the coordinates $`(X,W)`$, where $`X`$ is the coordinate defined in (11) and $`W`$ is defined in terms of the $`Y`$ coordinate there by

$$W=\frac{1}{4}(1Y^2).$$
(50)
Then, we can show that the metric is written in terms of functions $`a(X,W)`$, $`b(X,W)`$, $`c(X,W)`$, $`h(X,W)`$ and $`\mathrm{\Omega }(X,W)`$ that are regular analytic in a neighbourhood of $`W=0`$ i$`Y=\pm 1`$) as

$$ds^2=aW^2\left(dt^2h^2\frac{dW^2}{W^4}\right)+bX^2(d\varphi \mathrm{\Omega }dt)^2+cdX^2.$$
(51)
The asymptotic behavior of these functions at $`W=0`$ are given by

$`a={\displaystyle \frac{2(p^2+q^2X^4)}{(p+1)(1+X^2)^2}}+{\displaystyle \frac{8(p^2q^2X^6)}{(p+1)(1+X^2)^3}}W+\mathrm{O}\left(W^2\right),`$ (52a)
$`b={\displaystyle \frac{8\sigma ^2(p+1)}{p^2+q^2X^4}}+\mathrm{O}\left(W\right),`$ (52b)
$`c={\displaystyle \frac{8\sigma ^2(p+1)(p^2+q^2X^4)}{p^4(1+X^2)^4}}+\mathrm{O}\left(W\right),`$ (52c)
$`h^2={\displaystyle \frac{(p+1)^2}{p^4}}\left\{1+4W+4(p+4)W^2\right\}`$
$`{\displaystyle \frac{4(p+1)^2\left\{2p^3+p^2(2p1)X^2+q^2X^4\right\}X^2}{p^4(1+X^2)^2}}W^2+\mathrm{O}\left(W^3\right),`$ (52d)
$`\sigma \mathrm{\Omega }={\displaystyle \frac{pq}{p+1}}W{\displaystyle \frac{q(p^2+q^2X^4)X^2}{p(p+1)(1+X^2)^2}}W^2+\mathrm{O}\left(W^3\right).`$ (52e)
From this behavior, if we put

$$h_0(W)=h(X=0,W),\mathrm{\Omega }_0(W)=\mathrm{\Omega }(X=0,W),$$
(53)
$`h^2`$ and $`\mathrm{\Omega }`$ can be expressed in terms of functions $`k(X,W)`$ and $`\mathrm{\Omega }_1(X,W)`$ that are regular analytic at $`W=0`$ as

$`h^2=h_0(W)^2+k(X,W)W^2,`$ (54a)
$`\mathrm{\Omega }=\mathrm{\Omega }_0(W)+\mathrm{\Omega }_1(X,W)W^2.`$ (54b)
Hence, in the advanced/retarded coordinates

$`u_\pm :=t\pm {\displaystyle ^W}h_0(W){\displaystyle \frac{dW}{W^2}},`$ (55a)
$`\varphi _\pm :=\varphi \pm {\displaystyle ^W}\mathrm{\Omega }_0(W)h_0(W){\displaystyle \frac{dW}{W^2}},`$ (55b)
the metric takes the form that is regular analytic at $`W=0`$:

$$ds^2=aW^2du_\pm ^2\pm 2h_0du_\pm dW+akdW^2+bX^2(d\varphi _\pm \pm \mathrm{\Omega }_1h_0dW\mathrm{\Omega }du_\pm )^2+cdX^2.$$
(56)
Since $`W=0`$ is a null hypersurface and the Killing vector $`\xi =_t`$ is its null tangent, the hypersurface $`W=0`$ is a Killing horizon. Because $`W=0`$ is the double root of $`\rho ^2/R^2`$, the horizon is degenerate, i.e., the surface gravity vanishes. Further, since $`\mathrm{\Omega }`$ vanishes at $`W=0`$, the horizons are non-rotating.

### III.3 Area and shape of the horizons

It is easy to estimate the area of the horizons. From (52) and (56), we obtain

$$A=2\times \underset{W0}{lim}2\pi _0^{\mathrm{}}\sqrt{bc}X𝑑X=\frac{16\pi \sigma ^2(p+1)}{p^2},$$
(57)
which coincides with (31) for the non-rotating limit $`p=1`$. Since it can be written as

$$A=4\pi M(M+\sqrt{M^2a^2})$$
(58)
with $`a=J/M`$, it is half of the area of the Kerr black hole with the same mass $`M`$.

Although the horizon area was finite in spite of its non-compactness in the case of ZVW2, the horizons become compact in the case of TS2. In order to see this, let us rewrite the metric of the two-dimensional section $`W=0,u_\pm =`$const of the horizon as

$$ds_H^2=\frac{2\sigma ^2(p+1)}{p^2\mathrm{}(\theta )}\left[d\theta ^2+\mathrm{}^2(\theta )\mathrm{sin}^2\theta d\varphi _\pm ^2\right],$$
(59)
where

$`X=\mathrm{tan}{\displaystyle \frac{\theta }{2}}(0\theta \pi ),`$ (60)
$`\mathrm{}(\theta )={\displaystyle \frac{4p^2}{(\mathrm{cos}\theta +p^2q^2)^2+4p^2q^2}}.`$ (61)
For $`0<p,q<1`$ with $`p^2+q^2=1`$, $`\mathrm{}(\theta )`$ is positive and finite for $`0\theta \pi `$. Hence, each horizon at P± is compact and homeomorphic to $`S^2`$. Further, since $`\mathrm{}(0)=1`$, it is smooth for $`0\theta <\pi `$. However, since $`\mathrm{}(\pi )=p^2/q^2`$, it has a conic singularity at $`\theta =\pi `$, i.e., at the point where it contacts with the segment I, except for the case $`p=q=1/\sqrt{2}`$, as illustrated in Fig. 10. The deficit angle is the same as that at the open segment I discussed in III.1.3; it is positive for $`q^2>p^2`$ and negative for $`q^2<p^2`$. In particular, in the limit $`q0`$, the deficit angle tends negative infinity, which is consistent with the fact that each regular portion of the horizon is a half sphere for ZVW2.

## IV Summary and Discussions

We have investigated the structure of the domains of the outer communication and their boundaries for the Zipoy-Voorhees-Weyl solution with arbitrary $`\delta `$ and for the Tomimatsu-Sato solution with $`\delta =2`$ and $`0<p,q<1`$(TS2). We have shown that three types of singularities arise for these spacetimes. One is the naked curvature singularity. It appears on the segment I: $`\rho =0,\sigma <z<\sigma `$ for the ZVW spacetime with $`\delta 0,1`$. We have shown that this singularity is point-like for $`\delta <0`$, string-like for $`0<\delta <1`$ and ring-like for $`1<\delta `$. The last feature has not been noticed so far, and has a close relation with the existence of the ring singularity for TS2. Another is an inextendible quasi-regular singularity at the segment I for TS2. We have confirmed that this singularity is conical locally, but has a more complicated helical structure globally, as was pointed out by Gibbons and Russell-ClarkGibbons.G&Russell-Clark1973 . The third is the directional quasi-regular singularities at the two points P±: $`(\rho ,z)=(0,\pm \sigma )`$ for the ZVW spacetime with $`\delta 2`$ and for TS2. We have proved that these points are degenerate Killing horizons for the ZVW spacetime with $`\delta =2,3`$ or $`\delta 4`$ and for TS2, by constructing regular extensions across these points explicitly.

An interesting feature of the horizon for the ZVW spacetime is that the two regular horizons at P± share the ring singularity at I as the common boundary and these pieces altogether form a single spherical horizon with a ring singularity at the equator. In contrast, although the two horizons of TS2 also contact with the conical segment I, each of them is a compact surface which is homeomorphic to $`S^2`$ and has a conical singularity at the point where it contacts I for $`pq`$. Thus, it is natural to regard that the horizon of TS2 has two disconnected components. This result is rather surprising because TS2 can be obtained from the Kramer-Neugebauer solution representing a superposition of two Kerr solutionsKramer.D&Neugebauer1980 , as the limit when the centers of two black holes coincide Oohara.K&Sato1981 . This may indicates a new possibility for the final states of gravitational collapse, if the cosmic censorship does not hold.

In this connection, it should be noted that the horizon of the ZVW solution with $`\delta =2`$(ZVW2) practically has a single component, as pointed above, although this solution can be also obtained as a similar limit from the Israel-Khan solutionIsrael.W&Kahn1964  with two black holes with the same mass. Further, we have pointed out that the naked curvature singularity of the ZVW spacetime can be regarded to have a positive gravitational mass. This should be contrasted with the fact that the ring singularity of TS2 has zero mass.

## Acknowledgement

H.K. would like to thank Misao Sasaki for valuable comments. This work was partly supported by the Grand-In-Aid for scientific researches by JSPS (Number. 11640273).

## Appendix A Geodesic equation in the ZVW spacetime

Geodesics in the ZVW spacetime are determined as solutions to the following set of equations:

$`\dot{t}=E\left({\displaystyle \frac{x+1}{x1}}\right)^\delta ,`$ (62a)
$`\dot{\varphi }=L{\displaystyle \frac{(x1)^{\delta 1}}{(x+1)^{\delta +1}(1y^2)}},`$ (62b)
$`{\displaystyle \frac{\dot{x}^2}{x^21}}+{\displaystyle \frac{\dot{y}^2}{1y^2}}={\displaystyle \frac{(x^2y^2)^{\delta ^21}}{(x^21)^{\delta ^2}}}F,`$ (62c)
$`F{\displaystyle \frac{d^2y}{dx^2}}=F{\displaystyle \frac{(\delta ^21)y(1y^2)}{(x^21)(x^2y^2)}}L^2{\displaystyle \frac{y(x1)^{2\delta 2}}{(1y^2)(x+1)^{2\delta +2}}}`$
$`+[ϵ{\displaystyle \frac{(x1)^{\delta 1}}{(x+1)^{\delta +1}}}F{\displaystyle \frac{2(x\delta )(x^2y^2)+(\delta ^21)x(1y^2)}{(x^21)(x^2y^2)}}`$
$`L^2{\displaystyle \frac{(x2\delta )(x1)^{2\delta 2}}{(1y^2)(x+1)^{2\delta +2}}}]{\displaystyle \frac{dy}{dx}}`$
$`\left[F{\displaystyle \frac{y\{x^2y^2+(1\delta ^2)(1y^2)\}}{(1y^2)(x^2y^2)}}+L^2{\displaystyle \frac{y(x1)^{2\delta 1}}{(1y^2)^2(x+1)^{2\delta +1}}}\right]\left({\displaystyle \frac{dy}{dx}}\right)^2`$
$`+[{\displaystyle \frac{\delta ϵ}{1y^2}}\left({\displaystyle \frac{x1}{x+1}}\right)^\delta F{\displaystyle \frac{x(x^21)+\delta ^2x(1y^2)2\delta (x^2y^2)}{(1y^2)(x^2y^2)}}`$
$`L^2{\displaystyle \frac{(x2\delta )(x1)^{2\delta 1}}{(1y^2)^2(x+1)^{2\delta +1}}}\left]\right({\displaystyle \frac{dy}{dx}})^3,`$ (62d)
where the dot denotes the derivative with respect to an affine parameter, and $`F`$ is given by

$$F:=E^2L^2\frac{(x1)^{2\delta 1}}{(1y^2)(x+1)^{2\delta +1}}ϵ\left(\frac{x1}{x+1}\right)^\delta .$$
(63)
$`E`$ and $`L`$ are the specific energy and the specific angular momentum of a particle, and $`ϵ`$ is equal to $`1`$ for time-like geodesics and $`0`$ for null geodesics.


text:
# Revisiting the Schrödinger Probability Current

## I Introduction

In this letter we revisit the definition of the probability current for the Schrödinger equation. We discuss the probability current of the hydrogen atom and the isotrop harmonic oscillator for which the classical results can be expressed in the form: ”the probability current of the wave eigenfunctions $`\psi _{nlm}`$ of Schrödinger with $`m=0`$ is nil”, cf. for example CohenTannoudji . Indeed, as the function $`\psi _{nl0}=R_{nl}(r)Y_{l0}(\theta ,\phi )`$ is real, the classical definition of the Schrödinger probability current leads to

$$𝐉_\mathrm{𝟏}=\frac{i\mathrm{}}{2m}(\psi \psi ^{}\psi ^{}\psi )=0.$$
(1)
First, we prove that for the Dirac equation, the probability current of stationary wave functions of the hydrogen atom and of the isotrop harmonic oscillator are not nil and correspond to a circular rotation of the probability.

Then, we recall how a good approximation of the Dirac probability current is obtained in the nonrelativist case. It is necessary to add in equation (1) the Gordon current Gordon28-1  which reads

$$𝐉_\mathrm{𝟐}=\frac{\mathrm{}}{2m}rot(\psi ^{}\sigma \psi )$$

in the Pauli approximation Gurtler Landau89-1  and

$$𝐉_\mathrm{𝟐}=\frac{1}{m}\rho \times 𝐬$$
(2)
in the Schrödinger approximation Holland99 Colijn , where $`𝐬`$ corresponds at a constant spin vector. Consequently, we get a circular probability current in the Schrödinger approximation for the hydrogen atom and the isotrop harmonic oscillator.

## II The Dirac probability current of the hydrogen atom

In the presence of electromagnetic couplings, the Dirac equation reads

$$i\gamma ^\mu \frac{\psi }{x_\mu }+\frac{e}{\mathrm{}}\gamma ^\mu A_\mu \psi +\frac{mc}{\mathrm{}}\psi =0$$
(3)
where the $`\gamma ^\mu `$ ($`\mu =0,1,2,3`$) are the Dirac matrices.

The probability density $`\rho =\psi ^{}\psi `$ and the three components of the probability current $`𝐉`$:

$$J^i=c\psi ^{}\gamma ^0\gamma ^i\psi $$
(4)
satisfy the continuity equation

$$\frac{\rho }{t}+.𝐉=0.$$
(5)
In the case of an electron in the hydrogen atom with a spherically symmetric potential $`V=\frac{e}{4\pi r}`$ and a potential vector $`𝐀=0`$, the Dirac eigenvalues equation is then

$$\gamma _0(\frac{E}{\mathrm{}c}+\frac{eV}{\mathrm{}c})\psi (𝐫)i𝚪.\psi (𝐫)+\frac{mc}{\mathrm{}}\psi (𝐫)=\mathrm{𝟎}.$$

The eigenfunctions $`\psi _{n,l,j,m}(𝐫)`$ depend on the quantum numbers $`n,l,j`$ and $`m`$ with $`n`$ integer ($`n1`$), $`l`$ integer ($`0ln1`$), $`j=\frac{1}{2}`$ if $`l=0`$ and $`j=l+m`$ if $`l1`$, $`m=\pm \frac{1}{2}`$. Then $`\psi _{n,l,j,m}(𝐫)`$ is given Landau89-1  by

$$\psi _{n,l,j,m}(𝐫)=\left(\begin{array}{c}f(r)\mathrm{\Omega }_{jlm}\\ ig(r)\mathrm{\Omega }_{jl^{}m}\end{array}\right)$$

with $`l^{}=j+\frac{1}{2}`$ if $`l=j\frac{1}{2}`$ and $`l^{}=j\frac{1}{2}`$ if $`l=j+\frac{1}{2}`$, and where $`f(r)`$ and $`g(r)`$ are real functions of $`r`$ and where

$$\mathrm{\Omega }_{l\pm \frac{1}{2},l,m}=\left(\begin{array}{c}\pm \sqrt{\frac{l+\frac{1}{2}\pm m}{2l+1}}Y_{l,m\frac{1}{2}}\\ \sqrt{\frac{l+\frac{1}{2}m}{2l+1}}Y_{l,m+\frac{1}{2}}\end{array}\right).$$

We obtain four classes of wave eigenfunctions: $`m=\frac{1}{2}`$ or $`\frac{1}{2}`$, $`l=0`$ or $`l1`$. For $`m=\frac{1}{2}`$

$$\psi _{n,l,j,m=+\frac{1}{2}}(r,\theta ,\phi )=\left(\begin{array}{c}f(r)\left(\begin{array}{c}a(\theta )\\ b(\theta )\mathrm{exp}(i\phi )\end{array}\right)\\ ig(r)\left(\begin{array}{c}c(\theta )\\ d(\theta )\mathrm{exp}(i\phi )\end{array}\right)\end{array}\right)$$
(6)
with for $`l=0`$: $`a(\theta )=Y_{00}`$, $`b(\theta )=0`$, $`c(\theta )=Y_{00}\mathrm{cos}\theta `$, $`d(\theta )=Y_{00}\mathrm{sin}\theta `$ and, for $`l>0`$: $`a(\theta )=\sqrt{\frac{l+1}{2l+1}}Y_{l,0}`$, $`b(\theta )=\sqrt{\frac{l}{2l+1}}Y_{l,1}\mathrm{exp}(i\phi )`$, $`c(\theta )=\sqrt{\frac{l+1}{2l+3}}Y_{l+1,0}`$ and $`d(\theta )=\sqrt{\frac{l+2}{2l+3}}Y_{l+1,1}\mathrm{exp}(i\phi )`$, and for $`m=\frac{1}{2}`$

$$\psi _{n,l,j,m=\frac{1}{2}}(r,\theta ,\phi )=\left(\begin{array}{c}f(r)\left(\begin{array}{c}b(\theta )\mathrm{exp}(i\phi )\\ a(\theta )\end{array}\right)\\ ig(r)\left(\begin{array}{c}d(\theta )\mathrm{exp}(i\phi )\\ c(\theta )\end{array}\right)\end{array}\right)$$
(7)
with for $`l=0`$: $`a(\theta )=0`$, $`b(\theta )=Y_{00}`$, $`c(\theta )=Y_{00}\mathrm{sin}\theta `$, $`d(\theta )=Y_{00}\mathrm{cos}\theta `$ and, for $`l>0`$: $`a(\theta )=\sqrt{\frac{l+1}{2l+1}}Y_{l,1}\mathrm{exp}(i\phi )`$, $`b(\theta )=\sqrt{\frac{l}{2l+1}}Y_{l,0}`$, $`c(\theta )=\sqrt{\frac{l1}{2l1}}Y_{l1,1}\mathrm{exp}(i\phi )`$ and $`d(\theta )=\sqrt{\frac{l}{2l1}}Y_{l1,0}`$.

For a wave function $`\psi =\left(\begin{array}{c}\chi _1\\ \chi _2\end{array}\right)`$ and the Dirac matrices

$$\gamma ^0=\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right)and\gamma ^i=\left(\begin{array}{cc}0& \sigma ^i\\ \sigma ^i& 0\end{array}\right)$$

the probability density is equal to: $`\rho =\chi _1^{}\chi _1+\chi _2^{}\chi _2`$, and the probability current is given by: $`J^i=c(\chi _1^{}\sigma ^i\chi _2+\chi _2^{}\sigma ^i\chi _1)`$.

For the functions (6) and (7), we find : $`\rho =f^2(r)[a^2(\theta )+b^2(\theta )]+g^2(r)[c^2(\theta )+d^2(\theta )]`$, $`J^1=2cf(r)g(r)[adbc]\mathrm{sin}\phi `$, $`J^2=2cf(r)g(r)[adbc]\mathrm{cos}\phi `$, $`J^3=0`$, and as probability current velocity:

$$𝐯=\frac{𝐉}{\rho }=\frac{2cf(r)g(r)[adbc]}{f^2(r)[a^2+b^2]+g^2(r)[c^2+d^2]}𝐮_\phi $$
(8)
which corresponds to circular probability currents.

We have the same result (8) for all the Dirac equations with a spherically symetric potential; in particular for the isotrop harmonic oscillator.

The Dirac eigenfunctions are not stationary. All the electron wave eigenfunctions of the hydrogen atom correspond to circular probability currents. It is not the case for the Schrödinger equation with the current given by (1), since we have then stationary solutions for the eigenfunctions $`\psi _{nl0}`$. In the fondamental state $`1s_{\frac{1}{2}}`$ with the energy $`E_{n,\frac{1}{2}}=m_ec^2\sqrt{1\alpha ^2}`$, we have: $`f(r)=a_0(\frac{r}{r_0})^{\sqrt{1\alpha ^2}1}e^{\frac{r}{r_0}}`$, $`g(r)=a_0\frac{1\sqrt{1\alpha ^2}}{\alpha }(\frac{r}{r_0})^{\sqrt{1\alpha ^2}1}e^{\frac{r}{r_0}}`$, where $`r_0=\frac{\mathrm{}^2}{m_ee^2}=0.53\AA `$ is the Bohr radius. From (8), we deduce $`\rho (r,\theta ,\phi )=\rho (r)=2a_0^2\frac{1\sqrt{1\alpha ^2}}{\alpha ^2}(\frac{r}{r_0})^{2(\sqrt{1\alpha ^2}1)}e^{2\frac{r}{r_0}}`$, and

$$𝐯=\alpha c\mathrm{sin}\theta 𝐮_\phi .$$
(9)
From (7) and (8), we deduce $`𝐯=\alpha c\mathrm{sin}\theta 𝐮_\phi `$. There is a big difference between (9) and the current nil given by (1) with the Schrödinger equation. The lesson of this example is that, although Schrödinger equation is a good approximation of the Dirac equation, the Schrödinger probability current (1) is not the good approximation of the Dirac probability current (4). This conclusion appears clearly, as we recall now Gordon28-1 Gurtler Landau89-1 , in the approximations to transform first the bispinor of the Dirac equation to the spinor of the Pauli equation, then the Pauli spinor to the wave function of the Schrödinger equation.

## III The Pauli and Schrödinger probability currents

After the change of variable $`\psi =\psi ^{}e^{\frac{i}{\mathrm{}}mc^2t}`$, the Dirac equation (3) for the bispinor $`\psi =(\chi _1,\chi _2)`$ can be written

$$[i\mathrm{}\frac{}{t}eV]\chi _1=c\sigma (i\mathrm{}\frac{e}{c}𝐀)\chi _2$$
(10)
$$[i\mathrm{}\frac{}{t}eV+2mc^2]\chi _2=c\sigma (i\mathrm{}\frac{e}{c}𝐀)\chi _1.$$
(11)
Small electron velocities allow us to neglect the two first terms of the left part of equation (11) and to eliminate $`\chi _2`$ in (10). So we obtain the Pauli equation for the spinor $`\chi _1`$

$$i\mathrm{}\frac{\chi _1}{t}=[\frac{1}{2m}(i\mathrm{}\frac{e}{c}𝐀)^2+eV\frac{e\mathrm{}}{2mc}\sigma .𝐇]\chi _1$$

where $`𝐇=\times 𝐀`$ is the magnetic field.

As $`\chi _2\chi _1`$, the Dirac density $`\rho =\psi ^{}\psi =\chi _1^{}\chi _1+\chi _2^{}\chi _2`$ is equal in first approximation to the density $`\rho =\chi _1^{}\chi _1`$ of the Pauli equation. Then the Dirac current density can be written

$$𝐉=i\frac{\mathrm{}}{2m}(\chi _1\chi _1^{}\chi _1^{}\chi _1)\frac{e}{mc}𝐀\chi _1\chi _1^{}+\frac{\mathrm{}}{2m}\times (\chi _1^{}\sigma \chi _1)$$
(12)
Thus, to obtain in Pauli a good approximation of the Dirac probability current, it is necessary to add the Gordan current Gordon28-1

$$𝐉_\mathrm{𝟐}=\frac{\mathrm{}}{2m}\times (\chi _1^{}\sigma \chi _1)$$

to the classical Pauli current

$$𝐉_\mathrm{𝟏}=i\frac{\mathrm{}}{2m}(\chi _1\chi _1^{}\chi _1^{}\chi _1)\frac{e}{mc}𝐀\chi _\mathrm{𝟏}\chi _\mathrm{𝟏}^{}.$$

Then, the current $`𝐉=𝐉_\mathrm{𝟏}+𝐉_\mathrm{𝟐}`$ verifies, as the current $`𝐉_\mathrm{𝟏}`$, the continuity equation(5).

To obtain the Schrödinger equation from Pauli equation, consider the case where we have no magnetic field Holland99  and where the system is in a spin eigenstate

$$\chi _1(𝐫,t)=\phi (𝐫,t)\chi $$

with a constant spinor $`\chi `$ such as $`\chi ^{}\chi =1`$. Then the function $`\phi (𝐫,t)`$ verifies the Schrödinger equation

$$i\mathrm{}\frac{d\phi }{dt}=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\phi +eV\phi .$$
(13)
However, in this case, the current (12) does not reduce Holland99  to classical Schrödinger probability current. Indeed, writing $`\phi (𝐫,t)=\sqrt{\rho }e^{i\frac{S}{\mathrm{}}}`$, we have

$$𝐉=\frac{\rho }{m}S+\frac{1}{m}\rho \times 𝐬$$

with

$$𝐬=\frac{\mathrm{}}{2}(\chi ^{}\sigma \chi ).$$

Indeed, to obtain a good approximation of the Dirac probability current for the Schrödinger equation, it is necessary to add the spin-dependent current

$$𝐉_\mathrm{𝟐}=\frac{1}{m}\rho \times 𝐬$$

to the classical Schrödinger current

$$𝐉_\mathrm{𝟏}=\frac{\rho }{m}S.$$

Therefore the probability current depends on a constant spin vector $`𝐬`$. The add of a constant spin vector to the Schrödinger equation is not new, and has been used before by Landau Landau89-2 .

This current $`𝐉_\mathrm{𝟐}`$ gives also a simple explanation of the value 2 of the gyromagnetic ratio of the Dirac electron; Indeed, this contribution to the particle orbital angular momentum is

$$\text{L}_\text{2}=\text{r}\times m\text{J}_\text{2}=\text{r}\times \rho \times \text{s}$$

and this mean angular momentum is then

$$<\text{L}_\text{2}>=\text{L}_\text{2}d^3r=(\text{r}\text{s})\rho d^3r(\text{r}\rho )\text{s}d^3r=\text{s}+3\text{s}=2\text{s}.$$

We can verify the precedent theoritical approach with the computation of the current (13) of the hydrogen atom wave eigenfunctions in the approximation of Schrödinger.

For the wave functions $`\psi _{nlm}(r,\theta ,\phi )=R_{nl}(r)P_{lm}(\theta )\mathrm{exp}(im\phi )\mathrm{exp}(i\frac{E_nt}{\mathrm{}})`$, the classical Schrödinger current is

$$𝐉_\mathrm{𝟏}=\frac{\rho }{m_e}\frac{m\mathrm{}}{r\mathrm{sin}\theta }𝐮_\phi $$

with $`\rho =R_{nl}^2(r)P_{lm}^2(\theta )`$. The current is nil for the wave functions with $`m=0`$, in particular for the states $`ns`$. If $`m0`$, we take as spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$ and $`𝐬=\frac{\mathrm{}}{2}𝐤`$ if $`m<0`$. Then, this vector can be written

$$𝐬=sgn(m)\frac{}{}(\mathrm{cos}\theta 𝐮_𝐫\mathrm{sin}\theta 𝐮_\phi )$$

and the additional spin-dependent current (Gordon current) $`𝐉_\mathrm{𝟐}`$ is then equal to

$$𝐉_\mathrm{𝟐}=sgn(m)\frac{\mathrm{}}{m_e}(R_{nl}R_{nl}^{}P_{lm}^2\mathrm{sin}\theta +\frac{R_{nl}^2}{r}P_{lm}P_{lm}^{}\mathrm{cos}\theta )𝐮_\phi .$$

Then, the probability current velocity is

$$𝐯=\frac{𝐉}{\rho }=\frac{\mathrm{}}{m_e}[sgn(m)(\frac{R_{nl}^{}(r)}{R_{nl}(r)}\mathrm{sin}\theta +\frac{P_{lm}^{}(\theta )}{rP_{lm}(\theta )}\mathrm{cos}\theta )+\frac{m}{r\mathrm{sin}\theta }]𝐮_\phi $$

which corresponds to circular probability currents exactly as for the Dirac equation (8). In the case of the fondamental state $`1s`$, we find again exactly the velocity (9) of the Dirac probability current

$$𝐯_{1s}=\alpha c\mathrm{sin}\theta 𝐮_\phi .$$

For the states $`2s`$, $`2p_0`$, $`2p_1`$ and $`2p_1`$, we obtain for the current velocity

$$𝐯_{2s}=\frac{\alpha c}{2}(1+\frac{1}{1\frac{r}{2r_0}})\mathrm{sin}\theta 𝐮_\phi ,$$

$$𝐯_{2p0}=𝐯_{2p_1}=\frac{\alpha c}{2}\mathrm{sin}\theta 𝐮_\phi and𝐯_{2p_1}=\frac{\alpha c}{2}\mathrm{sin}\theta 𝐮_\phi .$$

## IV Conclusion

As we have shown, although the Schrödinger equation is a good approximation of the Dirac equation, the Schrödinger probability current (1) is not the good approximation of the Dirac probability current (4). The main conclusion is that it is necessary to add to the classical Schrödinger current the spin-dependent current

$$𝐉_\mathrm{𝟐}=\frac{1}{m}\rho \times 𝐬$$

which corresponds to a constant spin vector $`𝐬`$.

In the case of the isotrop harmonic oscillator with $`V(𝐫)=\frac{1}{2}m\omega ^2𝐫^2`$, the fondamental wave function $`1s`$

$$\psi _{0,0,0}(𝐫)=(\frac{m\omega }{\pi \mathrm{}})^{\frac{3}{4}}\mathrm{exp}(\frac{m\omega }{2\mathrm{}}𝐫^2)\mathrm{exp}(\frac{3i\omega t}{2})$$

gives a Schrödinger classical current $`𝐉_\mathrm{𝟏}`$ nil. If we consider an oscillator with a constant spin $`𝐬=\frac{\mathrm{}}{2}𝐤`$, then the equation (13) gives a circular probability current $`𝐉`$ with the velocity

$$𝐯_{1s}=\frac{𝐉}{\rho }=\omega r\mathrm{sin}\theta 𝐮_\phi .$$
(14)
It is the same velocity as in the classical case. For the other eigenfunctions of the harmonic oscillator, we find also a non nil probability current, but non circular as in (14).