An analytic spacetime exterior to neutron stars

# An all-purpose metric for the exterior of any kind of rotating neutron star

## 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 , the orbital frequency of circular geodesics, and their epicyclic frequencies , as well as the change of the energy of circular orbits per logarithmic change of orbital frequency . 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:
gravitation – stars: neutron – equation of state – methods: analytical – relativity – accretion.
12

## 1 Introduction

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, . 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.

## 2 The two-soliton solution

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 and are functions of the Weyl-Papapetrou coordinates (). By introducing the complex potential , 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 is the metric function , 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 are polynomials of of order 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 are the roots of the polynomial and are complex coefficients appropriately chosen so that the latter form of (eq. 4) is equal to the former one (eq. 3). Subsequently one determines the roots of equation

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

where denotes complex conjugation. These roots are denoted as , with and from these one defines the complex functions . All these functions and roots are then used as building blocks for the following determinants

 E±=∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣11⋯1±1⋮±1C0⋮0∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣ (6)
 G=∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣0R1+ξ1−z⋯R2n+ξ2n−z−1⋮−1C0⋮0∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣ (7)
 H=∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣z1⋯1−β1⋮−βnCe∗1⋮e∗n∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣ (8)
 K0=∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣1ξ1−β1⋯1ξ2n−β1⋮⋱⋮1ξ1−βn⋯1ξ2n−βne∗1ξ1−β∗1⋯e∗1ξ2n−β∗1⋮⋱⋮e∗nξ1−β∗n⋯e∗nξ2n−β∗n∣∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣∣ (9)

where is the matrix

 C=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝R1ξ1−β1⋯R2nξ2n−β1⋮⋱⋮R1ξ1−βn⋯R2nξ2n−βne∗1ξ1−β∗1⋯e∗1ξ2n−β∗1⋮⋱⋮e∗nξ1−β∗n⋯e∗nξ2n−β∗n⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (10)

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 and their complex conjugates that appear in the determinants cancel out (the 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 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 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 corresponds to the mass monopole of the space-time, the parameter is the reduced angular momentum, is the deviation of the reduced quadrupole from the corresponding Kerr quadrupole (the one that has the same and ), and is associated with the deviation of the current octupole moment from the current octupole of the corresponding Kerr.

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 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.

1. This case is characterized by two real roots . The Kerr family of solutions, which corresponds to is definitely not included in this family of solutions.

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

2. This is a degenerate case where the roots coincide. The degeneracy is due to being zero which corresponds to the parameter constraint . 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 , the previous constraint corresponds to a two-dimensional surface.

2. In this case, the roots are either complex or imaginary, since . Furthermore this means that there are nonvanishing values of where the functions assume zero value, which then leads to singularities at the corresponding points.

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

2. 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 (). Similarly to Case Ib this solution is also represented by a two dimensional surface in the reduced-parameter space.

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

3. 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 plane. The Kerr and the Schwarzschild solutions lie entirely within this family of solutions.

4. 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 (), while subcase IVb has a quadruple root () and the computation of the metric functions needs special treatment. We should also note that cases IVb and IVc include the extreme Kerr solution () 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.

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 , 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, , are related to the multipole moments (see eq. (2)) of the corresponding space-time in the following way: The first parameter is the spin parameter (where is the reduced angular momentum) which is the only parameter, besides the mass, that uniquely characterizes a Kerr space-time. The second parameter expresses the deviation of the quadrupole moment of the solution from the quadrupole moment of the corresponding Kerr (the one with the same value); an increase of the value of produces solutions that are less oblate than Kerr. The final parameter 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 , and . Of course the higher moments are also affected by these parameters.

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

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

Obviously in this case the parameter is redundant; thus each line , which is parallel to the axis on the plane , 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 (see Figure 1). The particular solution, depending on the values of , 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 (eq. (24)) in the formula for the quadrupole moment (2), the quadrupole moment takes the following value when ,

 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 the metric is oblate while for the metric is prolate.

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 . Since all the parameters are assumed real, then this corresponds to a hyperextreme Kerr metric (). In Figure 1 this is represented by the two hyperbolas that lay outside the strip on the Kerr plane (the plane ) 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 on the 3 independent parameters , 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 . In order for the four-velocity to be timelike, should satisfy the equation

 γ−2=−gtt−2Ωgtϕ−Ω2gϕϕ, (27)

and it should be , which corresponds to an taking values between the two roots

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

This condition can not be satisfied when . Thus the condition defines the horizon. In the case of the two-soliton, expressed in the Weyl-Papapetrou coordinates, this condition corresponds to having , since . 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 go to zero, or where the determinant goes to zero, or where 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 ) and regions with closed timelike curves (CTCs) (i.e., regions where ). 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.

## 3 Matching the analytic to the numerical solution

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 and current moments 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 , , , 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 is uniquely determined by the quadrupole moment ,while the parameter is uniquely determined by the current octupole . 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 and angular momentum with the numerical model, while the quadrupoles and the current octupoles of each single two-soliton space-time take the values and respectively with various and values. The quantities and 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 is the radius at the surface of the star, and have thus constructed contour plots of on the plane of and , 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 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 combined with the contours for 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 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 , while 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 (the contours of are almost orthogonal to the contours of ).

We should note here that the exact position of the optimum point in the contour plots of did not deviate from 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 component is related to the gravitational redshift of a photon and the injection energy of a particle. The component is related to the frame dragging effect, and the angular velocity of the zero angular momentum observers (ZAMOs). Finally the component is related to the circumference of a circle at a particular radial distance and defines the circumferential radius . Also together with 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 and are the conserved energy and angular momentum parallel to the axis of symmetry, per unit mass. is an effective potential for the radial motion and in the case of orbits that are circular, we additionally have the conditions and , which are equivalent to the conditions for a local extremum of the potential, i.e., , and . The radius of the ISCO is evaluated if we further demand the constraint , 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 and from that the , 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 . 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, is the same effective potential which was defined in the second equation of (31) the 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 takes either the value or to obtain the frequency of the radial or the vertical perturbation respectively. These expressions are evaluated on the equatorial plane (at ); 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 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 is given by the expression

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

The quantity 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 -fold. The quantity, , 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 , 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)).

## 5 Results of the comparison

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 , the angular momentum and the reduced quadrupole 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 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 and on the equatorial plane, as well as the function 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 for all the metric functions (there is an exception to that for the comparison of the metric component right at the pole where for some models it’s fractional difference is a bit smaller than ). 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 component of the metric, where the two-soliton is usually clearly better. This was anticipated because the component of the metric is, as we have shown in section 3, more sensitive to the value of , 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 component of the metric which is consistent with Hartle-Thorne’s vanishing spin octupole (in the Appendix we show that 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 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 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.

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 ). 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., , and . 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 . This result is very important, since the orbital frequency together with the 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.

With regard to the comparison of the radial and vertical perturbation frequencies , and , respectively, and the corresponding precession frequencies , and , 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 , the values are very close to the corresponding values of ; consequently there is low accuracy in the calculation of the precession frequency 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 , instead of and . 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 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 , 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 ’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 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 . 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 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 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 , the efficiency of a thin accretion disc (if there was one around the particular neutron star), the orbital frequency at the position of the 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 (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.

## 6 Conclusions

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).

## 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., , and , are evaluated from the first four multipole moments () 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:

1. a sequence corresponding to a neutron star of in the non-rotating limit,

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

3. 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 , where the functions 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 can be substituted as , where is given by the expression,

 L = (r−r++R−R+)L1+(r−+R−+r++R+)L2 (44) − (r−−R−+r+−R+)L3+(r−R−−R+r+)L−5 + (r−−r++R−−R+)L−4+(r−R+−R−r+)L+5 + (r−−r+−R−+R+)L+4+(R−r++r−R+)L+6 + (R−r−+r+R+)L−6,
 L1 = −16dk(ia−ib+M−z)M2, (45) L2 = i4dM(a+iM)[M4−2(a2−b2−2k)M2 (46) + (a−b)2((a+b)2−4k)], L3 = 4dM[ia3−Ma2−i(b2+M2+2k)a+M3 (47) + 2ibk+b2M]κ−κ+,
 L±4 = 4κ∓M{ibM5+abM4 (48) − i(2ba2−ka−b(b2±d+2k))M3 − (2ba3−ka2−b(b2±d+3k)a + k(2b2∓d+k))M2 + (2k2−(a2+3ba∓2d)k+b(a+b)(a2∓d)) × i(a−b)