Abstract

We have tested the appropriateness of two-soliton analytic metric to describe
the exterior of all types of neutron stars, no matter what their equation of state or
rotation rate is. The particular analytic solution of the vaccuum Einstein equations
proved quite adjustable to mimic the metric functions of all numerically constructed
neutron-star models that we used as a testbed. The neutron-star models covered
a wide range of stiffness, with regard to the equation of state of
their interior, and all rotation rates up to the maximum possible rotation rate allowed for
each such star. Apart of the metric functions themselves,
we have compared the radius of the innermost stable circular orbit
RISCO, the orbital frequency Ω≡dϕdt of
circular geodesics, and their epicyclic frequencies Ωρ,Ωz, as well as the
change of the energy of circular orbits per logarithmic change of orbital frequency Δ~E.
All these quantities, calculated by means of the two-soliton analytic metric, fitted with good accuracy
the corresponding numerical ones as in previous analogous comparisons
(although previous attempts were restricted to neutron star models with either high or low rotation rates).
We believe that this particular analytic solution could be considered as an analytic faithful representation of the
gravitation field of any rotating neutron star with such accuracy,
that one could explore the interior structure of a neutron star by using this space-time to interpret observations
of astrophysical processes that take place around it.

keywords:

The amount and accuracy of modern observations in various parts of
the electromagnetic spectrum has increased dramatically. In order
to give astrophysically plausible explanations of the various
problems related to the observations we have to rely on
theoretical assumptions that are at least as accurate as the data
we are trying to analyze. There is a large class of observations
(see e.g., van der Klis (2006)) that are related to the astrophysical environment
of compact relativistic objects (AGNs, LMXBs, etc).
Furthermore, the anticipated successful gravitational wave detection will open a new window to
observe such objects. In order to understand these phenomena, one has to have a sufficiently
accurate analytic description of the spacetime around such compact objects.
If the central object is a black hole, there is a unique
choice in the framework of general relativity: the Kerr spacetime.
On the other hand the geometry around a rotating neutron star is
much more complicated, since it depends on many parameters related
to the internal structure of the neutron star and the way it rotates.

The assumption that the geometry around such an object is approximately that of
a Schwarzschild, or a Kerr metric (see for example van der Klis (2006))
is very simplistic and it may lead to erroneous conclusions about the actual
astrophysical processes that take place in the close neighborhood
of the star itself (cf. Pappas (2012); Pachón, Rueda & Valenzuela-Toledo (2012)).

One can alternatively rely on numerical codes that are able to
describe the geometry around a realistic neutron star in a tabular
form on a given grid with sufficiently high accuracy. There are various groups
(see Stergioulas & Friedman (1995), and for an extended list of numerical schemes
see Stergioulas (2003)), that have acquired expertise in building relativistic
models of astrophysical objects with adjustable physical
characteristics and constructing the metric inside and outside
such objects by solving numerically the full Einstein equations in stationary, axisymmetric cases.

Although studying astrophysical phenomena in a geometric
background that has been constructed numerically is plausible,
there are certain drawbacks in using such metrics: (i)
Computing various physical quantities of a system, like the
orbital frequencies, or the innermost circular orbit, from a
metric that is given in a tabular form is not very practical and
is often plagued by numerical errors.
(ii) Astrophysical observations from the environment of a compact
object could be used to read the physical parameters that are
related to the structure of the compact object such as its mass,
its equation of state, its rotation, or to obtain the law of its differential
rotation, etc. This would be very difficult to achieve with a
metric that is numerically constructed. Solving the inverse
problem by means of a numerical metric is a blind process that can not
be easily led by physical insight. Instead an analytic expression for the
corresponding spacetime would be much more preferable.

There are various analytic metrics that have been used in the past
to describe the exterior geometry of a neutron star. As mentioned
above, the Schwarzschild metric is not accurate enough for
rotating neutron stars, while the Kerr metric is good only for a
collapsed object (a black hole) but it fails to describe the
exterior of a neutron star, as comparisons of Kerr with numerical
geometries of rotating neutron stars by Berti & Stergioulas (2004) have shown.

The Hartle-Thorne metric of Hartle & Thorne (1968), which has been constructed as
an approximate solution of the vacuum Einstein equations (VEE)
for the exterior of a slowly rotating star, has been extensively used by
various authors to describe neutron stars of low rotation rate (see for example
Berti et al. (2004)). Finally, various other analytic solutions of VEE have been constructed and some of them have been
used, especially during the last decade, to describe the exterior geometries
of neutron stars (see Stute & Camenzind (2002); Berti & Stergioulas (2004); Pachón, Rueda & Sanabria-Gómez (2006); Pappas (2009); Teichmüller, Fröb & Maucher (2011)).
Such solutions are based on the formalism developed by Ernst (1968a, b) which reformulates
Einstein equations in the case of axisymmetric, stationary
space-times. Manko et al. and Sibgatullin (see the articles of
Sibgatulin (1984); Manko & Sibgatulin (1993); Manko, Martin & Ruiz (1995); Ruiz, Manko & Martin (1995); Manko, Mielke & Sanabria-Gómez (2000))
have used various analytic methods to produce such space-times
parameterized by various parameters that have different physical
context depending on the type of each solution.

Such an analytic solution, with its parameters appropriately adjusted to
match numerical models of neutron stars, could then be used to describe the
stationary properties of the space-time around the neutron star itself;
that is, study the geodesics in the exterior of the neutron star.
More specifically, from the analytic solution we could obtain
bounds of motion for test particles orbiting the neutron star,
find the location of the innermost stable circular orbit (ISCO), compute the orbital frequency of the
circular orbits on the equatorial plane as well as the epicyclic
frequencies around it and perform any sort of dynamical analysis on the
geodesics (for example see Lukes-Gerakopoulos (2012)). These properties of the space-time could be used
to study quantitatively astrophysical phenomena that take place in the vicinity of
neutron stars, such as accretion discs. Inversely one could use the astrophysical observations
related to such phenomena to determine the parameters describing the analytic space-time and from
that acquire information for the central object.

The central issue with analytic metrics is whether one can find
solutions that are able to describe with sufficient faithfulness all kinds
of rotating neutron stars; either slowly or rapidly rotating ones, or even
differentially rotating ones.

One solution that has been recently used by Stute & Camenzind (2002) and later by
Berti & Stergioulas (2004) to describe the exterior space-time of rotating
neutron stars is the three-parameter solution of Manko, Mielke & Sanabria-Gómez (2000) (also mentioned as Manko et al.).
Although this solution was shown to match quite well
the space-time of highly rotating neutron stars, it failed to match
the slowly rotating ones. The reason for failing to describe slow rotation is that in the zero angular momentum
limit, this particular solution has a non vanishing quadrupole
moment, while one would expect slowly rotating neutron stars
to be approximately spherically symmetric. This problem of the Manko et al.
solution was not considered disappointing by Berti and Stergioulas, since
the space-time around slowly rotating stars could be described equally well by the Hartle-Thorne approximation.

The three-parameter solution of Manko et al. is a special case of the
so called two-soliton solution, which was constructed by
Manko, Martin & Ruiz (1995). The two-soliton is a four-parameter analytic metric which,
contrary to the previous one, can be continuously reduced to a
Scwarzschild or a Kerr metric while it does not suffer from the
problematic constraints of the Manko et al. solution with respect to the
anomalous behavior of its quadrupole moment. Actually the first four multiple moments of the
two-soliton solution can be freely chosen. Of course the analytic form of the two-soliton solution is not as
compact as the Manko et al. solution, but this is the price one has to pay in
order to cover the whole range of the physical parameters of a neutron star with a single analytic metric.

In the present work, that constitutes the extension and completion
of preliminary results presented by Pappas (2009), that were suitably
corrected with respect to the right extraction of the multipole moments of the
numerical space-time as was recently demonstrated by Pappas & Apostolatos (2012),
we are using this two-soliton solution to describe the space-time
around a wide range of numerically constructed rotating neutron
stars. We use the numerical multipole moments to set the multipole
moments of the analytic space-time. Then we examine how well the
two metrics match each other. Moreover we have performed comparisons between
astrophysicaly relevant geometric quantities produced from the
numerical and the analytic space-times, like the position of the
innermost stable circular orbit (ISCO), the orbital frequencies,
the epicyclic frequencies that are related to slightly
non-circular and slightly non-equatorial orbits, and the change of
energy of the circular orbits per logarithmic change of the
orbital frequency, Δ~E. The overall picture is that the new metric matches the numerical one
with excellent accuracy for all rotation rates and all equations of state.

The rest of the paper is organized as follows: In Sec. 2
the proposed analytic solution (two-soliton) is briefly presented
and some of its properties are thoroughly analyzed. The parameter space of
the two-soliton is investigated and it is shown how to obtain the
limiting cases of Schwarzschild, Kerr and Manko et al.
A brief discussion of the physical properties of the space-time such as
the presence of singularities, horizons, ergoregions and regions
of closed timelike curves (CTCs) is also given. In Sec. 3 we
show how we match the analytic two-soliton solution to a specific numerical one
by matching the first four multipole moments, and show why this is generally the best choice.
In Sec. 4 we discuss various criteria that could be used to compare the two metrics.
Finally, in Sec. 5 the final comparison criteria and the results of the corresponding
comparisons are presented. In Sec. 6 we give an overview of the conclusions obtained by our study.

The vacuum region of a stationary and axially symmetric space-time can be described by the Papapetrou line element,
which was first used by Papapetrou (1953),

ds2=−f(dt−ωdϕ)2+f−1[e2γ(dρ2+dz2)+ρ2dϕ2],

(1)

where f,ω, and γ are functions of the
Weyl-Papapetrou coordinates (ρ,z). By introducing the complex
potential E(ρ,z)=f(ρ,z)+ıψ(ρ,z),
Ernst (1968a) reformulated the Einstein field equations for this type of
space-times in a concise complex equation

(Re(E))∇2E=∇E⋅∇E.

(2)

The real part of the Ernst potential E is the metric
function f, which is also the norm of the timelike Killing
vector related with stationarity, while ψ is a scalar
potential related with the twist of the same vector.

A general procedure for generating solutions of the Ernst
equations was developed by Sibgatulin (1984); Manko & Sibgatulin (1993); Ruiz, Manko & Martin (1995); Manko, Martin & Ruiz (1995). Each
solution of the Ernst equation is produced from a choice of the Ernst potential along
the axis of symmetry of the metric in the form of a rational function

E(ρ=0,z)=e(z)=P(z)R(z),

(3)

where P(z),R(z) are polynomials of z of order n
with complex coefficients in general. The algorithm developed by
Ruiz, Manko & Martin (1995) works as follows: First the Ernst potential along the axis is expressed in the form

e(z)=1+n∑k=1ekz−βk,

(4)

where βk are the roots of the polynomial R(z) and ek
are complex coefficients appropriately chosen so that the latter
form of e(z) (eq. 4) is equal to the former one (eq. 3).
Subsequently one determines the 2n roots of equation

e(z)+e∗(z)=0,

(5)

where ∗ denotes complex conjugation. These roots are denoted as ξk, with k=1,2,…,2n
and from these one defines the 2n complex functions Rk=√ρ2+(z−ξk)2. All these
functions and roots are then used as building blocks for the following determinants

The Ernst potential and the metric functions are finally expressed in terms
of the determinants given above as:

E(ρ,z)

=

E+E−,

(11)

f(ρ,z)

=

E+E∗−+E∗+E−2E−E∗−,

(12)

e2γ(ρ,z)

=

E+E∗−+E∗+E−2K0K∗0∏2nk=1Rk,

(13)

ω(ρ,z)

=

2I[E−H∗−E∗−G]E+E∗−+E∗+E−.

(14)

We should note that due to the form of the metric functions, the
parameters ek and their complex conjugates e∗k that appear in the
determinants cancel out (the ∏nk=1eke∗k is a common
factor of all products of determinants that show up in the metric
functions), so they do not affect the final expressions.

The vacuum two-soliton solution (proposed by Manko, Martin & Ruiz (1995)) is
a special case of the previous general axisymmetric solution that
is obtained from the ansatz (see also Sotiriou & Pappas (2005))

e(z)=(z−M−ia)(z+ib)−k(z+M−ia)(z+ib)−k,

(15)

where all the parameters M,a,k,b are real. From the
Ernst potential along the axis one can compute the mass and
mass-current moments of this space-time. Particularly, for the
two-soliton space-time the first five non-vanishing moments are:

M0

=

M,M2=−(a2−k)M,

M4

=

[a4−(3a2−2ab+b2)k+k2−17kM2]M

J1

=

aM,J3=−[a3−(2a−b)k]M.

(16)

The mass moments of odd order and the mass-current moments of even
order are zero due to reflection symmetry with respect to the
equatorial plane (z=0) of the space-time (this is actually
ensured by restricting all parameters of eq. (15) to assume real values).
From the moments we see that the parameter M corresponds to the mass monopole of the
space-time, the parameter a is the reduced angular momentum, k
is the deviation of the reduced quadrupole from the corresponding
Kerr quadrupole (the one that has the same M and a),
and b is associated with the deviation of the current octupole moment from the current octupole of the
corresponding Kerr.

Case

d2

κ2+

κ2−

ξ+

ξ−

R+

R−

r+

r−

Ia

>0

>0

>0

R

R

R

R

R

R

Ib∗

>0

>0

=0

R

R

R

R

R

R

=ξ+

=R+

=R−

IIa

>0

>0

<0

C

C

C

C

C

C

=(ξ+)∗

=(R+)∗

=(R−)∗

IIb∗

>0

=0

<0

I

I

C

C

C

C

=(ξ+)∗

=(R+)∗

=R−

=R+

IIc

>0

<0

<0

I

I

C

C

C

C

=(R+)∗

=(r+)∗

III

<0

C

C

R

I

R

R

C

C

=(κ2+)∗

=(r+)∗

IVa∗

=0

>0

=κ2+

R

=0

R

R

R

R

=κ+

=r+

IVb∗

=0

=0

=0

=0

=0

-

-

-

IVc∗

=0

<0

=κ2+

I

=0

C

C

R

R

=κ+

=r+

Table 1: Classification of the various two-soliton solutions depending on the values of the parameters
d,κ±,ξ±,R± and r±. The table also shows the various conjugation relations
between the parameters. The types of solutions that have degeneracies are indicated with an asterisk (∗).
R means real, I means imaginary, and C
means complex.

For the two-soliton ansatz (15), the characteristic equation (5) takes the form,

z4−(M2−a2−b2+2k)z2+(k−ab)2−b2M2=0.

(17)

Since the coefficients of the polynomial are real, the roots can be either real, or conjugate pairs. The four roots of
(17) can be written as,

ξ1=−ξ3=ξ+,ξ2=−ξ4=ξ−,

(18)

where,

ξ±=12(κ+±κ−),

(19)

with

κ±=√M2−a2−b2+2k±2d,

(20)

and

d=√(k−ab)2−b2M2.

(21)

Using these symbols for the four roots we redefine the four corresponding functions Rk as

R±=√ρ2+(z±ξ+)2,r±=√ρ2+(z±ξ−)2.

(22)

Next we proceed to classify the various types of solutions depending on whether the four roots
have real, purely imaginary or complex values. This classification is outlined in Table 1.

This case is characterized by two real roots ξ±.
The Kerr family of solutions, which corresponds to
k=0 is definitely not included in this family of solutions.

This subfamily of case I is the simplest to compute, since all
functions R±,r± are real.

This is a degenerate case where the roots ξ± coincide.
The degeneracy is due to κ− being zero
which corresponds to the parameter constraint M2−a2−b2+2k−2d=0.
In such degenerate cases, the computation of the metric function is not straightforward since
the expressions for the metric become
indeterminate, and a limiting procedure should then be applied. In the
reduced-parameter space (a/M,b/M,k/M2), the previous
constraint corresponds to a two-dimensional surface.

In this case, the roots ξ± are either complex or
imaginary, since κ2−<0. Furthermore this means that there
are nonvanishing values of (ρ,z) where the functions R±,r±
assume zero value, which then leads to singularities at the corresponding points.

This sub-case, as with Case I, belongs to a class of solutions that
cannot have a vanishing parameter k.

Here, a degeneracy shows up again as in Case Ib, which admits
the same treatment (limiting procedure) as in the former situation. Contrary to all previous cases, Case IIb admits
a Kerr solution that belongs to the hyper-extreme branch (|a|>M).
Similarly to Case Ib this solution is also represented
by a two dimensional surface in the reduced-parameter space.

This case is similar to the previous one, without the degeneracy
in the roots ξ±. It also includes hyper-extreme Kerr solutions.

In this case one of the ξ+,ξ− is real while the other one is imaginary.
Thus the same type of singularity issues,
as in Case II arise. In particular such problematic behavior shows up on the z=0 plane.
The Kerr and the Schwarzschild solutions lie entirely within this family of solutions.

All types of solutions belonging to this case are degenerate
(there is a special constraint between the parameters) and
as such are probably of no interest to realistic neutron stars. Subcases IVa and IVc have one double root
(ξ−=0), while subcase IVb has a quadruple root (ξ+=ξ−=0)
and the computation of the metric functions needs special treatment. We should also note that
cases IVb and IVc include the extreme Kerr solution (|a|=M) as a special case.

As we can see from the classification, the two-soliton solution
can produce a very rich family of analytic solutions with the
classical solutions of Schwarzschild and Kerr being special cases
of the general solution. Also the Manko et al. solution of Manko, Mielke & Sanabria-Gómez (2000) that has been
used previously by Berti & Stergioulas (2004); Stute & Camenzind (2002) to match the exterior
space-time of rotating neutron stars is a special case of the two-soliton solution
as we will see next.

Figure 1: The first two plots show the parameter space (a,b,k) of
the two soliton solution for a particular mass (M=1) from two
different view points so that the intricate foldings of the surface are better
understood. The 2D surface plotted is the constraint of
k which corresponds to the solution of Manko et al. used by
Berti & Stergioulas (2004). The third plot is the k=0 plane of the
parameter space which corresponds to all Kerr solutions.
The plots clearly show that the Manko et al. solution has no set of parameters to describe the case k=0
which corresponds to the Kerr and the Schwarzschild (for a=0) solutions, since there is no intersection of
the constraint of k and the Kerr plane in the appropriate range of parameters.
The two hyperbolas plotted on the plane k=0 are the only points were the constraint of k touches the plane
tangentially. As we can see these hyperbolas correspond to |a|>M, i.e., to hyper-extreme Kerr
space-times.

All types of solutions discussed above can be represented in a
three dimensional parameter space, the reduced-parameter space that was
mentioned in Case IIb. Although the two-soliton solution is
characterized by four parameters, one of them, the monopole mass
M, is simply a scaling parameter which can be used to reduce the
rest of the parameters to dimensionless ones. The three
dimensionless parameters thus formed,
(a/M,b/M,k/M2), are related to the
multipole moments (see eq. (2)) of the corresponding
space-time in the following way: The first parameter a/M
is the spin parameter (where a is the reduced angular momentum) which is the only
parameter, besides the mass, that uniquely characterizes a Kerr
space-time. The second parameter k/M2 expresses the
deviation of the quadrupole moment of the solution from the
quadrupole moment of the corresponding Kerr (the one with the same
a/M value); an increase of the value of k/M2
produces solutions that are less oblate than Kerr. The final
parameter b/M controls in a linear fashion the current
octupole moment. The actual deviation of the two soliton octupole
moment from the Kerr octupole moment depends on all three
parameters a/M, k/M2 and b/M.
Of course the higher moments are also affected by these parameters.

In this three dimensional parameter space, the plane
k/M2=0 corresponds to all types of Kerr solutions. This
is clear from the form of the Ernst potential along the axis,
where if one sets k=0 it reduces to the Ernst potential of the
Kerr solution,

e(z)=z−M−iaz+M−ia.

(23)

Obviously in this case the parameter b/M is
redundant; thus each line a/M=const, which is
parallel to the b/M axis on the plane k/M2=0,
corresponds to a single Kerr (modulo the mass of the black hole).

As mentioned in the introduction, the solution of Manko, Mielke & Sanabria-Gómez (2000)
has been used to describe the exterior of rotating neutron stars.
As it was briefly discussed above this solution is included in the two-soliton solution and
can be obtained by imposing a specific constraint on the
two-soliton parameters. The Manko et al. solution is obtained by
setting

k=−14[M2−(a−b)2]−M2b2M2−(a−b)2+ab.

(24)

This constraint defines a surface in the three parameter space
(a/M,b/M,k/M2) (see
Figure 1). The particular solution, depending on the
values of a,b, falls under either Case Ib or Case IIb, where either
κ− or κ+ is equal to zero, respectively. We should note
that the Manko et al. solution is the union of these two
cases. By substituting the above expression for k (eq. (24))
in the formula for the quadrupole moment (2), the
quadrupole moment takes the following value when a=0,

M2=−M4(M2+b2)2M2−b2.

(25)

This is why the quadrupole moment of the Manko et al. solution does
not vanish in the limit of zero rotation. From the above
expression one can see that the metric is not spherically
symmetric as one would expect for a non rotating object.
Especially for |b|<M the metric is oblate while for |b|>M the
metric is prolate.

Figure 2: Typical types of the surfaces that define the static
limit and the regions of CTCs for the various cases of the
two-soliton. On the left panels we have two typical figures that
correspond to Case Ia, in the middle panels the figures correspond to
Case IIa and on the right panels the figures correspond to Case III. The
solid curves correspond to the static limit and the dotted curves
correspond to the boundary of the regions with CTCs. For the Case
Ia, the region of CTCs is in contact with the axis of symmetry while
the same also applies to the surface of the static limit. For the Case
IIa the static limit is in contact with the axis, but the region
of CTCs has been detached. Finally for the Case III the inverse is true, that is
the region of CTCs is in contact with the axis of symmetry while the
static limit is detached from the axis. In all cases the upper plots correspond to slow rotation
while the lower plots correspond to faster rotation. As the rotation rate
increases the upper configurations are continuously deformed to the lower ones.

This anomalous behavior of the quadrupole moment is an important drawback for using the Manko et al. solution
to describe every rotating neutron star and it was pointed out by
Berti & Stergioulas (2004). In fact this analytic metric is good to match only rapidly rotating
neutron stars.

As shown by Manko, Mielke & Sanabria-Gómez (2000) this particular metric turns into a
Kerr metric if b2=a2−M2. Since all the parameters are
assumed real, then this corresponds to a hyperextreme Kerr metric (a≥M). In Figure 1 this is represented by the two
hyperbolas that lay outside the strip |a|<M on the Kerr plane (the plane
k=0) along which the intricate surface of the Manko et al. solution
tangentially touches the corresponding plane.

The two-soliton solution, which we will thoroughly study later on,
is a much better metric to describe the
exterior of an arbitrary rotating neutron star than the Manko et al. solution because
(i) the former has 4 independent parameters
(compared to the 3 independent parameters of the latter one) that
offer more flexibility to adjust the metric, and (ii) these 4
parameters are able to cover the whole space of the first 4
moments of the space-time, while the first 4 moments of the latter
metric are actually correlated with each other through the
dependence of k on the 3 independent parameters M,a,b,
that was mentioned previously.
Exactly this restriction renders the Manko et al. solution
inappropriate to describe slowly rotating neutron stars.

Before closing this section we will give a brief description
of the space-time characteristics of the two-soliton solution for the
range of parameter that we are going to use.

A horizon of a space-time is the boundary between the region where
stationary observers can exist and the region where such
observers cannot exist. For a stationary and axially symmetric
space-time, the stationary observers are those that have a
four-velocity that is a linear combination of the timelike and the
spacelike Killing vectors that the space-time possesses, i.e.

uμ=γ(ξμ+Ωημ)

(26)

where ξμ,ημ are the timelike and spacelike Killing fields
respectively and Ω is the observer’s angular velocity. The factor
γ is meant to normalize the four-velocity so that
gμνuμuν=−1. In order for the four-velocity to be
timelike, γ should satisfy the equation

γ−2=−gtt−2Ωgtϕ−Ω2gϕϕ,

(27)

and it should be γ−2>0, which corresponds to an Ω
taking values between the two roots

Ω±=−gtϕ±√(gtϕ)2−gttgϕϕgϕϕ.

(28)

This condition can not be satisfied when
(gtϕ)2−gttgϕϕ≤0. Thus the condition
(gtϕ)2−gttgϕϕ=0 defines the horizon. In the
case of the two-soliton, expressed in the Weyl-Papapetrou
coordinates, this condition corresponds to having ρ=0, since
ρ2=(gtϕ)2−gttgϕϕ. Thus the issue of
horizons is something that we will not have to face; in these
coordinates the whole space described corresponds to the exterior of any possible horizon.

Another issue is the existence of singularities.
Singularities might arise where the metric functions have
infinities. From equations (11-14)
one can see that singularities might exist where the functions
R±=√ρ2+(z±ξ+)2,r±=√ρ2+(z±ξ−)2 go to zero, or where the determinant E− goes to zero,
or where E+E∗−+E∗+E− goes to zero. Whether or
not some of these quantities vanish depends on which Case
the solution belongs to. A thorough investigation of the
singularities of the two-soliton is out of the scope of this
analysis; so we should only point out that for all the neutron star
models that we have studied and the corresponding parameters
of the two-soliton solution, any such singularities, when present, are always confined in the
region covered by the interior of the neutron star and thus
they do not pose any computational problems in our analysis.

The final issue is with regard to the existence of ergoregions
(i.e., regions where gμνξμξν=gtt<0) and regions with
closed timelike curves (CTCs) (i.e., regions where gabηaηb=gϕϕ<0). In
Figure 2 we have plotted the boundary surfaces of
such regions for the two-soliton metric. One can see that there are three distinct topologies
for these surfaces that are observed for the different two-soliton Cases.
In any case though, for all the models used here,
these surfaces are again confined at regions where the interior of the neutron star lays.

In the following sections we analyze the method that we are going to use to
obtain the right values for the parameters of the two-soliton
solution for each neutron star model and compare its properties
with the corresponding numerical metric. In order to do that, we
have constructed several sequences of numerical neutron star
models with the aid of the RNS numerical code of Stergioulas & Friedman (1995). The
numerical neutron star models used are the same models used by
Pappas & Apostolatos (2012) for demonstrating how to correct the numerical multipole
moments. They are produced using three equations of state (EOSs),
i.e., AU, FPS and L. The scope of using these models is twofold. First we are using them
to provide the appropriate parameters describing realistic
neutron stars in order to build the corresponding analytic
metrics. Then we use them as a testbed against which to
compare the analytic metrics and thus test their accuracy. As we
have already mentioned, we will use as matching conditions between the analytic
and the numerical metrics the first four non-zero
multipole moments. For the neutron star models that we have
studied, the corresponding analytic space-times
that are produced belong to three of the Cases of the
aforementioned classification, i.e., to Cases Ia, IIa and III.

Figure 3: Contour plots that point out what is the best choice for the
parameters of an analytic metric so that it matches well a numerical one. The
left plot shows the contour plots of the overall mismatch σij
between the analytic and the numerical metric for the tt (black curves) and
tϕ components (grey curves), respectively, as a function of the fractional deviation of
the quadrupole, δM2, and the current octupole, δS3, of the analytic metric
from those calculated directly from the
numerical metric (assuming the same mass and angular momentum though).
Since the contours of the σtt (almost horizontal) are
orthogonal to the ones of σtϕ (almost vertical), the
combination indicates an optimum choice for the multipole
moments of the analytic space-time. That choice is the moments that
have zero deviation from the moments of the numerical space-time.
The next two contour plots are similar contour plots but the contours
correspond to different quantities. The middle one is for the relative difference
between the analytic and numerical RISCO, while the
right one is for the overall difference between the analytic
and numerical orbital frequency Ω (defined analogously to the
overall mismatch of the metric components). Both these latter plots are
consistent with the first one. All plots correspond to the model
#15 for the AU EOS of the models presented by Pappas & Apostolatos (2012) and they give a representative
picture of what’s happening with all numerical models that we have computed.

When one attempts to match an analytic solution to a numerical one, it is
desirable to find a suitable matching criterion
that would be characteristic of the whole structure of the particular
numerical space-time, instead of just a finite region of it.
That is, the matching should be global and not local.
Berti & Stergioulas (2004) have argued that a suitable
global condition should be the matching of the first few multipole moments.
Indeed, the full set of multipole moments (as defined relativistically
by Geroch (1970); Hansen (1974); Fodor, Honselaers & Perjes (1989)) of a stationary and axially symmetric
space-time can fully specify the Ernst potential on the axis of symmetry. On the
other hand, when the Ernst potential along the axis of symmetry is given,
there is a space-time which is unambiguously specified by that Ernst
potential as it was shown by Xanthopoulos (1979, 1981); Hauser & Ernst (1981). Thus, the
full set of multipole moments are uniquely characterizing a space-time and they can be used as a
global matching condition.

When the space-time of a neutron star model is constructed from a numerical algorithm,
one can evaluate its mass moments M,Q,... and current moments J,S3,...
with an accuracy depending on the grid used to present the
numerical metric (for further discussion see Berti & Stergioulas (2004), and
Pappas & Apostolatos (2012)). Practically, the first few numerically evaluated moments can be
used as matching conditions to the analytic
space-time. The first four nonzero multipole moments of the two-soliton solution
as a function of its parameters M, a, b, k are shown in eq. (2)
from which it is clear that once we specify the mass and the angular
momentum of the space-time, the parameter k is uniquely
determined by the quadrupole moment Q≡M2 ,while the parameter b
is uniquely determined by the current octupole S3≡J3. Thus,
having constrained the four parameters of the two-soliton, we have
completely specified an analytic space-time that could be used to
describe the exterior of the particular neutron star model. What remains to be seen
is how well do the properties of the analytic space-time compare to those of the numerical one.

At this point there is an issue that should be addressed. Having
specified the first four non-zero moments of the two-soliton
metric, we have fixed all the higher moments of the space-time in a
specific manner related to the particular choice of the
analytic metric. These higher moments will probably
deviate from the ones of the numerical space-time. So the question
is, could one make a better choice when trying to match the analytic to
the numerical space-time than the one of setting the first four
analytic moments exactly equal to the first four numerical moments? To answer
this question we have performed the following test. For several
numerical models of uniformly rotating neutron stars that we have constructed, we
formed a set of two-soliton space-times for each neutron star model
that have the same mass M and angular momentum J with the numerical model, while
the quadrupoles and the current octupoles of each single two-soliton space-time
take the values M(a)2=M(n)2(1−δM2) and
S(a)3=S(n)3(1−δS3) respectively with various δM2 and
δS3 values. The quantities δM2 and δS3 denote the fractional differences of
the corresponding analytic moments of each two-soliton space-time from the numerical one. Then for
each one of these sets of moments we calculated the overall
mismatch between the analytic and the numerical metric functions,
which are defined (see Pappas & Apostolatos (2012)) as

Missing or unrecognized delimiter for \right

(29)

where RS is the radius r at the surface of the star,
and have thus constructed contour plots of σij on the plane
of δS3 and δM2, like the ones shown in Figure 3.
The same type of contour plots were drawn for other quantities as
well, like the relative difference of the RISCO and
the overall difference between the analytic and numerical orbital
frequency Ω for circular, equatorial orbits (defined in the
same fashion as it was defined the overall metric mismatch in eq. (29)).
The contour plots for the particular neutron star model that is shown in
Figure 3 are typical of the behavior that we observed in
all models. An important result is that the contours for the
overall mismatch σtt combined with the contours for
σtϕ give us the best choice for matching a numerical
to an analytic space-time; namely, the equation of the first four multipole
moments between the two space-times. That is because gtt seems to be
sensitive mainly in deviations from the numerical quadrupole and thus the
contours appear to be approximately horizontal and parallel to the axis of
δS3, while gtϕ seems to be sensitive mainly in deviations
from the numerical current octupole and thus the contours appear
to be approximately vertical and parallel to the axis of δM2
(the contours of σtϕ are almost orthogonal to the contours of σtt).

We should note here that the exact position of the optimum point in the contour plots of
σtt,σtϕ did not deviate from (0,0) by more than 3-4 per cent
in all cases studied. The largest deviations showed up in some of the fastest rotating models.

The conclusion is that what we expected to be true based on
theoretical considerations turns out to be exactly the case after the
implementation of the aforementioned test. Therefore in what follows, we
will set the first four non-zero multipole moments of the analytic space-time
equal to that of the corresponding numerical space-times.

4 Criteria for the comparison of the analytic to the
numerical space-time

Once we have constructed the analytic metric, appropriately matched to the
corresponding numerical one, we proceed to thoroughly compare the
two space-times. In order to do that, we should try again
to use criteria that are characteristic of the geometric structure
of the whole space-time, and if possible coordinate independent.
It would be preferable if these criteria are also
related to quantities that are relevant to astrophysical observations. Thus, if two space-times are in good
agreement, with respect to these criteria, they could be considered
more or less equivalent.On the other hand, such criteria, as well as possible
observations associated to them, could be used to distinguish
different space-times and consequently different compact
objects that are the sources of these space-times.

As a first criterion of comparing the metrics we will use the
direct comparison between the analytic and the numerical metric
components themselves. Although the metric components are quantities that are not
coordinate independent, they have specific physical meanings and can
be related to observable quantities. Thus, the gtt component
is related to the gravitational redshift of a photon and the
injection energy of a particle. The gtϕ component is
related to the frame dragging effect, and the angular velocity
ω=−gtϕ/gϕϕ of the zero angular momentum observers (ZAMOs).
Finally the gϕϕ component is related to the circumference of a circle at a particular radial
distance and defines the circumferential radius
Rcirc=C/2π=√gϕϕ. Also gϕϕ together with grr
are used to measure surface areas. So, if the relative difference
between the numerical and the analytic metric components

gtt=−f,gtϕ=fω,

gϕϕ=f−1ρ2−fω2,gρρ=gzz=f−1e2γ

(30)

is small, then one could consider the analytic metric as a good approximation of the numerical metric.

Another criterion for comparing an analytic to a numerical space-time is the location of the innermost stable
circular orbit (ISCO). Particles moving on the equatorial plane are governed by the equation of motion
(see for example Ryan (1995))

−gρρ(dρdτ)2=1−~E2gϕϕ+2~E~Lgtϕ+~L2gttρ2≡V(ρ),

(31)

where ~E and ~L are the conserved energy and
angular momentum parallel to the axis of symmetry, per unit mass.
V(ρ) is an effective potential for the radial motion and in
the case of orbits that are circular, we additionally have
the conditions dρ/dτ=0 and d2ρ/dτ2=0,
which are equivalent to the conditions for a local
extremum of the potential, i.e., V(ρ)=0, and dV(ρ)/dρ=0.
The radius of the ISCO is evaluated if we further demand the constraint
d2V(ρ)/dρ2=0, the physical meaning of which is
that the position of the circular orbit is also a turning point
of the potential. From these three conditions we can evaluate a
specific ρISCO and from that the
RISCO=√gϕϕ(ρISCO),
which we then compare to the corresponding numerical one. The position
of the ISCO is of obvious astrophysical interest since it is the
inner radius of an accretion disc and recently it
has been used to evaluate the rotation parameter of black holes from
fitting the continuous spectrum of the accretion disc
around them (see work by Shafee et al. (2006)).

Another criterion for comparing the metrics can be the orbital
frequency of circular equatorial orbits Ω. The orbital
frequency is given by the equation

Ω(ρ)=−gtϕ,ρ+√(gtϕ,ρ)2−gtt,ρgϕϕ,ρgϕϕ,ρ.

(32)

Apart from the orbital frequency one could also use the precession
frequencies of the almost circular and almost equatorial
orbits, i.e., the precession of the periastron Ωρ and
the precession of the orbital plane Ωz. These frequencies
are derived from the perturbation of the equation of motion

−gρρ(dρdτ)2−gzz(dzdτ)2=V(ρ,z)

(33)

around the circular equatorial orbits. In this
expression, V(ρ,z) is the same effective potential which was defined in
the second equation of (31) the z dependence of which now
has not been omitted as in eq. (31). The perturbation frequencies derived from the above
equation are then given with respect to the metric functions as

κ2a

=

−gaa2{(gtt+gtϕΩ)2(gϕϕρ2),aa

(34)

−

2(gtt+gtϕΩ)(gtϕ+gϕϕΩ)(gtϕρ2),aa

+

(gtϕ+gϕϕΩ)2(gttρ2),aa},

where the index a takes either the value ρ or z to
obtain the frequency of the radial or the vertical perturbation
respectively. These expressions are evaluated on the equatorial
plane (at z=0); thus they are functions of ρ alone. The corresponding
precession frequencies are given by the difference between the orbital frequency and the perturbation
frequency:

Ωa=Ω−κa.

(35)

These quantities are quite interesting with respect to
astrophysical phenomena as well. More specifically, they can be
associated to the orbital motion of material accreting onto a
compact object through an accretion disc. These very frequencies have
been proposed to be connected to the observed quasi-periodic
modulation (QPOs) of the X-ray flux of accretion discs that are present in
X-ray binaries (see Stella (2000); van der Klis (2006); Boutloukos et al. (2006); Lamb (2007)).

Finally, the last criterion that we will use to compare metrics is
the quantity Δ~E of circular orbits, which expresses
the energy difference of the orbits per logarithmic orbital
frequency interval as one moves from one circular orbit to the
next towards the central object. This quantity is defined as

Δ~E=−Ωd~EdΩ,

(36)

where the energy per unit mass ~E is given by the
expression

~E=−gtt−gtϕΩ√−gtt−2gtϕΩ−gϕϕΩ2.

(37)

The quantity Δ~E is a measure of the energy that a
particle has to lose in order to move from one circular orbit to
another closer to the central object so that the frequency increases by one e-fold. The quantity,
Δ~E, is associated to the emission of gravitational
radiation and was used by Ryan (1995) to measure the multipole
moments of the space-time from gravitational waves emitted by test particles orbiting
in that background. The same quantity can also be associated to accretion discs and in
particular, in the case of thin discs, it would correspond to the amount of
energy that the disc will radiate as a function of the radius
from the central object and thus it will be related to the
temperature profile of the disc and consequently to the total
luminosity of the disc (for a review on accretion discs see Krolik (1999)).

The last set of criteria, i.e., the frequencies
and Δ~E, are related to specific observable
properties of astrophysical systems, in particular of accretion discs around
compact objects; thus they are very useful and
relevant to astrophysics (for an application see the work by
Pappas (2012)).

Figure 4: These are the plots of the logarithm of the relative
difference of the analytic to the numerical metric components
gtt, gtϕ, √gϕϕ=Rcirc, and gzz
on the equatorial plane, as well as gtt on the axis of
symmetry. The plots shown here are indicative of all comparisons
performed for all the neutron star models that we have constructed. These
particular plots are drawn for the model #10 of the AU EOS, the characteristics of which are
presented in Table II of the supplement of Pappas & Apostolatos (2012). Different curves correspond to
different metrics, i.e., the two-soliton (solid curve), the Manko et al. solution
with the negative root for b (dashed-dotted curve) and the one with
the positive root (dashed-double-dotted), and finally the Hartle-Thorne (dotted) metric. In most cases shown here, the Manko et al. curves are on top of the two-soliton curve.
The case of the gtϕ component, where the two-soliton
curve is almost an order of magnitude below the Manko et al. one
for a large interval of radii, is a notable exception.
The Hartle-Thorne’s characteristic failure to describe gtϕ is also evident.

The results of the comparison between the analytic and the
numerical metrics, describing the exterior of realistic neutron
stars, that are presented here is indicative of all comparisons
performed for every numerically constructed neutron star model.
The models that we have used as a testbed of comparison are briefly presented in the
Appendix and are the same models used in Pappas & Apostolatos (2012).

For illustrative reasons
we have also plotted comparisons between the numerical and the Manko et al. solution (Manko, Mielke & Sanabria-Gómez (2000))
which was used by Berti & Stergioulas (2004), as well as
comparisons between the numerical and the Hartle-Thorne metric
(Hartle & Thorne (1968)). The reason for using these two metrics is that on the
one hand the Hartle-Thorne metric is considered to be a good
approximation of slowly rotating relativistic stars and on the
other hand the Manko et al. metric has been shown by Berti & Stergioulas (2004) to be a good
approximation for relativistic stars with fast
rotation. In cases with slow rotation rates, for which corresponding models of the
Manko et al. metric cannot be constructed, we have used only the
Hartle-Thorne metric to compare, though. In the case of models with fast rotation
besides the Manko et al. metric, we have used the
Hartle-Thorne metric as well. for these cases we treated the Hartle-Thorne metric as a three parameter exterior
metric, where the three parameters are the mass M, the angular
momentum J and the reduced quadrupole q=M2/M3 of the neutron star. It should be noted that this is not
a consistent way to use the Hartle-Thorne metric, since the
quadrupole and the angular momentum in the Hartle-Thorne cannot
take arbitrary values, while the metric is essentially a two
parameter solution (parameterized by the central density of the
corresponding slowly rotating star and a small parameter
ε that corresponds to the fraction of the angular velocity
of the star relative to the Keplerian angular velocity of the
surface of the star) that has to be properly matched to an
interior solution following the procedure described by
Berti et al. (2004). Here we are taking some leeway in using
Hartle-Thorne for fast rotation since it is used simply for
illustrative purposes and not in order to draw any conclusions from it.

In Figure 4 we present the comparison of the
various analytic metric functions (using the two-soliton, the
Manko et al., and the Hartle-Thorne solutions) to the corresponding
numerical ones for a single model constructed using the EOS AU
(model #10 of the AU EOS the characteristics of which are
presented in Table II of the supplement of Pappas & Apostolatos (2012)). The
figures display the relative difference between the various analytic and
the numerical metric functions gtt,gtϕ,√gϕϕ and gzz=gρρ
on the equatorial plane, as well as the function gtt on the axis of
symmetry.

The general picture we get from these figures is typical for all models
constructed using all three equations of state, that is AU, FPS and
L. The overall comparison of the two-soliton to the numerical
metrics shows that this analytic metric is an excellent substitute
of the numerical space-time both for slow and fast rotating models,
with an accuracy that is everywhere outside the neutron star
always better than about 1/1000 for all the metric functions (there is an exception to that for the
comparison of the gtt metric component right at the pole where for
some models it’s fractional difference is a bit smaller than
1/100). In comparison to the other two analytic metrics discussed above, we
see that for the models that a Manko et al. metric can be found,
that is for the rapidly rotating neutron stars, both this metric and the two-soliton metric perform
very well (actually the Manko et al. solution performs better than
Berti & Stergioulas (2004) had initially found as was shown by
Pappas & Apostolatos (2012)) and there are only tiny
differences between the two-soliton and the Manko et al. analytic
metrics. The greatest difference between the two-soliton and the
Manko et al. appears to be in the gtϕ component of the
metric, where the two-soliton is usually clearly better. This was anticipated
because the gtϕ component of the metric is, as we have
shown in section 3, more sensitive to the value of S3, which can be
suitably adjusted in the two-soliton metric, but not in the Manko et al. solution. For
the rapidly rotating models the Hartle-Thorne is not such a good
representation of the numerical metric as the other two metrics. That also was
expected since the Hartle-Thorne metric is not suitable for
fast rotation. Hartle-Thorne’s failure is more evident in the
gtϕ component of the metric which is consistent with
Hartle-Thorne’s vanishing spin octupole S3 (in the Appendix we show
that S3=0 for the Hartle-Thorne metric).

For the slowly rotating models, there are no Manko et al.
solutions to compare to the numerical metric, so in these cases
the only alternative is the Hartle-Thorne metric. We should say
again that the consistent way to calculate the Hartle-Thorne
parameters is the one described by Berti et al. (2004), but as it is
discussed by Pappas & Apostolatos (2012) the parameters of the
Hartle-Thorne metric (specifically the parameter q which is the
reduced quadrupole) consistently calculated are in very good
agrement with the numerical multipole moments, so these moments
were used straightforwardly for the construction of the
Hartle-Thorne metric.
Again we saw that the two-soliton performs better when compared
to the numerical space-time than the Hartle-Thorne
metric. We should note though that the problem of
Hartle-Thorne’s metric to acurately describe the gtϕ metric component is present
even at slow rotation.

Having demonstated the overall superiority of two-soliton,
compared to Manko et al. and Hartle-Thorne, to accurately
describe the metric functions of any numerical neutron star model,
in the following comparisons we will only compare the two-soliton
quantities to the corresponding numerical ones.

Figure 5: These figures show the relative difference between the
numerical and the analytic RISCO for the
neutron star models calculated using the AU EOS. The relative
difference is shown only for the models that have their ISCO
outside the surface of the neutron star. The negative values of
j correspond to counter-rotating orbits (retrograde). The
data for the co-rotating orbits can be found in Table 2.

The next quantity we have used for comparison is the position of the ISCO. In
Figure 5 we present the relative difference between
the numerical and the analytic ISCO for all neutron star
models constructed with the AU EOS for both prograde and
retrograde orbits (the latter are indicated by negative parameter
j≡J/M2). The general conclusion is that for all models
constructed using all three EOSs the ISCO of the analytic metric
does not deviate by more than 4 per cent from the ISCO of the
corresponding numerical model and such deviations are observed for
the prograde orbits of the fastest rotating models. We should note
that we don’t perform any comparison with numerical models
the radius of the surface of which on the equatorial plane is further than
the radius of the corresponding ISCO; therefore the points that would
correspond to these models are missing from the plots. At this
point we should mention that apart from the position of the
marginally stable circular orbit for the particles, there is also
the position of the unstable photon circular orbit that could also be
used as a criterion for comparison. This orbit though
is usually, for the prograde case, below the surface of the
neutron star (while for the retrograde it is usually outside
the star) and it doesn’t have an immediately measurable effect3.

We continue to the results of the comparison of the various
frequencies associated to the circular orbits on the equatorial
plane, i.e., Ω,Ωρ,κρ,Ωz, and
κz. The analytic orbital frequency compares very well to
the numerical orbital frequency for all the models with the
relative difference being in all cases smaller than ∼10−3.
This result is very important, since the orbital frequency
together with the RISCO are relevant to
observations from accretions discs. Typical plots of Ω and
the relative difference between the analytic and the numerical one
as functions of the logarithm of the distance from the central
object are shown in Figure 6.

Figure 6: The plots show the various frequencies
and the relative difference of the analytic to the corresponding
numerical ones. The top row of plots shows on the left the numerical (dotted)
and the analytic (solid) orbital frequency Ω for comparison reasons (they are hardly distinguishable),
while on the right it is plotted the logarithm of their relative difference. The
middle row shows on the left the numerical (dotted) and analytic
(solid) precession frequencies Ωρ, in the middle the
corresponding radial oscillations frequencies κρ
(same correspondence of lines) and on the right the logarithm of the
relative difference of the latter oscillation frequencies. Finally the bottom row
shows the corresponding plots for the plane precession (left) and the
vertical oscillation frequencies (middle), Ωz, and κz respectively, as
well as the logarithm of the relative difference of the latter frequencies (right). All
plots are for the model #10 of the AU EOS presented by Pappas & Apostolatos (2012).

With regard to the comparison of the radial and vertical
perturbation frequencies κρ, and κz,
respectively, and the corresponding precession frequencies
Ωρ, and Ωz, things are more complicated. These
frequencies include second derivatives of the metric functions in
their calculation. Consequently the results of the corresponding numerical
calculations are plagued by accuracy problems. In particular
numerical calculation of the second derivatives induces
artificial oscillations in the results. These issues have been
discussed earlier by Berti et al. (2004). Following the suggestions of
N. Stergioulas, we found two ways to mitigate the problem. The first
one is to calculate the frequencies directly in the coordinates
that the RNS code produces the metric functions so as to avoid any
numerical errors caused by the transformations of the coordinates
and the metric functions themselves. Then one would only have to
identify the coordinates of the points at which the frequencies
are calculated with the corresponding Weyl-Papapetrou coordinates
(the calculated frequencies themselves are the frequencies that
static observers at infinity measure so they don’t depend on the
coordinate system used). The second one is to smooth out these artificial oscillations by taking a three point average of the
frequencies. The efficiency of this technique has been tested in the
case of the non-rotating models (the exterior of which are described by a
Schwarzschild metric) and has been verified to give trustworthy
results. Another thing that we should also consider is that in the
case of κz, the values are very close to the corresponding values of
Ω; consequently there is low accuracy in the
calculation of the precession frequency Ωz=Ω−κz
in some cases. That is why we consider as a better indicator of
the actual ability of the analytic metric to capture the behavior
of the numerical metric, the deviations in the oscillation frequencies
κρ and κz, instead of Ωρ and Ωz.
Nevertheless, we present both
the precession and the oscillation frequencies of the analytic
metric in comparison to those of the numerical metric, together with
the relative difference between the analytic and the numerical
oscillation frequencies, in order to get a clearer picture of the
comparison. All these plots, again for the model #10 of the AU
EOS, the characteristics of which are
presented in Table II of the supplement of
Pappas & Apostolatos (2012), are shown in Figure 6.

Generally, the relative differences in κρ between the
numerical and the analytic metrics is small and in some cases they
could climb up to 10 per cent near the ISCO. This is due to the
fact that the radial oscillation frequency tends to zero as the
ISCO is approached, causing an increase of the relative difference.
In contrast, for κz, the relative differences are always below 1 per cent
at the ISCO. Now the picture is inverted for the precession
frequencies, which is related to the fact that in this case the
small quantities are the Ωz’s. The overall picture we obtain
is that the analytic frequencies capture quite well the behavior
of the numerical frequencies both qualitatively and
quantitatively. This is especially evident in the case of the
κz frequency where in some cases it becomes greater than
the orbital frequency, an effect more prominent in the models of
EOS L, which the two-soliton metric captures quite well. The
importance of this effect and its relevance to QPOs is further
discussed by Pappas (2012).

The final comparison criterion is the quantity Δ~E.
The numerical computation of this quantity has similar difficulties
with the precession frequencies; these issues could be fixed by performing the
same tricks to avoid numerical oscillations.
In Figure 7 we show for the same model of a
rotating neutron star as in the previous cases, the quantity Δ~E computed from the
numerical and the analytic metric on the left, and their relative
difference on the right. Again we see that the two-soliton metric
describes with high accuracy the Δ~E which we obtain
numerically from the numerical models.
We remind here that this quantity is relevant for the emitted
spectrum of a thin accretion disc and its temperature profile, as
well as for the efficiency of the disc, i.e., the amount of
kinetic energy transformed to radiation.

We close this section with Table 2, where we
present for all numerical models of EOS AU (the multipole
moments of which are given in the Tables of Pappas & Apostolatos (2012)) the
parameters and the type of the two-soliton metric along with a
few quantities of astrophysical interest, and specifically their
comparison between the ones calculated using the analytic and the
numerical space-times. These quantities are, the circumferential
radius at the ISCO RISCO, the efficiency
η=1−~EISCO of a thin accretion disc (if there was one
around the particular neutron star), the orbital frequency
at the position of the ISCO ΩISCO (this is a
frequency expected to show up in QPOs if the latter are related to
the orbital motion), and finally the vertical oscillation frequency at the
ISCO (κz)ISCO (this could also be related to QPOs). The Table shows that the relative
differences between the numerical and the analytic quantities is of the order
of 1 per cent or lower.

Figure 7: Plots showing the analytic (solid) and
numerical (dotted) Δ~E (left) as well as the logarithm of their
relative difference (right) for the model #10 of AU EOS presented in
Pappas & Apostolatos (2012).

model

Case

j

b

k

ΔRISCO

Δη

ΔΩISCO

(Δκz)ISCO

(km)

(km2)

(per cent)

(per cent)

(per cent)

(per cent)

1

III

0.

-3.

0

0.003

0.032

0.006

-

2

Ia

0.2015

-0.0784

-0.5271

0.013

0.019

0.003

0.024

3

IIa

0.3126

-0.1305

-1.2503

0.283

0.142

0.027

0.141

4

IIa

0.414

-0.1858

-2.1664

-

-

-

-

5

IIa

0.4749

-0.224

-2.8362

-

-

-

-

6

IIa

0.5297

-0.2626

-3.517

-

-

-

-

7

IIa

0.5789

-0.3014

-4.1972

-

-

-

-

8

IIa

0.617

-0.3351

-4.7738

-

-

-

-

9

IIa

0.651

-0.368

-5.3277

-

-

-

-

10

IIa

0.6618

-0.3793

-5.5114

-

-

-

-

11

III

0.

-3.

0

0.003

-0.021

-0.004

-

12

III

0.194524

-0.0915

-0.1528

-0.021

0.022

0.002

-0.034

13

III

0.309849

-0.1669

-0.4288

-0.018

0.025

0.003

-0.071

14

III

0.406932

-0.2417

-0.8163

0.019

0.03

0.005

0.113

15

III

0.485572

-0.3158

-1.2698

0.148

0.114

0.019

0.04

16

III

0.550214

-0.3897

-1.7658

0.391

0.204

0.035

0.217

17

III

0.603381

-0.4628

-2.279

0.778

0.463

0.077

0.477

18

III

0.645447

-0.5317

-2.7705

1.274

0.83

0.137

0.815

19

III

0.676639

-0.5916

-3.197

1.798

1.205

0.197

1.15

20

III

0.706299

-0.6585

-3.6602

2.453

1.831

0.297

1.638

21

III

0.510282

-0.3726

-0.7958

0.043

0.026

0.005

0.15

22

III

0.510617

-0.3717

-0.8179

0.054

0.046

0.008

0.058

23

III

0.514032

-0.3731

-0.8721

0.063

0.05

0.009

0.031

24

III

0.520506

-0.3789

-0.9409

0.083

0.016

0.004

0.065

25

III

0.547452

-0.4083

-1.1827

0.164

0.072

0.013

0.153

26

III

0.587439

-0.4607

-1.5504

0.346

0.206

0.034

0.164

27

III

0.626593

-0.5214

-1.9574

0.63

0.391

0.064

0.362

28

III

0.659098

-0.5806

-2.3502

0.991

0.645

0.104

0.64

29

III

0.694585

-0.6577

-2.8456

1.551

1.018

0.162

0.948

30

III

0.713165

-0.7054

-3.1406

1.95

1.349

0.214

1.221

Table 2: For each neutron star model, that was constructed by the RNS code using the AU EOS, we have computed a
number of parameters that are related to the particular analytic two-soliton metric which better approximates
the numerical metric. These parameters are (i) the type (Case) of the
two-soliton space-time that each particular model corresponds to, (ii) the parameters k and b of the
two-soliton as well as the spin parameter j and, (iii) the relative
difference between the following analytic and numerical quantities: RISCO, η,
ΩISCO and (κz)ISCO, for the models that the ISCO lies outside the surface
of the star. The models shown here are the same models presented in Pappas & Apostolatos (2012). The rest of the
physical parameters of these models, such as the mass and all the other multipole moments of
the models, can be found there. All the relative differences are given as a percentage.

In this work we have tested whether the four-parameter two-soliton analytic metric,
which was derived by Manko, Martin & Ruiz (1995), can be used as a trustworthy approximation for the space-time around
all kind of neutron stars.
To match the particular analytic metric to a specific neutron star model,
which was produced through the numerical code RNS of Stergioulas & Friedman (1995), we have
used as matching conditions the first four non-zero multipole
moments. Our choice was justified, apart from theoretical reasoning, by comparing the numerical metrics with
different analytic metrics produced by slightly varying two of their moments (quadrupole and spin-octupole).
The comparison showed clearly that the best matching comes from imposing the
condition that the parameters of the analytic metric should be such
that the analytic space-time acquires the first four non-vanishing moments of the numerical metric.

Having demonstrated the appropriateness of the matching conditions, we proceeded to compare the various numerical
neutron star space-times with the corresponding analytic space-times. To perform the comparison, we have assumed several
criteria having in mind that they should correspond to geometric and physical properties of the space-time,
with a special interest in physical quantities that could be associated to astrophysical
processes that are usually observed from the vicinity of neutron stars.

The result of these comparisons is that the two-soliton space-time
can reproduce the properties of the space-time around realistic
neutron stars, and in particular it can reproduce all
astrophysically interesting properties. Probably the most
important fact is that the analytic metric can capture properties
of the neutron star space-time that a corresponding Kerr space-time
could not, such as the behavior of the precession frequencies of
almost circular and almost equatorial orbits. A typical
example is shown in Figure 8, where we present the analytic and
numerical frequencies of the precession of the orbital plane for a
model constructed using the L EOS. The possible importance and
implication of this, i.e,. the capability of the two-soliton to
capture this particular behavior in contrast to the Kerr geometry, was
further discussed in Pappas (2012).

Generally the two-soliton metric can be a very useful tool for studying
phenomena that happen around all kind of neutron stars and are quite sensitive
to more realistic and accurate geometries than the ones used
so far. Relying on a single analytic metric for all neutron stars is
practically more favorable than using numerical space-times, or more than one analytic
metrics depending on the type of the neutron star. Thus,
the two-soliton metric can be further used for more elaborate applications
such as the ones explored by Psaltis (2008); Bauböck, Psaltis, Özel & Johannsen (2012); Bauböck, Psaltis & Özel (2012).

Figure 8: Plot showing the numerical (dotted) and analytic (solid)
precession frequencies Ωz for the model #10 of the L
EOS. One can see how the frequencies that are calculated from the analytic
metric capture the qualitative behavior of the corresponding numerical
frequencies. The parameters (i.e., the multipole moments) for this
model can be found in Pappas & Apostolatos (2012).

Acknowledgments

We would like to thank Kostas Kokkotas and Kostas Glampedakis
for many useful discussions, and Nikos Stergioulas for
providing us access to his RNS numerical code as well as solid advice.
This work has been supported by the I.K.Y. (IKYDA 2010).
G.P. would also like to acknowledge DAAD scholarship number A/12/71258.

Appendix A Neutron star models

In order to construct the analytic space-time exterior to a compact
object, one has to choose the appropriate multipole moments. For
neutron stars, these moments could be computed from numerical models that are
constructed with realistic equations of state. There are several
schemes developed for numerically integrating stellar models (see
Stergioulas & Friedman (1995), and for an extended list of numerical schemes see
Stergioulas (2003)). We have used Stergioulas’s RNS code for the
construction of the models.

In order to cover more space on the “neutron
star parameter space”, we constructed numerical
neutron star models with equations of state of varying
stiffness. For that purpose we have chosen AU as a typical soft
EOS, FPS as a representative moderate stiff EOS, and L to describe
stiff EOS. Having the numerical models ready, we proceeded in
evaluating their multipole moments according to the algorithm described in
Pappas & Apostolatos (2012). The parameters then used to construct the analytic space-time models,
i.e., M,a,k, and b, are evaluated from the first four
multipole moments (M,J,M2,S3) of each model by inverting equations
(2).

For the specifics of the various models chosen here, we have followed
Berti & Stergioulas (2004). We have constructed the same constant rest-mass
sequences as the ones presented in Berti & Stergioulas (2004) for the
corresponding equations of state. For every equation of state, 3
sequences of 10 models were constructed, that corresponded to:

a sequence corresponding to a neutron star of 1.4M⊙ in the non-rotating limit,

a sequence terminating at the maximum-mass neutron star in
the non-rotating limit,

a supermassive sequence that does not terminate at a
non-rotating model at its lower rotation limit.

All the sequences end at the mass-shedding limit on the side of
fast rotation, i.e., at the limit were the angular velocity of a
particle at the equator is equal to the Keplerian velocity at that
radius. These sequences are the so called evolutionary
sequences.

All the parameters for the computed models can be found in
Pappas & Apostolatos (2012).

Appendix B The functions of the two-soliton

In this Appendix we present the full expressions for writing the
metric functions of the two-soliton. The determinants that are
given in Sec. 2 and appear in the formulas for the metric
functions, can be substituted with the following expressions,
starting with E±=A∓B, where the functions A,B are
given as:

A

=

−16dk(r−r++R−R+)M2

(38)

−

[A−1(R−r++r−R+)−A−2(R−r+−r−R+)]κ2+

+

[A+1(r−R−+r+R+)−A+2(r−R−−r+R+)]κ2−,

A±1

=

(4(a−b)2(ab±d−k)−4((a−2b)b±d)M2),

(39)

A±2

=

4i((a−b)(ab±d−k)−bM2)κ±,

(40)

B

=

2κ−κ+M[(R+−R−)B−1+(r−−r+)B+1

(41)

+

(R−+R+)B−2−(r−+r+)B+2],

B±1

=

i(2k(a−b)+b(M2−a2+b2±κ−κ+))

(42)

×

(κ+±κ−),

B±2

=

2d(M2−a2+b2±κ−κ+).

(43)

The determinant H can be substituted as H=−L, where L is
given by the expression,