Hyperbolic Relaxation Method for Elliptic Equations
Hannes R. Rüter David Hilditch Marcus Bugner Bernd Brügmann



 




^1Theoretical Physics Institute, University of Jena, 07743 Jena, Germany.

^2CENTRA, University of Lisbon, 1049 Lisboa, Portugal.




 
We show how the basic idea of parabolic Jacobi relaxation can be
modified to obtain a new class of hyperbolic relaxation schemes that
are suitable for the solution of elliptic equations. Some of the
analytic and numerical properties of hyperbolic relaxation are
examined. We describe its implementation as a first order system in a
pseudospectral evolution code, demonstrating that certain elliptic
equations can be solved within a framework for hyperbolic evolution
systems. Applications include various initial data problems in
numerical general relativity. In particular we generate initial data
for the evolution of a massless scalar field, a single neutron star,
and binary neutron star systems.




  02.60.Cb,    
  02.70.Hm,    
  04.25.D-     

































Hyperbolic Relaxation Method for Elliptic Equations
    Hannes R. Rüter^1, David Hilditch^1,2, Marcus
Bugner^1, and Bernd Brügmann^1
    December 30, 2023
================================================================================








§ INTRODUCTION



The solution of elliptic partial differential equations (elliptic
PDEs) is an important problem in many areas of physics, including
fluid dynamics, quantum mechanics, general relativity, and many more
besides. Correspondingly large is the variety of analytic and
numerical methods dealing with the solution of elliptic PDEs. The
starting point for many methods is a discretization and (if required)
a linearization, which for typical problems arising in physics leads
to a sparse system of linear equations for a large but finite number
of degrees of freedom.

A key role in the solution of linear systems is played by iterative
methods, <cit.>. Among the basic iterative methods
are relaxation methods, in particular Gauss-Seidel and Jacobi
relaxation methods, and the family of Krylov subspace methods. Closely
connected are strategies to accelerate the convergence of such
methods, such as preconditioners and multigrid methods.

In this paper we study a modification of the classic Jacobi method,
which is motivated by its origin in physical relaxation problems. For
concreteness, consider as a minimal example the Laplace equation,

    Δϕ = 0,

for a function ϕ(x,y,z) on a regular subset of
ℝ^3 together with appropriate boundary conditions. The
Jacobi relaxation scheme can be obtained by introducing a time
parameter t and considering instead of (<ref>) the parabolic
diffusion equation

    ∂_t ϕ = Δϕ.

As time approaches infinity, any initial data for ϕ “relaxes”
to a stationary state, where ∂_t ϕ = 0 and hence
Δϕ=0 (<ref>) is satisfied as well. The Jacobi
iteration method is obtained by discretizing the diffusion equation
(<ref>). In essence, we introduce an “unphysical” time
dependence, which is not part of the original problem, and obtain the
solution to the time-independent problem by means of a fixed point
iteration.

In this paper we investigate a similar strategy, which however relies
on a different type of evolution equation. Instead of replacing the
elliptic equation (<ref>) by the parabolic equation
(<ref>), we consider a hyperbolic wave equation with
damping,

    ∂_t^2 ϕ + ∂_t ϕ = Δϕ.

Combining time derivatives as in (<ref>) adds
strong diffusion to the pure wave equation while maintaining the
hyperbolic character of the PDE. The idea is that deviations from the
stationary state satisfying Δϕ=0 are damped to zero or are
propagated away, and furthermore it can be advantageous to perform
hyperbolic as opposed to parabolic evolutions.

In the limit of vanishing damping we obtain the wave equation

    ∂_t^2 ϕ = Δϕ.

If a stationary state is reached, we again have
solved (<ref>). Experimenting with (<ref>) we
found that the damping is the main desirable feature, while
propagating waves off the grid is far less relevant for the reduction
of the residual.  Therefore, the proposal is to study
(<ref>) for a variety of elliptic operators with strong
damping.

Hyperbolic equations containing diffusion or viscosity terms are a
well-known topic, <cit.>, which is relevant
to our discussion, see Sec. <ref>. The
topic revolves around applying relaxation methods to hyperbolic
equations, in particular hyperbolic conservation laws with diffusion.
Although the type of equations that are considered are similar, our
perspective is different.  The elliptic equation is for us the
fundamental problem, and we add a hyperbolic, damped time-dependence
to obtain an iterative scheme for the solution of the elliptic
equation.  Since there does not seem to be an established name for
this idea, we refer to the method as hyperbolic relaxation for
  elliptic equations (HypRelax), as opposed to parabolic relaxation
that is at the heart of the Jacobi method.

With regard to previous literature on hyperbolic relaxation for
elliptic equations, some aspects have been explored
in <cit.> in the context of “gauge drivers” for
numerical relativity.  In particular, <cit.> introduced
one of the most used gauge conditions for certain black hole
evolutions, the Gamma-driver for the shift vector, which employs a
hyperbolic equation related to the elliptic equation for a minimal
distortion shift. Also see <cit.> on gauge drivers, where
however only parabolic relaxation is considered.

The goal of the present paper is to develop hyperbolic relaxation
given by the prototype in (<ref>) into a method to solve
a general class of second order, non-linear elliptic equations. The
basic observation that elliptic equations can be related to the
stationary end points of evolutions is well known. However, the
challenge is to develop a concrete formulation of a hyperbolic
relaxation method that can solve non-trivial problems.[ We
  are not aware of other work on hyperbolic relaxation introduced
  explicitly for the solution of non-trivial elliptic equations, but
  given the simplicity of the idea, it would be surprising if there is
  none. The authors welcome any pointers to existing literature.]  The
problem of immediate interest to us is defined by the constraint
equations of general relativity, which we solve as a system of
non-linear elliptic equations to obtain initial data for evolution in
numerical relativity. However, the formalism is quite independent of
this particular problem.

The main result of this paper is that hyperbolic relaxation can be
formulated for non-trivial equations, and numerical experiments for
some specific test problems were indeed successful. Specific examples
include the Poisson equation, the conformal-thin-sandwich equations,
scalar field initial data, and some simple configurations of single
and binary neutron stars.

Considering (<ref>), let us collect some basic
observations here in order to introduce the main questions we want to
address. First of all, we have to address the well-posedness of the
hyperbolic PDEs. Given a self-adjoint, elliptic operator, the
hyperbolicity of equations of type (<ref>) should be
clear. We demonstrate this below for a general class of
equations. There exists a rich theoretical background regarding
well-posedness and numerical stability for hyperbolic
PDEs <cit.>, which helps to find relaxation
schemes that are well suited for numerical applications. However,
evolutions of hyperbolic PDEs are not trivial, so we should expect
that some elliptic problems are not amenable to hyperbolic relaxation
while others are, which is why we include a non-trivial set of test
problems below.

Second, in addition to the boundary conditions of the original
elliptic equation we have to choose boundary conditions for the
hyperbolic equations that are compatible with the asymptotic elliptic
problem This choice is not unique, but of great importance to obtain
successful evolutions. In particular, we consider maximally
dissipative boundary conditions.

Third, assuming feasibility and stability of hyperbolic relaxation, a
key question concerns the efficiency of the method. In both parabolic
and hyperbolic relaxation methods the time parameter is unrelated to
the elliptic equation, the time evolution is of no interest as
long as the stationary state is reached efficiently. This is the basis
for different acceleration strategies. For hyperbolic relaxation,
there is a finite propagation speed, and in contrast to the diffusion
equation it is not clear how to by-pass that speed to accelerate the
method.

Beyond the intrinsic interest in a new method, we have to ask whether
hyperbolic relaxation, after some significant further development that
is beyond the scope of this paper, might become an interesting
alternative to the highly developed standard methods.

As it stands, there are pragmatic considerations that can make
hyperbolic relaxation methods interesting, in particular when solving
elliptic equations as part of a larger project. For example, elliptic
PDEs are often solved to provide initial data for evolution systems
that are subject to certain constraint equations, the Maxwell
equations or the Einstein equations. However, the main work load is
the actual evolution of the data by integrating a hyperbolic PDE.  In
such a case the hyperbolic relaxation method does not have to compete
with optimized standard methods in terms of efficiency as long as
solving the elliptic equation is only a small part of the entire work
load.  On the other hand, a hyperbolic relaxation method may be easy
to implement using the existing infrastructure of a numerical
evolution code, avoiding the need for and the complications of an
external elliptic solver.  Using the same infrastructure also has the
advantage that interpolation errors can be avoided by using the same
grid discretization. Considering our research in numerical relativity,
a sophisticated infrastructure for evolutions is indeed available, but
we were looking for alternative elliptic solvers. Hence we implemented
hyperbolic relaxation in the pseudospectral hyperbolic evolution code
 <cit.>, which only required minor
modifications once the formalism itself was established.

The paper is structured as follows. In Sec. <ref>
we derive and motivate the evolution equations on which our relaxation
method is based and discuss some of its properties. In
Sec. <ref> we give details on the numerical
implementation and state the methods used in the  code.
In Secs. <ref> and <ref> we
apply our hyperbolic relaxation method to some test cases and to the
construction of initial data for numerical relativity. We conclude in
Sec. <ref>.

Throughout the paper we use the Einstein summation convention, we
sum over indices that occur once as an upper index and once as a lower
index, s^i t_i = ∑_i s^i t_i.  Latin letters i, j, k, ...
denote coordinate components and they are lowered and raised by an
arbitrary metric with positive signature.  An index s denotes a
contraction with a vector s_i, in particular ∂_s = s^i
∂_i = s^i ∂/∂ x^i.  Greek letter indices
α, β, γ denote components of a field and they are
lowered and raised by the Euclidean metric.  We also use Latin letters
a, b, c to denote the spacetime components in general relativity and
which are lowered and raised by the spacetime metric g_ab.




§ THE HYPERBOLIC RELAXATION EQUATIONS






 §.§ Evolution System



In the following we present the principal ideas of the hyperbolic
relaxation method and derive the equations that follow for the
iteration scheme. In this paper we consider only systems of elliptic
equations given in second order form, 
    (ψ)_α = a(x^k)^i^j^β_α∂_i
    ∂_j ψ_β + F_α(x^k,ψ_β, ∂_i
    ψ_β) = 0   ,

where the ψ_α are the N unknown solution variables and F
is a continuous function of the solution variables, their derivatives
and the D coordinates x^k. We take a^i^j^β_α
to be a smooth function of the coordinates and suppress this
dependence in the notation. In the following we consider
classically elliptic systems <cit.> only, systems
with

    (a(x^k)^i^j^β_α s_i s_j) ≠ 0     ∀ s
    ∈ℝ^D∖{0}  ,

where the determinant is understood to be taken on the indices
α and β.

Every elliptic system will be accompanied by a set of boundary
conditions on the variables ψ_α and we discuss their
treatment in section <ref>. We can
reduce the second order elliptic system to first order by introducing
the reduction variables r_iα:

    0    = a^i^j^β_α∂_i r_jβ +
    F_α(ψ_β, r_iβ)   ,
    
    0    = ∂_i ψ_α - r_iα  .

To solve the second order equation (<ref>) one
could employ the Jacobi method, which can be motivated by evolving the
parabolic partial differential equation:

    ∂_t ψ_α = (ψ)_α  ,

where t is some parameter that plays the role of time.

For a classically elliptic system with constant coefficients the
Jacobi method can only converge if a^j^i^β_α is
positive definite on the whole domain, there exists an
ϵ>0

    a^j^i^β_α t_jβ t_i^α≥ϵ t^iα t_iα    ∀ t ∈ℝ^D
      × N∖{0}  ,

which we assume in the rest of this paper. This condition corresponds
to the notion of strong ellipticity, which defines an important
subclass of classically elliptic systems <cit.>.  Note that we
have the freedom to multiply the elliptic equation with an invertible
matrix d^β_α, yielding (ψ)_α =
d^β_α (ψ)_β = 0, which has the same
solutions as the original equation. This freedom allows us to
transform some systems into a strongly elliptic system, that
originally were not.

In analogy to the Jacobi method (<ref>) we evolve
ψ_α by taking Eq. (<ref>) as the
right-hand side, yielding

    ∂_t ψ_α = a^i^j^β_α∂_i
    r_jβ + F_α(x^i,ψ_β, r_iβ)

and we proceed similarly with the equations for the reduction
variables r_i:

    ∂_t r_iα = b^j_i^β_α ( ∂_j
    ψ_β - r_jβ )   ,

where b^j_i^β_α is arbitrary under the
requirement of positive definiteness, meaning in analogy to
Eq. (<ref>)

    b^j_i^β_α t_jβ t^iα > ϵ
    t_iα t^iα    ∀ t ∈ℝ^D ×
      N∖{0}  .

The system of Eqs. (<ref>)
and (<ref>) forms a first order hyperbolic
differential equation which we refer to as the the hyperbolic
relaxation system. Clearly the reduction constraint
Eq. (<ref>) is not enforced at all times and
will indeed be violated during the relaxation process, however we
are only interested in the steady state, which fulfills the
reduction constraint, because Eq. (<ref>)
constantly drives the reduction variable r_i towards ∂_i
ψ. To see this, let us assume for arguments sake that
∂_t ∂_i ψ_α = 0. Then the solution for
r_iα has the form

    r_iα(t) = ∑_l=1^n e^- λ_l t∑_k_l=0^m_lx^k_l h_k_l i α + ∂_i ψ_α  ,

where h is constant and the λ_l are the n eigenvalues
defined by the eigenvalue equation b^j_i^β_α
t_l^i^α = λ_l t_l^j^β and m_l
depends on the geometric multiplicity of λ_l.  From the
positive definiteness of b we know that all the eigenvalues have
positive real part and it follows immediately that r_iα
approaches ∂_i ψ_α exponentially. We emphasize
however that in some cases, if the elliptic system has no
solution, ∂_i ψ_α can grow faster than r_i α
and thus the reduction constraints cannot be satisfied asymptotically
in time.

Assuming that the solution to the elliptic system exists and is unique
under the provided boundary conditions, then it is obvious that if the
hyperbolic relaxation system reaches a steady state, then we must have
the solution to the first order elliptic system
Eqs. (<ref>)
and (<ref>), and also of the original elliptic
equation (<ref>).

The relationship between hyperbolic equations and parabolic diffusion
equations has already been investigated in some special
cases <cit.>.  In particular it can be shown that
for large times t the solution of the hyperbolic equation

    ∂_t^2 ϕ + ∂_t ϕ - Δϕ = 0

will tend towards the solution of the parabolic PDE

    ∂_t ϕ - Δϕ = 0   .

The hyperbolic equation can be cast in first order form by introducing
the reduction variables ψ = ∂_t ϕ + ϕ and r_i =
∂_i ϕ, yielding the system

    ∂_t ϕ   = ψ - ϕ  , 
    ∂_t ψ   = δ^ij∂_i r_j   , 
    ∂_t r_i    = ∂_i ψ - r_i   .

Clearly the first of these equations is an ordinary differential
equation that has no direct counterpart in our hyperbolic relaxation
system, however it is directly evident that ϕ will tend towards
ψ exponentially. Thus it is plausible that for large t the
variable ψ of our hyperbolic relaxation system will behave
similarly to that of the Jacobi-type relaxation
equation (<ref>).




 §.§ Residual Evolution



The residuals of the first order system,
Eqs. (<ref>),(<ref>),
are given by

    R_α   = a^i^j^β_α∂_i r_jβ +
    F_α(x^i,ψ_β, r_iβ)   , 
    
    R_iα   = b^j_i^β_α ( ∂_j ψ_β -
    r_jβ )   .

A simple calculation shows that the residuals will evolve according to

    ∂_t R_α   = a^i^j^β_α∂_i
    R_jβ + ∂ F_α/∂ψ_β R_β +
    ∂ F_α/∂ r_i β R_i β  , 
    ∂_t R_iα   = b^j_i^β_α ( ∂_j R_β - R_jβ )   .

For a working relaxation scheme, we want the residual evolution system
to be stable, the first order residuals should converge to zero
for t →∞, for residuals that are sufficiently close
to zero.  Systems of this type and stability conditions are discussed
in detail in <cit.> and <cit.>. It is not
possible for us to give general results on the stability of the
hyperbolic relaxation scheme, as the multitude of possible systems is
too large to be covered in a closed form, especially for elliptic
systems with more than one variable.  A stability analysis must
therefore be done individually for the concrete problem.




 §.§ Mode Analysis



To shed some light on the behavior of solutions to the hyperbolic
relaxation equation (<ref>), we perform a simple
mode analysis, ignoring the issue of boundary conditions. Introduce
the plane-wave ansatz

    ϕ_pw(t,x) = e^i(k x - ω t),

where k and ω are constants. The wavenumber k is a real
number related to the wave length, k = 2π/λ, while ω
may be a complex number.  Inserted into the hyperbolic relaxation
equation (<ref>), we obtain

    ω^2 + i ω   = k^2,
    ω_±(k)    = - 1/2 (i ±√(4k^2 - 1)).

Recall that for the wave equation ω_±(k) = ± k, while for
the heat equation ω(k) = -i k^2. For hyperbolic relaxation,
there is a further case distinction for the sign under the square root
√(4k^2 - 1).

For sufficiently large wavenumber,

    ϕ_pw =
      e^-1/2t e^i(k x ±1/2√(4k^2 - 1) t),
         k ≥1/2,

which is a damped wave with phase velocity
v(k)=√(1-1/4k^2).  The damping is independent of k (as
opposed to the heat equation with e^-k^2 t). The phase velocity
approaches v=1 for large k, but for k approaching the critical
value 1/2 from above the phase velocity tends towards v=0.

For sufficiently small wavenumber,

    ϕ_pw = e^-1/2(1±√(1-4k^2))t e^i k x,    0
       ≤ k ≤1/2,

which is a non-moving wave profile e^i k x times a k-dependent
damping factor.  For k = 1/2, the damping is
e^-1/2t, while for k equal to zero there are two cases,
-1/2(1±√(1-4k^2))t = 0 or -t. For small k, the
worse (more weakly) damped case is -1/2(1±√(1-4k^2))t
≈ - k^2 t, which is the same damping as for the basic heat
equation.





Summarizing, the plane-wave mode analysis suggests that solutions to
the hyperbolic relaxation equation exhibit a mixture of relaxation and
wave propagation phenomena, see Fig. <ref>.  For wave
numbers larger than a critical value, k ≥ k_crit with k_crit
= 1/2, there is wave propagation with simultaneous damping.
Noteworthy is that the damping is independent of k,
e^-1/2t. This is a promising feature compared to parabolic
relaxation with e^-k^2 t for intermediate values of k.  For
large values of k parabolic relaxation has much stronger damping,
but the overall convergence rate is dominated by small k.  For
hyperbolic relaxation, there is no wave propagation for k ≤
k_crit, but the damping persists. Interestingly, the damping factor
asymptotes towards e^-k^2 t for k→0, and is never
worse than parabolic relaxation for small k.

The existence of a transition at a specific length scale
λ_crit=2π/k_crit = 4π signals that our choice
of hyperbolic relaxation equation has fixed a scale. Let us
generalize (<ref>) to

    ϵ∂_t^2 ϕ + η∂_t ϕ = Δϕ,

where ϵ and η are real, non-negative constants. A third
constant in front of Δϕ, the velocity squared v^2 in the
basic wave equation, has been rescaled to one without loss of
generality. With ϵ=η=1 as in (<ref>)
we fix the unit of time to be dimensionless (unity since [T^2]=[T]).
If we mimic (<ref>) and
(<ref>) by

    ∂_t ψ = ∂_j r^j,   ∂_t r_i = b (∂_i
    ψ - r_i),

with a real constant b, then we obtain (<ref>) with
ϵ = 1/b, η=1, and

    ω_±(k) = - 1/2 (i b ±√((4k^2 - b)b)).

Alternatively, we can take our motivation from the gauge driver
construction <cit.> and set ϵ = 1, η an
arbitrary non-negative constant, and obtain

    ω_±(k) = - 1/2 (i η±√(4k^2 - η^2)).

We also considered =η^2 and arbitrary η for uniform
scaling of time.

The bottom line is that the coefficients in (<ref>)
allow us to adjust the size of the critical length scale. For example,
for ϵ = 1 and η arbitrary the critical parameters change
to k_crit = η/2, while for ϵ = 1/b, η=1
we find k_crit = √(b)/2. The damping for large k
becomes e^-η/2 and e^-b/2, respectively.

To avoid small k that drop below (or too far below) k_crit, we
can adjust b or η such that the length scale of k corresponds
to the physical size L of the domain, say λ_crit = 2
L. This will slow down the convergence for large k, but will also
avoid the severe slow down when the damping approaches that of
parabolic relaxation.




 §.§ Hyperbolicity Analysis



In this section we introduce a notation for inverse tensors. These
have to be understood as the inverse of matrices with respect to the
field indices. For example the inverse of the tensor
c^β_α is (c^-1)^α_γ and
we have

    c^β_α (c^-1)^α_γ =
     δ^β_γ  .


For the hyperbolicity analysis we start with writing the hyperbolic
relaxation system in matrix form

    ∂_t  u =  P^k ∂_k  u + h(x^i, u)

with

    P^k =
    [         0 a^k^i^β_α; b^k_i^β_α         0 ]  ,
           u = 
    [  ψ_α; r_iα ]  .

The principal symbol of this system is then given by

    P^s =  P^k s_k = 
    [         0 a^s^i^β_α; b^s_i^β_α         0 ]  ,

where s_i is an arbitrary unit vector, s^i s_i = 1. Suppose
c^β_α = a^s^i^ϵ_α
b^s_i^β_ϵ has a complete set of eigenvectors
w_γ^α with c^β_α
w_γ^α = σ^α_γ
w_α^β, where σ^α_γ is
diagonal. If furthermore all the eigenvalues, the diagonal
elements of σ^α_γ, are positive then P^s has the
following left eigenvectors

    v^0_j γ   =  ( 0 , δ^α_γδ^i_j - b^s_j^ϵ_γ (
                b^s_s^-1 )^α_ϵ s^i
                 )   ,
     v^±_γ   =  ( w_γ^α, ±
          (ρ^-1)^δ_γ w_δ^ϵ
          a^s^i^α_ϵ )   ,

where ρ^β_γ is the root of σ^β_γ, it
is a positive diagonal tensor with ρ^α_γρ^β_α = σ^β_γ.  Note that of the
eigenvectors v^0_jγ only (D-1)N are linearly
independent, while the v^±_γ are 2 N independent
vectors. If there exists a constant K, independent of s^i, such
that V_2 +  V^-1_2 < K, where V is an,
in general s^i-dependent, square matrix constructed from a linearly independent 
set of the eigenvectors v, then the system is strongly
hyperbolic <cit.>. The characteristic
variables û and their characteristic speeds λ are thus

    û^0_jγ   = r_jγ - b^s_j^ϵ_γ (
         b^s_s^-1 )^α_ϵ s^i
         r_iα  ,         λ^0_iγ = 0   , 
    û^±_γ   = w_γ^αψ_α±
    (ρ^-1)^δ_γ w_δ^ϵ
    a^s^i^α_ϵ r_iα  ,    λ_γ^± = ±ρ^α_γ e_α  ,

where e_α denote the Cartesian basis vectors. Note that of the
characteristic variables û^0_jγ only (D-1)N are
linearly independent, while the û^±_γ are 2 N
independent vectors.  From this we can recover the evolved variables
in terms of the characteristics:

    ψ_α   = 1/2 ( w^-1 )_α^γ (û^+_γ + û^-_γ)   , 
    
    r_iα   = û^0_iα + b^s_i^ϵ_α ( c^-1
    )^β_ϵ
               · (
    (w^-1)_β^δρ^γ_δû^+_γ - û^-_γ/2 - a^s^j^γ_βû^0_j γ )   .


We can use the freedom in the choice of b^j_i^β_α
to impose certain properties on the hyperbolic relaxation system. In
the following we briefly discuss three interesting choices, that
fulfill the restrictions we have set for b.

1. b is the identity. A very easy and natural choice is
b^j_i^β_α = δ^j_i δ^β_α. With
this choice we have c^β_α =
a^s^s^β_α, which has only eigenvalues with
positive real part due to (<ref>). The imaginary part
however can be non-vanishing.  If however
a^s^s^α_γ = a^s^s_γ^α, then
c is guaranteed to have a complete set of eigenvectors with purely
real eigenvalues and thus system is strongly hyperbolic.  If we have
a^i^j^α_γ = a^i^j_γ^α, then
the system is even symmetric hyperbolic with symmetrizer:

    H = 
    [          δ^α_γ              0;              0 a^i^j^α_γδ_j l ]  .


2. b is the transpose of a. We can also make the system
trivially symmetric hyperbolic by choosing
b^j_i^β_α = a^j_i_α^β.  The
principal symbol of this system is symmetric and thus the system is
symmetric hyperbolic.

3. b is the inverse of a. We can choose b to be the
inverse of a in the sense that b fulfills
a_k^i^α_γ b^j_i^β_α =
δ_k^jδ^β_γ. This choice is particularly
interesting because we then have c^β_α =
δ^β_α and thus all the non-zero characteristic speeds
have values ± 1. Furthermore the eigenvectors of c become
trivial: w_γ^α = δ^α_γ.  A
symmetrizer for this system is

    H =
    [               δ^α_γ                   0;                   0 a_m^i^α_ω a^m_l_γ^ω ]  .

Since all traveling characteristic variables have the same speeds, we
consider this the best choice for b^j_i^β_α.  A
straightforward generalization of this choice allows
b^j_i^β_α to be scaled by a constant factor,
which will also uniformly scale the characteristic speeds.




 §.§ Boundary Conditions




The basic idea to impose boundary conditions in our method is to
modify the right hand side of the hyperbolic relaxation system. The
outward pointing unit normal covector s_i to the boundary surface is
naturally defined by taking the gradient of a scalar field which is
increasing across, but constant in the boundary, and then normalizing
this gradient to unit magnitude using our arbitrary but fixed
metric. This metric is subsequently used to raise the index and form
the outward pointing vector s^i. We restrict our attention to
strongly hyperbolic systems, for which a regular (s^i-dependent)
similarity transformation matrix T_s exists which transforms
between our evolved variables u and a linearly independent set
of the characteristic variables û given in
Eq. (<ref>) and (<ref>)

    u =  T_s û  .

We can then decompose our evolution
equations (<ref>) as

    ∂_t  u =  P^s ∂_s  u +  P^k q_k^i
    ∂_i  u + h(x^i, u)   ,

where q_k^i = δ_k^i - s_k s^i is the projector onto the
boundary surface. We now multiply by T_s^-1 and obtain

    d_t û =    T_s^-1 P^s ∂_s  u
    +  T_s^-1 ( P^k q_k^i ∂_i  u + h(x^i,
      u) ) 
    
    =    T_s^-1 P^s T_s  T_s^-1∂_s  u + T_s^-1 ( P^k q_k^i ∂_i  u + h(x^i,
      u) ) 
    
    =   Λ^s  d_s û +  T_s^-1
    ( P^k q_k^i ∂_i  u + h(x^i, u) )   .

Here the straight derivative symbol d denotes that the
transformation matrix stands outside of the derivative, d_i
û≡ T_s^-1∂_i  u, and Λ^s is a diagonal matrix containing the characteristic
speeds.  We can now impose boundary conditions on the incoming
variables, those with positive characteristic speeds, by modifying
their right hand sides.  After the right hand sides have been modified
we transform the system back by multiplying with T_s.




  §.§.§ Penalty Method




In the penalty method <cit.> the
boundary conditions are weakly imposed by modifying the right
hand sides of the incoming characteristic variables û^+_γ
in the following way

    d_t û^+_γ=̂ D_t û^+_γ + p (û^+
       BC_γ - û^+_γ)   ,

where p is the penalty parameter, u^+  BC_γ is some
given boundary data that we want û^+_γ to approach, D_t
û^+_γ is the unmodified right hand side and =̂
denotes equality at the boundary. The penalty parameter can not be
chosen arbitrarily and we refer the reader to <cit.> for a
detailed derivation of the penalty parameters used in .




  §.§.§ Maximally Dissipative Boundary Conditions

 


Maximally dissipative boundary
conditions <cit.> will allow us to
set boundary conditions of the form

    s^i ∂_i ψ_α|_∂Ω =
    ϕ_α(ψ_β, ∂_i ψ_β, q^i_j ∂_i
    ∂_k ψ_β )   ,

where ϕ is a function that is allowed to depend on the
coordinates x^i, the fields ψ_α, their derivatives and the
transverse projections (q^i_j = δ^i_j - s^i s_j) of their
second derivatives.  For brevity we suppress dependence on all the
arguments in the following.

Maximally dissipative boundary conditions are imposed by requiring

    (ρ^-1)^β_γ w_β^α∂_t
    ψ_α + w_γ^α∂_s ψ_α =
    w_γ^αϕ_α  .

This boundary condition is actually be different from
Eq. (<ref>) during the relaxation process. However,
again we are only interested in the steady state at the end of the
evolution, where we have ∂_t ψ_α = 0 and thus the
correct boundary condition will be imposed. For numerical stability
the functions ϕ_α must not depend on normal derivatives of
the evolved variables. Therefore in (<ref>) in the
arguments of ϕ_α we have to make the replacements
∂_i ψ_β→ r_i β and q^i_j ∂_i
∂_k ψ_β→ q^i_j ∂_i r_k β.  For
the normal derivatives of the incoming characteristic we obtain the
relation

    d_s û^+_γ   = w_γ^α∂_s
    ψ_α + (ρ^-1)^δ_γ w_δ^ϵ
    a^s^i^α_ϵ∂_s r_iα
       = w_γ^αϕ_α -
    (ρ^-1)^δ_γ w_δ^α ( ∂_t
    ψ_α - a^s^j^ϵ_α∂_s r_j
      ϵ )   ,

where in the actual implementation ∂_t ψ_α is to be
replaced by Eq. (<ref>).  This equation is now
used to impose the boundary condition by replacing the d_s
û^+_γterms in Eq. (<ref>),
yielding the modified right hand side

    d_t û^+_γ=̂    D_t û^+_γ
        - ρ^β_γ ( w_β^α∂_s
    ψ_α + (ρ^-1)^δ_β w_δ^α∂_t ψ_α - w_β^αϕ_α )   .

With the general expression at hand, we now discuss the implementation
of typical boundary conditions.

1. Dirichlet conditions. Dirichlet conditions are of the form
ψ_α |_∂Ω = g_α, where the g_α are
some function defined on the domain boundary ∂Ω. To
implement such a boundary condition ϕ_α has to take the form

    ϕ_α = s^i r_i α + e^β_α ( g_β - ψ_β)   ,

where e is positive definite, e^β_α t_β
t^α > 0.  In the steady state we have ∂_i ψ_α =
r_i α and thus Eq. (<ref>) becomes 0 =
e^β_α ( g_β - ψ_β), which is only
fulfilled for the requested boundary condition.  The positive
definiteness of e is important to guarantee stability at the
boundary. Suppose we have ∂_i ψ_α = r_i α
fixed, then Eq. (<ref>) has the form

    (ρ^-1)^β_γ w_β^α∂_t
    ψ_α = w_γ^α e^β_α (
    g_β - ψ_β)   ,

which would have solutions not asymptoting to g_α if e was
not positive definite.  Besides positive definiteness there are no
further restrictions apparent on e and therefore, it can be chosen
to be the identity e^β_α = δ^β_α,
which we use in our applications.

2. Neumann conditions. Neumann boundary conditions are of the
form s^i ∂_i ψ_α|_∂Ω = g_α.
Their implementation in our method is straightforward; one just has to
take ϕ_α = g_α.

3. Robin conditions. Robin boundary conditions are mixture of
Dirichlet and Neumann boundary conditions and can be written as s^i
∂_i ψ_α|_∂Ω = g_α + f^β_αψ_β, where the f^β_α are functions defined on the
domain boundary.  Their implementation is also straight forward
choosing ϕ_α = g_α + f^β_αψ_β.




§ NUMERICAL SETUP






 §.§ Grid Setup



We employ the pseudospectral hyperbolic evolution code 
and refer the reader to <cit.>, where the grid setup is
explained in detail.  Here we only give a short summary of the basic
grid setup and numerical method.  Our grid consists of different
coordinate patches, a cube patch in the center, transition shell
patches and outer shell patches. On each patch we have a mapping
between local Cartesian coordinates to global Cartesian coordinates,
where on shell patches we employ the “cubed sphere”
construction <cit.>.  The patches themselves can consist
of smaller subpatches, which are the smallest units we use for our
parallelization scheme. On each subpatch we approximate the fields by
a Chebyshev pseudospectral method, the subpatches are discretized
in every direction by the Gauss-Lobatto collocation points. It is then
possible to reconstruct the Chebyshev coefficients from the fields
values at the collocation points.

The  code is adapted to evolutions in three dimensions.
For axisymmetric and spherically symmetric problems we use the cartoon
method to reduce the computational domain to two or one dimensions
respectively <cit.>.




 §.§ Integration Method



The time integration for relaxation methods does not require a high
order of error convergence, since we are only interested in the steady
state at the end of the evolution.  More important are the efficiency
and stability of the integration algorithm.  For the time integration
we use the method of lines.  It is known that for linear hyperbolic
equations the simple forward Euler-method and also explicit
second-order Runge-Kutta methods are unstable <cit.> and
thus are not suited for the integration of the hyperbolic relaxation
equations.

In the applications presented in this paper we employ the popular
fourth-order Runge-Kutta scheme (RK4), which is stable for hyperbolic
equations.  This method needs four evaluations of the right-hand side
per time step, which appears to be not very efficient. After all, we
do not need a very accurate integrator, since we are only interested
in approaching the stationary state.  Therefore it is worthwhile to
investigate other classes of integrators, multistep methods like
the third- or fourth-order Adams schemes <cit.>, which
effectively only require one or two evaluations per time step and are
usually also stable for hyperbolic PDEs. Some simple experiments with
the Poisson equation indicated, however, that RK4 is more efficient
than RK3 or a fourth order Adams scheme since RK4 allows comparatively
large time steps.

Contrary to what is described in <cit.>, we neither use
nor need filtering to assure stability in the hyperbolic relaxation
method, since the system usually tends towards a stable static or
stationary solution automatically.

In our code we have two types of boundaries.  On the one hand we have
boundaries between different subpatches, and on the other hand the
boundaries of our computational domain, in particular the outer
boundaries.  To treat boundaries between subpatches we employ the
penalty method described in
Sec. <ref> setting the
boundary data to be the outgoing characteristic of the neighboring
subpatch.

Treating the outer boundary with this method we have to provide a
function g_γ equaling u^+  BC_γ at the boundary, 
    û^+  BC_γ = w_γ^αψ_α 
    + (ρ^-1)^δ_γ w_δ^ϵ
    a^s^i^α_ϵ r_iα  = g_γ  .

Such a boundary condition however is very unusual in practice for
elliptic equations and thus the penalty method is not suited for the
treatment of our outer boundaries. Instead we only use the maximally
dissipative boundary conditions described in
Sec. <ref>.




 §.§ Initial Guesses



To start the hyperbolic relaxation one has to provide an initial guess
to the solution. A suitable initial guess will always depend on the
specific form of the problem, in particular it should be chosen such
that in the course of the relaxation the variables do not have to
cross any points where the equations (terms in the non-principal
part) become singular.  In our tests we found the solver to be
particularly well behaved when starting with a guess that is
stationary in the interior, but not at the boundary.  The whole
solution then starts to relax from the boundary to the inside.  For
example in our applications to numerical relativity initial data,
taking the flat metric everywhere lead to stable relaxations, which
demonstrates a remarkably high robustness that can be achieved by the
method.  For the reduction variables we simply take the numerical
derivative of the initial guess, ^( ini)r_i α =
∂_i ^( ini)ψ_α




 §.§ Refinement Strategy



To speed up the relaxation process we employ a simple scheme of
successive refinement.  It can be assumed that the right-hand side of
the solution variables ∂_t ψ_α,
Eq. (<ref>), is a good approximation to the
residual of the elliptic equation (Lψ)_α,
Eq. (<ref>).  This however is only true until a
discretization limit is reached below which the norm of the residual
is no longer decreasing. The norm of the ∂_t ψ_α will
typically continue decreasing until machine precision is reached. This
makes it possible to construct an indicator signaling when the
discretization limit is reached and thus relaxation should be
continued on a higher resolution grid.  In particular we choose the
following criterion,

    ∫∑_α = 1^N | ∂_t ψ_α | dV
     < c ∫∑_α = 1^N | (Lψ)_α | dV   ,

where c is some constant smaller than one. For our applications we
found c = 0.1 to be a choice working reasonably well. We note
however that depending on the specific problem also smaller values
might be beneficial.  Additionally we increase the resolution when the
error of the elliptic equation reaches machine precision, when the
norm of (Lψ)_α is smaller than 10^-13 times the number
of grid points.

We start with relaxing the system on the coarsest grid and check every
1000 relaxation steps whether to proceed relaxation on a finer grid
based on whether one of the two criteria mentioned above is fulfilled.
The final resolution can be determined by an error bound on the
residual of the elliptic equation, or by some predetermined
resolution, which may be required for the evolution of the data.  For
the refinement we increase the resolution on every subpatch by two
collocation points in every direction, which is equivalent to adding
two Chebyshev modes in every direction.  We interpolate the coarse
steady state solution to the new subpatches and repeat the procedure
until we arrive at the desired resolution.  We also find it advisable
to use the interpolated values for the reduction variables instead of
taking the numerical derivative of the solution variables, since the
latter introduces new errors, which costs some extra effort to damp.




§ APPLICATION TO TEST CASES






 §.§ Poisson Equation – Finite Differencing 




To provide a reference point independent of the specific
pseudospectral methods of , we first discuss a minimal
implementation using a finite difference method to solve the Poisson
equation.  We consider the hyperbolic relaxation equation

    ∂_t^2 ϕ + η∂_t ϕ = Δϕ - ρ,

which we implement as a first order in time, second order in space system,

    ∂_t ϕ   = π - ηϕ,
    ∂_t π   = Δϕ - ρ.

We consider the fully first order version of this system in
Sec. <ref>. At the boundaries we use asymptotic
Dirichlet conditions analogous to (<ref>),
∂_t ϕ = g - ϕ and ∂_t π = g - π.

We choose centered, second order accurate finite differences in space,
and the default time integrator is the classic fourth-order
Runge-Kutta method.  The numerical domain is an equidistant grid of
points in [-L/2,L/2]^d, dimension d=1, 2, or
3, with Cartesian coordinates. There are N points in each of up to
three directions with a total of V=N^d points.

Let us discuss some results for vanishing source term, ρ = 0, and
vanishing Dirichlet boundary, g=0, where the method has to reduce an
initial guess of ϕ = 1 / (1+x_j x^j) and π = 0 at t=0 to
the asymptotic, late-time value ϕ = π = 0.





In Fig. <ref>, we show results for a box of size
L=20, damping parameter η=0, varying the number of points and
the number of dimensions. The norm is weighted by the grid spacing
Δ x to represent the integral of the residual, |f|_2 = (∑
f^2 Δ x^d)^1/2. Convergence is exponential in time, with two
distinct phases. Inspection of the evolution of ϕ and π shows
that the initial phase corresponds to the damping of short wavelengths
(in this example until t≈ 20), after which long wavelengths
dominate and the convergence is slower.  The convergence of the
(weighted) norm of the residual with time is quite independent of the
resolution.  In this example the time-step is Δ t = λΔ x for a fixed Courant factor λ, so the number of time
steps is proportional to the number N of grid points in one
direction. The work per right-hand-side evaluation is O(N^d), so the
total work to reach a final time T is O(N^d+1).





A key question is how efficient hyperbolic relaxation is compared to
other methods. In Fig. <ref>, we show a comparison of
different methods for a two-dimensional example with N = 101 points.
The methods considered are hyperbolic relaxation as above, the
standard Jacobi iteration <cit.>, and the BiCGSTAB method
as an example for a Krylov subspace method <cit.>. Also
included are two additional variants of hyperbolic relaxation. In
these examples Δ t = 1.0 Δ x for RK4 in 2d.

Referring to Fig. <ref>, the Jacobi method shows the
slowest convergence. Reducing the residuum of the 2d Poisson equation
by a factor 10^-p requires n ≈1/2 p N^2
iterations on a N× N grid <cit.>. For a 2d grid
with V=N^2 degrees of freedom, the operation count is therefore
O(V^2)=O(N^4), compared to O(N^3) for optimal SOR and O(Vlog V)
for multigrid methods. Hyperbolic relaxation with O(V)×
O(N)=O(N^3), as demonstrated in Fig. <ref>, is
therefore a reasonable candidate for further consideration.  In the
concrete example, the Jacobi method is significantly slower than
hyperbolic relaxation, but the Jacobi method is usually not considered
as a stand-alone method.

For this simple comparison, the BiCGSTAB method is used without a
preconditioner, but the Laplace operator leads to a sufficiently well
conditioned operator such that convergence is fast nevertheless,
compared to the other methods considered here. There is an initial
phase of relatively slow convergence, but once the trial solution is
sufficiently close to the final answer, convergence becomes much
faster.

Remarkably, hyperbolic relaxation does about as well as BiCGSTAB
during the first phase. However, convergence slows down after the
shorter wavelengths have been damped and errors due to larger
wavelength remain. We have considered three ideas to improve the
convergence of hyperbolic relaxation for long wavelengths.  Not shown
here is the multi-level refinement strategy which we employ in the
 code, see Sec. <ref>.

As an immediate application of the mode analysis of
Sec. <ref>, we introduced the damping parameter
η, which for the basic experiments so far was set to
η=1. Also shown in Fig. <ref> is the result for
η=0.4, which is slower than η=1 initially, but faster for
later iterations. This effect is related to the size L of the
box. It seems possible to construct a dynamically adjusted damping
η(t). Similar results hold for the parameter b in
(<ref>). The velocity associated with the largest wave
number k scales with √(b) but is independent of η, so for
optimal performance the Courant factor has to be adjusted for the
version with b but can be kept constant for the version with η.

We also experimented with a “one-step overrelaxation” method (as
opposed to successive overrelaxation).  This is based on the
observation that after the initial propagation/damping phase of
hyperbolic relaxation, the second time derivative of ϕ becomes
significantly smaller than the first time derivative,
∂_t^2ϕ≪∂_tϕ. Hence it seems promising to
attempt a linear extrapolation in time. The curve labeled
“overstep” in Fig. <ref> is obtained by searching
every few iterations for the time step Δ T = λΔ t
that minimizes the global residual of ϕ_new = ϕ + Δ T  
F(ϕ), where F is the update suggested by the time stepping
algorithm (RK4).  This is similar to various other 1d step-size
optimizations. For the example considered here (but also for
ρ≠0 as below), the late time solution of hyperbolic relaxation
is sufficiently regular that indeed an appropriate global Δ T
can be found. The overstep algorithm only accepts improvements by a
given factor, say 10 (we tried 2 to 1000). After the adjustment the
solution is disturbed but converges again with the typical speed for
shorter wavelengths to a new regular state, so in the optimal case the
overall convergence rate approaches that of the fast phase of
hyperbolic relaxation.

The main points regarding the convergence rate of hyperbolic
relaxation as shown in Fig. <ref> are that the method
works out-of-the-box and that its performance falls somewhere between
Jacobi and BiCGSTAB. There seems to be quite some potential for
accelerating the convergence rate of hyperbolic relaxation. From the
point of view of solving elliptic equations with a code designed for
hyperbolic equations, note that hyperbolic relaxation is “only”
slower by a factor of about 5 (to reach a residual of 10^-9 in
this example) than a standard method like BiCGSTAB, which however may
not be readily available.




 §.§ Poisson Equation – Pseudospectral Method




To test the hyperbolic relaxation elliptic solver we start by solving
Poisson's equation,

    Δψ - ρ = 0  ,

in spherical symmetry, ρ = ρ(r), r = √(x^i x_i). To
solve this equation we choose the hyperbolic relaxation system

    ∂_t ψ   = δ^ij∂_i r_j - ρ  ,
    ∂_t r_i     = ∂_i ψ - r_i   .

For our first test we take ρ to be smooth, it is infinitely
often continuous differentiable,

    ρ = ρ_0  ( -6/R^2 + 4 r^2/R^4 )
    e^-r^2/R^2  ,

where R and ρ_0 are non-zero parameters. For this ρ
Poisson's equation has the solution

    ψ_ analytic = ρ_0 e^-r^2/R^2  .

At the boundary a falloff in ψ compatible with this solution is
obtained by imposing the Robin boundary condition ∂_r ψ =
s^i ∂_i ψ = - 2 rψ/R^2.

For our second test we take a non-smooth ρ that corresponds to a
homogeneously charged sphere, which is like a toy model for stars.
The density ρ is then given by

    ρ = 
    ρ_0     if  r ≤ R 
     
    0      if  r > R   ,

for which the Poisson equation has the solution

    ψ_ analytic = ρ_0
    r^2/6 - R^2/2    if  r ≤ R 
     
      -R^3/3 r    if  r > R   .

Again we impose Robin boundary conditions according to the falloff of
this solution, ∂_r ψ = s^i ∂_i ψ = -ψ/r.

In our tests we place the outer boundary at radius of 10 and and we
divide the grid into a total of eight subpatches, where the inner five
extend over the interval [0,5] and the outer three, having a coarser
resolution, extend over [5,10].  The parameters determining ρ
are chosen to be R = 5 and ρ_0 = 1.  For the non-smooth case
special care has to be taken to ensure convergence.  In particular we
chose the grid such that the discontinuity lies at a boundary of
subpatch, ensuring second order convergence.  In both test cases the
relaxed solution converges with the number of grid points to the
analytical solution. To investigate the convergence we have to make
sure that the solution is completely relaxed on every resolution.
This is achieved by choosing in Eq. (<ref>) c =
0.0001.  In Fig. <ref> we report the absolute
difference between the analytical and numerical solution integrated
over the outermost subpatch. We note however that the convergence
behavior is the same on all other subpatches.  As expected we find the
error of the numerical solution to decrease exponentially with the
number of points for the smooth ρ from
Eq. (<ref>). For the non-smooth ρ of
Eq. (<ref>), we only get a convergence order of
approximately two, which is the expected convergence order for
discontinuous ρ.  Of course this is not very efficient for a
spectral method. For non-smooth right-hand sides it might be
preferable to increase the number of subpatches (h-refinement) instead
of the number of collocation points per subpatch.





In Fig. <ref> we investigate how the L1-norm of
different quantities, that can be used to approximate the error,
progresses during the relaxation process.  First we observe that the
difference to the analytical solution decreases even when the computed
residual, given by left-hand side of Eq. <ref>, is already
leveling off. This is especially remarkable for the non-smooth case,
where the residual itself is not converging at all.  For the smooth
case we secondly observe that after refining the grid the norm of
right-hand side of Eq. (<ref>) practically
continues at the same level as before. The norm of the residual on the
other hand drops quickly after refining, reaching the right-hand sides
level until again the discretization limit is reached.  These
observations suggest that for problems with smooth solutions it is
preferable to relax for longer on the coarse grid.  For problems with
non-smooth solutions, however, more new error develops during each
refinement and thus refining for longer on the coarse grid is not
paying off. Furthermore, it is preferable to increase the grid
resolution faster.
 




As a last simple test, we investigated the behavior in the case of
non-unique solutions. For this we took the smooth ρ from
Eq. (<ref>) and imposed the Neumann boundary
condition ∂_r ψ = 0, for which multiple solutions
differing only by an additive constant exist. We find that after some
relaxation the right hand side of Eq. (<ref>)
becomes approximately constant in space. From this point on the
solution is no longer improving, since only constant terms, which do
not improve the residual of Eq. (<ref>), are added.




§ APPLICATION TO INITIAL DATA FOR GENERAL RELATIVITY







 §.§ The Extended Conformal Thin-Sandwich Equations




In numerical relativity one usually decomposes the spacetime metric
g_ab into a temporal and spatial part in the form

    g_ab d x^a dx^b = -α^2 d t^2 + γ_ij (d x^i + β^i d
    t) (d x^j + β^j d t)   ,

where α is called the lapse, β^i the shift and
γ_ij the spatial metric.  The equations of motion in general
relativity are subject to constraint equations, which have to be
solved before the spacetime is evolved numerically.  A popular
formulation of the constraint equations is given by the extended
conformal thin-sandwich (XCTS) equations <cit.> in
which the spatial metric is decomposed into a conformal factor ψ
and a spatial conformal metric γ̅_ij as γ_ij =
ψ^4 γ̅_ij.  In the XCTS framework the constraint
equations take the form

    D̅^j D̅_j ψ =   ψ/8R̅ - ψ^5  ( 2
      πρ - K^2/12 + 1/8 A_ij A^ij )  ,
    D̅^j D̅_j β^i =    - 1/3D̅^i D̅_j β^j
      - R̅^i_jβ^j + 16 παψ^4 J^i 
        + (
      D̅^i β^j + D̅^j β^i - 2/3γ̅^ijD̅_k β^k ) 
         ·D̅_j ln(αψ^-6)
      
        - αψ^-6D̅_j ( α^-1ψ^6 ∂_t
      γ̅^ij ) + 4/3αD̅^i K   ,
     
    D̅^j D̅_j ( αψ ) =   αψ^5 
      (7/8 A_ij A^ij + 5/12 K^2 + 2 π (ρ + 2 J)
       ) 
        - ψ^5 (∂_t - β^j D̅_j ) K +
      1/8αψR̅  .

Here D̅_i is the covariant derivative compatible with the
conformal metric γ̅_ij, R̅_i_j is the
Ricci tensor of γ̅_ij and R̅ is the corresponding
Ricci scalar.  The tensor A_ij is the tracefree part of the
extrinsic curvature K_ij and K is the trace of K_ij.  The
matter source terms are defined as the following contractions of the
energy-momentum tensor T_ab: ρ = T_ab n^a n^b, J^i =
-T_abγ^i a n^b and J = γ^ab T_ab, where n^a is
the timelike vector normal to the spatial hypersurface n^a =
(1/α,-β^i/α). For conformal quantities (denoted by a
bar) the conformal spatial metric γ̅^ij lowers and raises
indices and for unbarred quantities the physical spatial metric
γ^ij is used.  In the XCTS equations γ̅_ij,
∂_t γ̅_ij, K and ∂_t K are given
functions, depending on the type of initial data you want to
construct.

In Eq. (<ref>) the product αψ is taken as
one variable. For our computations we rewrite this equation with the
help of Eq. (<ref>) as

    D̅^j D̅_j α =    -2/ψ ( D̅^j α) (D̅_j ψ) - ψ^4 (∂_t - β^j D̅_j ) K 
        + αψ^4  (A_ij A^ij + K^2/3 + 4 π (ρ + J)
       )   ,

so we can solve directly for α.  We have tested our elliptic
solver with both versions and found them to work equally well in our
applications. In the following we investigate the XCTS system with our
replacement for the lapse equation, but the analysis would be exactly
the same for the original system.

The principal part of the XCTS equations is given by

    [                       γ̅^kj                           0                           0;                           0 γ̅^kjδ^l_i + 1/3γ̅^l kδ^j_i                           0;                           0                           0                       γ̅^kj;                             ]∂_k ∂_j 
      [   ψ; β^i;   α;     ]  .

The principal part is coupled only between the components of
β^i. Thus we can carry out the hyperbolicity analysis
independently for ψ, α and β^i.

The metric to lower and raise spatial indices in the hyperbolic
relaxation method is in principal arbitrary, but for the XCTS
equations we use the conformal metric γ̅_ij, because this
choice simplifies the following formulas considerably. Another
peculiarity is the fact that spatial indices and field indices
“mix”, but in general they are lowered and raised by different
metrics. Therefore we do not lower and raise field indices and instead
write the Euclidian metric explicitly.

For the conformal factor part we have in the hyperbolic relaxation
system a^ij = γ̅^ij and we choose b^i_j =
δ^i_j. Here we have suppressed the field indices, because we
only have a single field ψ.  The characteristic variables and
speeds of the hyperbolic relaxation system of the ψ part are then
given by

    û^0_j    =  r_j - s_j s^i r_i^(ψ)  , 
               λ^0_i = 0   , 
    û^±   = ψ± s^i r_i^(ψ)  ,
                  λ^± = ± 1   ,

where r_i^(ψ) is the reduction variable for ψ.

The lapse part has exactly the same principal part, so we have
identical a^ij and b^i_j and the characteristics are

    û^0_j    =  r_j - s_j  s^i r_i^(α)  ,
               λ^0_i = 0   , 
    û^±   = α±s̅^i r_i^(α)  ,
                  λ^± = ± 1   ,

with r_i^(α) being the reduction variable for α.

In the derivation of the hyperbolic relaxation equations we labeled
the fields with lower indices. For the shift however we use upper
indices here and thus one has to be careful, not to confuse the
indices. Therefore, we introduce auxiliary fields ϕ_α with

    β^α = δ^αβϕ_β  .

Substituting β^α we obtain for the principal part of the 
shift equations

    a^i^j^β_α∂_i ∂_j ϕ_β = (
     γ̅^i jδ^β_α + 1/3γ̅^i
       ϵδ_ϵαδ^j β ) ∂_i
     ∂_j ϕ_β  .

We take b^i_j^β_α to be the inverse of a (as
defined in section <ref>),
b^i_j^β_α = δ^i_j δ^β_α -
1/6δ^i βδ_j α .  The characteristic
variables and speeds for this part of the hyperbolic relaxation system
are then given by

    û^0_jγ   = r_jγ - b^s_j^ϵ_γ (
         b^s_s^-1 )^α_ϵ s^i
         r_iα
       = r_jγ - s_j s^i r_i γ +
         1/5 (δ_j γ - s_j s^l δ_l γ) s_k δ^k α s^i r_i
           α  ,
    λ^0_iγ   = 0  , 
    û^±_γ   = ϕ_γ± a^s^i^α_γ
    r_iα      
       = ϕ_γ± (s^i r_iγ +
    1/3 s^ϵδ_ϵγδ^j α r_j
      α)   , 
    λ_γ^±   = ± 1   ,

where the r_i_α denote the reduction variables for
the auxiliary fields ϕ_α

At the domain boundary we want the solution to fall off like the
Schwarzschild solution, ψ = a/r + 1 and α =
b/r + 1.  This ansatz gives rise to the following Robin
boundary conditions

    ∂_s ψ |_∂Ω = 1-ψ/r  ,       ∂_s α |_∂Ω = 1-α/r  .

For the shift we likewise impose a radial falloff by the Robin
condition

    ∂_s β^i |_∂Ω = β^i/r  .


As an initial guess we always use the flat space solution, ψ =
1, α=1, β^i = 0.  Of course an initial guess, that is a
good approximation to the solution is always the preferred start for
the relaxation, since it will take less time to relax to the solution
or might be necessary to relax at all. However we find our simple
initial guess to work well and it demonstrates in a nice way the high
robustness of the hyperbolic relaxation method exhibited in our
experiments.




 §.§ Scalar Field



The energy-momentum tensor for a scalar field ϕ is given by 

    T_ab = ∇_a ϕ∇_b ϕ - 1/2 g_ab (
     ∇_c ϕ∇^c ϕ + m^2 ϕ^2  )   ,

where ∇ denotes the covariant derivative compatible with
g_ab.  We consider conformally flat moment-of-time-symmetry
initial data, γ̅_ij = δ_ij, n^a ∇_a ϕ
= 0, ∂_t γ̅_ij = 0 and maximal slicing K=0,
∂_t K = 0. This yields for the matter quantities

    ρ   = 1/2γ^ij (∂_i ϕ) (∂_j ϕ) +
     1/2 m^2 ϕ^2   , 
    
     J^i    = 0   , 
    
     J    = -1/2γ^ij (∂_i ϕ) (∂_j ϕ) -
     3/2 m^2 ϕ^2   .

For a massless scalar field (m=0) the solutions for lapse and shift
are trivially given by α = 1 and β^i = 0 and we only have
to determine the conformal factor ψ by solving
Eq. (<ref>), now taking the form

    0 = δ^ij∂_i ∂_j ψ + πψδ^ij
     (∂_i ϕ) (∂_j ϕ)   .

For a massive scalar field (m ≠ 0) one would have either have to
solve additionally for the lapse or one gives up the requirement on
∂_t K.  For the scalar field we choose radially symmetric,
smooth initial data of the form

    ϕ(r) = p (tanhr/σ - tanhr/σ)    ,

where σ and p are free parameters, that we choose for our
test to be σ = 1 and p=0.1.

The computational domain is divided into eight subpatches, where the
inner five subpatches are of smaller extent to improve the resolution
near the center. In Fig. <ref> we show the
numerical solution for the conformal factor and the absolute value of
the residual of Eq. (<ref>) for a resolution of 21
collocation points per subpatch.  A general feature that can be
observed in the residual are spikes at the boundary of subpatches,
which are expected for quantities involving first and second
derivatives in a discontinuous Galerkin approximation.

We also investigate how the solution converges with increasing
resolution. For this purpose we look at the Chebyshev coefficients and
investigate their convergence against the coefficients of a high
resolution solution. In Fig. <ref> we present the
convergence behavior of the lowest Chebyshev mode C(0,0,0) against
its value for a resolution of 31 collocation points. We observe an
exponential convergence until we hit machine precision at around a
resolution of 25 collocation points per subpatch.












 §.§ Tolman-Oppenheimer-Volkoff Star



The energy-momentum tensor of a perfect fluid is given by

    T_ab = (ϵ + p) u_a u_b + p g_ab  ,

where ϵ is the proper energy density, p the fluid pressure
and u_a the fluid four-velocity.  The Tolman-Oppenheimer-Volkoff
(TOV) solution <cit.> is a static
radially symmetric solution to general relativity, thus we have u_a =
(α, 0), γ̅_ij = δ_ij, ∂_t γ̅_ij = 0 and we consider again maximal slicing K=0,
∂_t K = 0.  The matter quantities then become

    ρ = ϵ  ,        J^i = 0   ,        J = 3 p     .

We assume for our tests a polytropic equation of state and express the
matter quantities in terms of the specific enthalpy h,

    ρ   =  ( 1 + n(h-1)/1+n )  (
     h-1/κ (1+n) )^n   , 
    
     J    = 3 κ (
     h-1/κ (1+n) )^n+1  ,

where κ is the adiabatic constant and n is the polytropic
index.  The Euler equation follows from energy-momentum conservation,
for a fluid with temperature T=0

    u^a ( ∇_a (h u_b) - ∇_b (h u_a) ) = 0   ,

which has to be satisfied in addition to the XCTS equations. For our
assumptions the Euler equation is satisfied for α h = const.
and thus specifying the specific enthalpy at the origin yields

    h(r) = 
     h(0) α(0)/α(r)    if α(r) < h(0) α(0) 
     
        1     else   .

We can immediately get the solution for the shift β^i = 0 and
are left with solving Eq. (<ref>) and (<ref>)
of the XCTS system.

For our test we choose an adiabatic constant of κ = 123.6489
and a polytropic index of n = 1 and the enthalpy in the star's
center is set to h(0) = 1.2.  We present the solution for the
conformal factor and the lapse in Fig. <ref> and
investigate the convergence of their lowest Chebyshev modes in
Fig. <ref>.  The residuals are large at the stellar
surface, where α(r) = h(0) α(0).  This is caused by
the fact that the matter terms are not smooth at this point, which can
be seen in the kink in the specific enthalpy h.  As for the
non-smooth right-hand sides discussed in
Sec. <ref>, we have to make sure that the kink lies
at a subpatch boundary. By trial and error we find the stellar surface
for the above parameters to be located at r = 9.7098 and we place a
subpatch boundary at this position “by hand”.  A more sophisticated
method is to fit coordinates automatically to the stellar
surface <cit.>. This however is beyond the scope of this paper
as we here want to focus on applications of the hyperbolic relaxation
method.  In contrast to what we observed in the non-smooth case for
the Poisson equation, here the Chebyshev coefficients converge
exponentially despite the kink in the specific enthalpy. The
convergence rate however is much smaller than that observed for the
initial data of the scalar field.












 §.§ Neutron Star Binaries



For the construction of neutron star binary initial data we follow the
scheme of <cit.> using the constant three-velocity
approximation. However, to solve the XCTS equations we do not rely on
iterating the solution of the equations for the conformal factor,
lapse and shift separately, but instead we solve for all variables
simultaneously relaxing the complete XCTS system, which accelerates
the solution process.  This is a feature that most solvers for this
type of initial data do not provide, and it could turn out to be an
advantage of our method.  Furthermore, we do not start with superposed
(boosted) TOV solutions, but instead start with a flat metric ψ =
1, α = 1, β^i, as discussed at the end of
Sec. <ref>.

For our test we consider equal mass neutron stars with a specific
enthalpy of h = 1.01 in each of their centers and a separation of
80. For the equation of state we choose again κ = 123.6489 and
n = 1. The stars' centers are located at the z-axis and their
velocities are parallel to the x-axis. We construct in initial data
for irrotational stars on a quasicircular orbit.  In
Fig. <ref> we present results for the conformal
factor, lapse, the x-component of the shift, the residual of the
conformal factor equation.  As for the TOV star initial data, the
residuals are biggest on those subpatches which contain the stellar
surfaces, where the matter fields are not smooth.

Because we are not using surface-fitted
coordinates <cit.> we can not place the
stellar surface at a subpatch boundary, and thus no high-order
convergence in the norm of the residuals can be seen with increasing
number of collocation points.  Although this may be a disadvantage for
studies of initial data per se, the situation changes if the goal is
evolution of the data. Since in an actual evolution of this data
surface-fitted coordinates are normally not retained, the high
accuracy of initial data with surface-fitted coordinates will be lost
relatively quickly anyway. On the other hand, methods like
<cit.> require expensive iterations to determine the surface
fitting coordinates as part of the solution process, so any method
which works without special coordinates is more efficient in that part
of the algorithm. In fact, part of the motivation for the multigrid
method in <cit.> was to construct a solver which works on
a general Cartesian grid without surface-fitting coordinates. The
hyperbolic relaxation method achieves the same goal.








§ CONCLUSIONS



The most common types of relaxation methods are based on the famous
Gauss-Seidel method, which can be motivated by rewriting the problem
as a parabolic diffusion equation that relaxes from an initial guess
to the solution of the elliptic PDE.  In this paper we investigated a
new class of relaxation methods, which are not based on parabolic, but
rather on hyperbolic PDEs. In the literature hyperbolic relaxation is
usually discussed from the point of view that a hyperbolic equation is
given which may contain physical or numerical dissipation terms.  In
this work we assume that an elliptic equation is given, which is
extended to a hyperbolic relaxation equation for the purpose of
solving the elliptic equation. In some respect, hyperbolic relaxation
might actually be as well suited for the solution of elliptic PDEs as
parabolic relaxation.

We investigated how a hyperbolic relaxation method (HypRelax) can be
constructed for a general class of second order, non-linear elliptic
PDEs and discussed its structure and properties.  A discussion of the
general hyperbolicity properties has been carried out and three
specific choices for the relaxation were discussed. For the special
case of the Laplace equation, a mode analysis revealed that there is a
critical wavenumber at which the qualitative behavior of the modes is
changing. It has also been seen that at low wavenumber the damping
rate approaches that of the Jacobi method from above.  It is an
interesting topic for the future how the specific choice of the
relaxation system (choice of b^j_i^β_α) affects
the behavior of the modes, how that choice can be optimized, and what
would be the modes for more general elliptic equations.

Furthermore, we have shown how the standard types of boundary
conditions – Dirichlet, Neumann and Robin boundary conditions – can
be implemented within the hyperbolic relaxation framework.

With regard to advantages, hyperbolic relaxation shares with several
other methods the feature that it is “matrix free”, it is
possible to avoid the construction of an explicit matrix form by
applying the differential operator directly. Also note that as in
Jacobi methods, non-linear equations can be treated without
linearization, avoiding the additional work of outer
Newton-Raphson iterations.

The ease of implementation is a key feature of hyperbolic relaxation.
In the past, we implemented a parallel geometric multigrid method for
the  and 
codes <cit.>.  We also
interfaced  to the 
package <cit.> for access to a variety of
elliptic solvers, in particular algebraic multigrid.  To avoid
complicated dependencies, some of the black hole initial data was
later computed with a stand-alone Newton-Raphson method, or even with
the direct solver in  for numerical
stability <cit.>. Our current production runs for a wide
range of binary neutron star initial data rely on the sophisticated
spectral solver  <cit.>.
Unsurprisingly, compared to the various approaches just mentioned,
hyperbolic relaxation is very straightforward to implement given a
hyperbolic evolution code.

We have implemented the hyperbolic relaxation method in our spectral
hyperbolic evolution code  <cit.> (and
in a finite differencing test code).  For numerical tests we have
applied it to Poisson's equation, and we also presented applications
to numerical relativity initial data, where we have seen a high
robustness with respect to the choice of an initial guess and
regarding the simultaneous solution of the XCTS equations.  Given the
complexity of binary neutron star initial data and the correspondingly
involved numerical methods that have been developed to solve these
elliptic equations, <cit.> and references therein, it is a
non-trivial result that the hyperbolic relaxation method results in a
robust and quite efficient elliptic solver.

We have seen that the damping rate at low wavenumbers is comparable to
that of the Jacobi method and thus very low.  In this study we were
most interested in the basic properties of the hyperbolic relaxation
method and thus we focussed on a simple scheme of successive
refinement, but there exist more sophisticated accelerators, in
particular multigrid methods, that could be implemented also for a
hyperbolic relaxation scheme. Another possible extension of our
hyperbolic relaxation implementation would be an adaptive mesh
refinement scheme, for which we see two main advantages.  First of
all, there are the usual savings due to optimized local resolution,
and secondly, since we want to use the solution of the elliptic
equation as initial data, this will provide us with a grid that should
already be well adapted to the evolution of the obtained data.
Another idea is to employ adaptive time stepping, for which however
special care would be required to keep the numeric scheme stable.

With regard to initial data for neutron stars, one of the potential
problems is the lack of differentiability at the surface of the stars.
We were prepared to evolve the neutron star data with a
high-resolution shock-capturing method, but in fact this was not
necessary given the strong diffusion of the equations. In applications
where the wave propagation feature is important, shock handling may be
a feature that comes at no extra cost assuming that there is an
evolution code providing the appropriate methods.

The investigations of this paper can only serve as a first step in the
exploration of this potentially promising branch of new relaxation
methods. Hyperbolic relaxation might become with further study, and in
particular in combination with acceleration methods, the basis of an
alternative numerical solution method for elliptic equations.





This work was supported in part by the DFG funded Graduiertenkolleg
1523/2, and the FCT (Portugal) IF Program IF/00577/2015. We thank
Niclas Moldenhauer and Tim Dietrich for fruitful discussions on the
construction of neutron star initial data.















