Realistic exact solution for the exterior field of a rotating neutron star
Leonardo A. Pachón, Jorge A. Rueda,  José D. Sanabria-Gómez, PHYSICAL REVIEW D 73, 104038 (2006)


A new six-parametric, axisymmetric, and asymptotically flat exact solution of Einstein-Maxwell field
equations having reflection symmetry is presented. It has arbitrary physical parameters of mass, angular
momentum, mass-quadrupole moment, current octupole moment, electric charge, and magnetic dipole, so
it can represent the exterior field of a rotating, deformed, magnetized, and charged object; some properties
of the closed-form analytic solution such as its multipolar structure, electromagnetic fields, and
singularities are also presented. In the vacuum case, this analytic solution is matched to some numerical
interior solutions representing neutron stars, calculated by Berti and Stergioulas [E. Berti and N.
Stergioulas, Mon. Not. R. Astron. Soc. 350, 1416 (2004)], imposing that the multipole moments be the
same. As an independent test of accuracy of the solution to describe exterior fields of neutron stars, we
present an extensive comparison of the radii of innermost stable circular orbits (ISCOs) obtained from the
Berti and Stergioulas numerical solutions, the Kerr solution [R. P. Kerr, Phys. Rev. Lett. 11, 237 (1963)],
the Hartle and Thorne solution [J. B. Hartle and K. S. Thorne, Astrophys. J. 153, 807 (1968)], an analytic
series expansion derived by Shibata and Sasaki [M. Shibata and M. Sasaki, Phys. Rev. D 58, 104011
(1998)], and our exact solution.We found that radii of ISCOs from our solution fits better than others with
realistic numerical interior solutions.

Off[General::spell]
conj[g_] := 
        g /. Complex[u_, v_] :> Complex[u, -v];
real[g_] := (g + conj[g])/2 // Together;
imag[g_] := (g - conj[g])/(2 I) // Together;

Axis data

Ms = 1.47663;

m = 1.405*Ms;    a = 0.8121;
k = -1.001;          s = 0.727;
q = 1 ;                      b = 
      10;

Ernst potential on the symmetry axis

EZ = ((I s - k z + (-I a - m) (z^2) + (z^3))/(I s - k z + (-I a + m) (z^2) + (z^3))) // Simplify;
FZ = ((I b z + q (z^2))/(I s - k z + (-I a + m) (z^2) + (z^3))) // Simplify;

Multipolar Moments

See Physical Review D 73 104038 (2006) for a physical disscussion of each term

Mass multipoles

({{M0, m}, {M1, 0}, {M2, -((a^2) - k) m}, {M3, 0}})

Current multipoles

({{J0, 0}, {J1, a m}, {J2, 0}, {J3, -m ((a^3) - 2 a k + s)}})

Electric multipoles

({{Q0, q}, {Q1, 0}, {Q2, -a b - (a^2) q + k q}, {Q3, 0}})

Magnetic multipoles

({{B0, 0}, {B1, b + a q}, {B2, 0}, {B3, -(a^2) b + b k - (a^3) q + 
            2 a k q - q s}})

Auxiliary Definitions

EZC = conj[EZ]; FZC = conj[FZ];
alphas[i_] := Root[EZ + EZC + 2 FZ FZC, i] // N // 
          Chop;

α1 = alphas[1];
α2 = alphas[2]; α3 = 
      alphas[3];
α4 = alphas[4]; α5 = alphas[5]; α6 = 
      alphas[6];

r1 = 
        Sqrt[ρ^2 + (z - α1)^2]; 
r2 = Sqrt[ρ^2 + (z - α2)^2];  
    r3 = Sqrt[ρ^2 + (z - α3)^2];
r4 = Sqrt[ρ^2 + (z - α4)^2]; 
    r5 = Sqrt[ρ^2 + (z - α5)^2]; 
    r6 = Sqrt[ρ^2 + (z - α6)^2];

g1 = r1 + α1 - z;
g2 = r2 + α2 - z; 
    g3 = r3 + α3 - z;
g4 = r4 + α4 - z; 
    g5 = r5 + α5 - z;
g6 = r6 + α6 - z;

beta[
          i_] := Root[I s - 
                k z + (-I a + m) (z^2) + (z^3), i] // 
            N // Chop;
β1 = beta[1]; β2 = beta[2]; β3 = beta[3];

e[1]=-((2 β(1^2) m)/((β1-β2) (β1-β3)))//Simplify;
e[2]=((2 ((β2)^2) m)/((β1-β2) (β2-β3)))//Simplify;
e[3]=-((2 ((β3)^2) m)/((β1-β3) (β2-β3)))//Simplify;

cβ1 = conj[β1]; cβ2 = conj[β2]; 
    cβ3 = conj[β3];
ec[1] = conj[e[1]]; ec[2] = conj[e[2]];  
    ec[3] = conj[e[3]];

f1=((I b β1+β1^2*q)/((β1-β2) (β1-β3)))//Simplify;
f2=((-I b β2-((β2)^2) *q)/((β1-β2) (β2-β3)))//Simplify;
f3=((I b β3+((β3)^2) q)/((β1-β3) (β2-β3)))//Simplify;
fc1 = 
        conj[f1] // Simplify;
fc2 = 
        conj[f2] // Simplify;
fc3 = 
        conj[f3] // Simplify;

f[1] = 
        f1/(α1 - β1)  +  f2/(α1 - β2) + 
          f3/(α1 - β3);
f[2] = 
        f1/(α2 - β1)  +  f2/(α2 - β2) + 
          f3/(α2 - β3);
f[3] = 
        f1/(α3 - β1)  +  f2/(α3 - β2) + 
          f3/(α3 - β3);
f[4] = 
        f1/(α4 - β1)  +  f2/(α4 - β2) + 
          f3/(α4 - β3);
f[5] = 
        f1/(α5 - β1)  +  f2/(α5 - β2) + 
          f3/(α5 - β3);
f[6] = 
        f1/(α6 - β1)  +  f2/(α6 - β2) + 
          f3/(α6 - β3);

h1[i_] := ec[1] + 2  fc1 f[i];
h2[i_] := ec[2] + 2  fc2 f[i];
h3[i_] := ec[3] + 2  fc3 f[i];

Ep=(({{1, 1, 1, 1, 1, 1, 1}, {1, (r1/(α1 - β1)), (r2/(α2 - β1)), (r3/(α3 - β1)), (r4/(α4 - β1)), (r5/(α5 - β1)), (r6/(α6 - β1))}, {1, (r1/(α1 - β2)), (r2/(α2 - β2)), (r3/(α3 - β2)), (r4/(α4 - β2)), (r5/(α5 - β2)), (r6/(α6 - β2))}, {1, (r1/(α1 - β3)), (r2/(α2 - β3)), (r3/(α3 - β3)), (r4/(α4 - β3)), (r5/(α5 - β3)), (r6/(α6 - β3))}, {0, h1[1]/(α1 - cβ1), h1[
                        2]/(α2 - cβ1), h1[
                        3]/(α3 - cβ1), h1[
                        4]/(α4 - cβ1), h1[
                        5]/(α5 - cβ1), h1[
                        6]/(α6 - cβ1)}, {0, h2[1]/(α1 - cβ2), h2[
                        2]/(α2 - cβ2), h2[
                        3]/(α3 - cβ2), h2[
                        4]/(α4 - cβ2), h2[
                        5]/(α5 - cβ2), h2[
                        6]/(α6 - cβ2)}, {0, h3[1]/(α1 - cβ3), h3[
                        2]/(α2 - cβ3), h3[
                        3]/(α3 - cβ3), h3[
                        4]/(α4 - cβ3), h3[
                        5]/(α5 - cβ3), h3[
                        6]/(α6 - cβ3)}}))//Simplify;
Em=(({{1, 1, 1, 1, 1, 1, 1}, {-1, (r1/(α1 - β1)), (r2/(α2 - β1)), (r3/(α3 - β1)), (r4/(α4 - β1)), (r5/(α5 - β1)), (r6/(α6 - β1))}, {-1, (r1/(α1 - β2)), (r2/(α2 - β2)), (r3/(α3 - β2)), (r4/(α4 - β2)), (r5/(α5 - β2)), (r6/(α6 - β2))}, {-1, (r1/(α1 - β3)), (r2/(α2 - β3)), (r3/(α3 - β3)), (r4/(α4 - β3)), (r5/(α5 - β3)), (r6/(α6 - β3))}, {0, h1[1]/(α1 - cβ1), h1[
                        2]/(α2 - cβ1), h1[
                        3]/(α3 - cβ1), h1[
                        4]/(α4 - cβ1), h1[
                        5]/(α5 - cβ1), h1[
                        6]/(α6 - cβ1)}, {0, h2[1]/(α1 - cβ2), h2[
                        2]/(α2 - cβ2), h2[
                        3]/(α3 - cβ2), h2[
                        4]/(α4 - cβ2), h2[
                        5]/(α5 - cβ2), h2[
                        6]/(α6 - cβ2)}, {0, h3[1]/(α1 - cβ3), h3[
                        2]/(α2 - cβ3), h3[
                        3]/(α3 - cβ3), h3[
                        4]/(α4 - cβ3), h3[
                        5]/(α5 - cβ3), h3[
                        6]/(α6 - cβ3)}}))//Simplify;
G=(({{0, g1, g2, g3, g4, g5, g6}, {-1, (r1/(α1 - β1)), (r2/(α2 - β1)), (r3/(α3 - β1)), (r4/(α4 - β1)), (r5/(α5 - β1)), (r6/(α6 - β1))}, {-1, (r1/(α1 - β2)), (r2/(α2 - β2)), (r3/(α3 - β2)), (r4/(α4 - β2)), (r5/(α5 - β2)), (r6/(α6 - β2))}, {-1, (r1/(α1 - β3)), (r2/(α2 - β3)), (r3/(α3 - β3)), (r4/(α4 - β3)), (r5/(α5 - β3)), (r6/(α6 - β3))}, {0, h1[1]/(α1 - cβ1), h1[
                        2]/(α2 - cβ1), h1[
                        3]/(α3 - cβ1), h1[
                        4]/(α4 - cβ1), h1[
                        5]/(α5 - cβ1), h1[
                        6]/(α6 - cβ1)}, {0, h2[1]/(α1 - cβ2), h2[
                        2]/(α2 - cβ2), h2[
                        3]/(α3 - cβ2), h2[
                        4]/(α4 - cβ2), h2[
                        5]/(α5 - cβ2), h2[
                        6]/(α6 - cβ2)}, {0, h3[1]/(α1 - cβ3), h3[
                        2]/(α2 - cβ3), h3[
                        3]/(α3 - cβ3), h3[
                        4]/(α4 - cβ3), h3[
                        5]/(α5 - cβ3), h3[
                        6]/(α6 - cβ3)}}))//Simplify;
H=(({{z, 1, 1, 1, 1, 1, 1}, {-β1, (r1/(α1 - β1)), (r2/(α2 - β1)), (r3/(α3 - β1)), (r4/(α4 - β1)), (r5/(α5 - β1)), (r6/(α6 - β1))}, {-β2, (r1/(α1 - β2)), (r2/(α2 - β2)), (r3/(α3 - β2)), (r4/(α4 - β2)), (r5/(α5 - β2)), (r6/(α6 - β2))}, {-β3, (r1/(α1 - β3)), (r2/(α2 - β3)), (r3/(α3 - β3)), (r4/(α4 - β3)), (r5/(α5 - β3)), (r6/(α6 - β3))}, {ec[
                    1], h1[1]/(α1 - cβ1), h1[
                        2]/(α2 - cβ1), h1[
                        3]/(α3 - cβ1), h1[
                        4]/(α4 - cβ1), h1[
                        5]/(α5 - cβ1), h1[
                        6]/(α6 - cβ1)}, {ec[
                    2], h2[1]/(α1 - cβ2), h2[
                        2]/(α2 - cβ2), h2[
                        3]/(α3 - cβ2), h2[
                        4]/(α4 - cβ2), h2[
                        5]/(α5 - cβ2), h2[
                        6]/(α6 - cβ2)}, {ec[
                    3], h3[1]/(α1 - cβ3), h3[
                        2]/(α2 - cβ3), h3[
                        3]/(α3 - cβ3), h3[
                        4]/(α4 - cβ3), h3[
                        5]/(α5 - cβ3), h3[
                        6]/(α6 - cβ3)}}))//Simplify;
K=(({{(1/(α1 - β1)), (1/(α2 - β1)), (1/(α3 - β1)), (1/(α4 - β1)), (1/(α5 - β1)), (1/(α6 - β1))}, {(1/(α1 - β2)), (1/(α2 - β2)), (1/(α3 - β2)), (1/(α4 - β2)), (1/(α5 - β2)), (1/(α6 - β2))}, {(1/(α1 - β3)), (1/(α2 - β3)), (1/(α3 - β3)), (1/(α4 - β3)), (1/(α5 - β3)), (1/(α6 - β3))}, {h1[1]/(α1 - cβ1), h1[2]/(α2 - 
                        cβ1), h1[3]/(α3 - 
                        cβ1), h1[4]/(α4 - 
                        cβ1), h1[5]/(α5 - 
                        cβ1), h1[6]/(α6 - cβ1)}, {h2[1]/(α1 - cβ2), h2[2]/(α2 - 
                        cβ2), h2[3]/(α3 - 
                        cβ2), h2[4]/(α4 - 
                        cβ2), h2[5]/(α5 - 
                        cβ2), h2[6]/(α6 - cβ2)}, {h3[1]/(α1 - cβ3), h3[2]/(α2 - 
                        cβ3), h3[3]/(α3 - 
                        cβ3), h3[4]/(α4 - 
                        cβ3), h3[5]/(α5 - 
                        cβ3), h3[6]/(α6 - cβ3)}}))//Simplify;
F=(({{0, f[1], f[2], f[3], f[4], f[5], f[
                    6]}, {-1, (r1/(α1 - β1)), (r2/(α2 - β1)), (r3/(α3 - β1)), (r4/(α4 - β1)), (r5/(α5 - β1)), (r6/(α6 - β1))}, {-1, (r1/(α1 - β2)), (r2/(α2 - β2)), (r3/(α3 - β2)), (r4/(α4 - β2)), (r5/(α5 - β2)), (r6/(α6 - β2))}, {-1, (r1/(α1 - β3)), (r2/(α2 - β3)), (r3/(α3 - β3)), (r4/(α4 - β3)), (r5/(α5 - β3)), (r6/(α6 - β3))}, {0, h1[1]/(α1 - cβ1), h1[
                        2]/(α2 - cβ1), h1[
                        3]/(α3 - cβ1), h1[
                        4]/(α4 - cβ1), h1[
                        5]/(α5 - cβ1), h1[
                        6]/(α6 - cβ1)}, {0, h2[1]/(α1 - cβ2), h2[
                        2]/(α2 - cβ2), h2[
                        3]/(α3 - cβ2), h2[
                        4]/(α4 - cβ2), h2[
                        5]/(α5 - cβ2), h2[
                        6]/(α6 - cβ2)}, {0, h3[1]/(α1 - cβ3), h3[
                        2]/(α2 - cβ3), h3[
                        3]/(α3 - cβ3), h3[
                        4]/(α4 - cβ3), h3[
                        5]/(α5 - cβ3), h3[
                        6]/(α6 - cβ3)}}))//Simplify;
Q =(({{f1 + f2 + f3, 0, f[1], f[2], f[3], f[4], f[5], f[
                    6]}, {z, 1, 1, 1, 1, 1, 1, 1}, {-β1, -1, (r1/(α1 - β1)), (r2/(α2 - β1)), (r3/(α3 - β1)), (r4/(α4 - β1)), (r5/(α5 - β1)), (r6/(α6 - β1))}, {-β2, -1, (r1/(α1 - β2)), (r2/(α2 - β2)), (r3/(α3 - β2)), (r4/(α4 - β2)), (r5/(α5 - β2)), (r6/(α6 - β2))}, {-β3, -1, (r1/(α1 - β3)), (r2/(α2 - β3)), (r3/(α3 - β3)), (r4/(α4 - β3)), (r5/(α5 - β3)), (r6/(α6 - β3))}, {ec[1], 0, h1[1]/(α1 - cβ1), h1[
                        2]/(α2 - cβ1), h1[
                        3]/(α3 - cβ1), h1[
                        4]/(α4 - cβ1), h1[
                        5]/(α5 - cβ1), h1[
                        6]/(α6 - cβ1)}, {ec[2], 0, h2[1]/(α1 - cβ2), h2[
                        2]/(α2 - cβ2), h2[
                        3]/(α3 - cβ2), h2[
                        4]/(α4 - cβ2), h2[
                        5]/(α5 - cβ2), h2[
                        6]/(α6 - cβ2)}, {ec[3], 0, h3[1]/(α1 - cβ3), h3[
                        2]/(α2 - cβ3), h3[
                        3]/(α3 - cβ3), h3[
                        4]/(α4 - cβ3), h3[
                        5]/(α5 - cβ3), h3[
                        6]/(α6 - cβ3)}}))//Simplify;

EpD = Det[Ep] // Simplify;
EmD = Det[Em] // Simplify;
GD = Det[G] // Simplify;
FD = Det[F] // Simplify;
HD = Det[H] // Simplify;
KD = Det[K] // Simplify;
ID = Det[Q] // Simplify;

EpC = conj[EpD] // Simplify;
EmC = conj[EmD] // Simplify;
GC = 
        conj[GD] // 
          Simplify;               
HC = conj[HD] // Simplify;
KC = 
        conj[KD] // 
          Simplify;               
FC = conj[FD] // Simplify; IC = conj[ID] // Simplify;

Ernst Potentials on all the planex plane

\(\(\( (*\ 
      Ersnt\ potential\ Φ\ *) \)\(
\)\(\(PPhi = 
        FD\/EmD // 
          Simplify;\)\(\ \ \ \ \ \ \ \ \ \)\(
\)\( (*\ 
      Ersnt\ potential\ \cal {Ε}\ *) \)\(
\)
    \(PE = EpD\/EmD // Simplify;\)\(\ \ \ \ \ \ \ \ \ \ \)\)\)\)

\(\(\( (*\ 
      Imaginary\ Part\ of\ \cal {Ε}\ Potentail\ *) \)\(\ \)\(
\)\(Omega = imag[PE];\)\(\ \ \)\)\)

phi = 
        real[PPhi];

(* 
      Imaginary Part of Φ Potentail *)  
Psi = imag[PPhi] // Simplify;

Kinnersley's Potential

(* 
      Kinnersley potential *) 
Phi2 = -I ID/EmD // 
          Simplify;       
 (* 
      Magnetic potential *) 

    PotMag = real[Phi2];

Metric Functions

f = (EpD*EmC + EpC*EmD + 2  FD FC)/(2*EmD*EmC);
e2γ = (EpD*EmC + EpC*EmD + 2  FD FC)/(2*KD*KC*r1*r2*r3*
              r4*r5*r6);
ω = 
        2*imag[EmD*HC - EmC*GD - FD IC]/(EpD*EmC + EpC*EmD + 
                2  FD FC);

Field Equations

f Integral Condition

left = f;

right = real[PE] +  PPhi conj[PPhi];

side1 = 
      left /. {ρ -> 2.1,  z -> 3.5} // Chop

side2 = 
      right /. {ρ -> 2.1,  z -> 3.5} // Chop

(side1 - side2) // Chop

ω Integral Conditions

Condition I

left = D[ω, ρ];

right = -(ρ/(f^2))*(D[Omega, z]  +  
              2 (phi D[Psi, z]  -  Psi D[phi, z]));

side1 = 
      left /. {ρ -> 2.1,  z -> 3.5} // Chop

side2 = 
      right /. {ρ -> 2.1,  z -> 3.5} // Chop

(side1 - side2) // Chop

Condition II

left = D[ω, z];

right = (ρ/(f^2))*(D[Omega, ρ]  +  
              2 (phi D[Psi, ρ]  -  
                    Psi D[phi, ρ]));

side1 = 
      left /. {ρ -> 2.1,  z -> 3.5} // Chop

side2 = 
      right /. {ρ -> 2.1,  z -> 3.5} // Chop

(side1 - side2) // Chop

γ Integral Conditions

Condition I

left = D[e2γ, ρ] /(2*e2γ);

right = 
        0.25*(ρ/(f^2))*(D[PE, ρ] + 
                2  conj[PPhi]  D[PPhi, ρ] )*(D[
                  conj[PE], ρ]  + 2 PPhi D[conj[PPhi], ρ]) - 
          0.25*(ρ/(f^2))*(D[PE, z]  + 
                2 conj[PPhi]  D[PPhi, z])*(D[conj[PE], z]  + 
                2 PPhi D[conj[PPhi], z]) - (ρ/
                f)*(D[PPhi, ρ]  D[conj[PPhi], ρ]  - 
                D[PPhi, z]  D[conj[PPhi], z]);

side1 = 
      left /. {ρ -> 2.1,  z -> 3.5} // Chop

side2 = 
      right /. {ρ -> 2.1,  z -> 3.5} // Chop

(side1 - side2) // Chop

Condition II

left = D[e2γ, z] /(2*e2γ);

right = 
        0.25 (ρ/(f^2))*((D[PE, ρ]  + 
                      2 conj[PPhi]  D[PPhi, ρ] )*(D[conj[PE], 
                        z]  + 2 PPhi D[ conj[PPhi], z]) + 
                conj[(D[PE, ρ]  + 
                        2 conj[PPhi]  D[PPhi, ρ] )*(D[conj[PE], 
                          z]  + 
                        2 PPhi D[ conj[PPhi], z])]) - (ρ/
                f)*( 
              D[ conj[PPhi], ρ]  D[PPhi, z] + 
                conj[ D[ conj[PPhi], ρ]  D[PPhi, 
                      z]]);

side1 = 
      left /. {ρ -> 2.1,  z -> 3.5} // Chop

side2 = 
      right /. {ρ -> 2.1,  z -> 3.5} // Chop

(side1 - side2) // Chop