ANALYTIC EXTERIOR SOLUTION APPROXIMATING 
                RAPIDLY ROTATING RELATIVISTIC STARS
                 
                BASED ON MANKO ET AL. PRD, 61, 081501(R)


                      AUTHOR:   Nikolaos Stergioulas
                      
                      ADDRESS: Department of Physics	
                      		     Aristotle University of Thessaloniki
                      		     54124 Thessaloniki, Greece
                      		     niksterg@astro.auth.gr
                      		    
                      PUBLICATION:   Berti, E. & Stergioulas, N., "Approximate Matching of Analytic and 
                      			   Numerical Solutions for Rapidly Rotating Neutron Stars ", gr-qc/0310061

Off[General::spell]
$DefaultFont={"Times",10};
plotoptions={Frame->True, Axes->None,AspectRatio->1};

DEFINE METRIC AND MULTIPLE MOMENTS

Define D for metric

Dme = (4*(kd^2*x^2 - δ*y^2)^2 + 
                2*kd*M*x*(2*
                        kd^2*(x^2 - 1) + (2*δ + a*b - 
                            b^2)*(1 - y^2)) + (a - 
                      b)*((a - b)*(d - δ) - 
                      M^2*b)*(y^4 - 1) - 4*d^2)^2 + 
          4*y^2*(2*kd^2*(x^2 - 1)*(kd*x*(a - b) - M*b) - 
                  2*M*b*δ*(1 - 
                        y^2) + ((a - b)*(d - δ) - 
                        M^2*b)*(2*kd*x + M)*(1 - 
                        y^2))^2;

Define E for metric

Eme = (4*(kd^2*(x^2 - 1) +  δ*(1 - 
                              y^2))^2 + (a - 
                      b)*((a - b)*(d - δ) - 
                      M^2*b)*(1 - y^2)^2)^2 - 
          16*kd^2*(x^2 - 1)*(1 - 
                y^2)*((a - b)*(kd^2*(x^2 - y^2) + 
                        2*δ*y^2) + M^2*b*y^2)^2;

Define F for metric

\(\(Fme = 
        8*kd^2*\((x^2 - 
                1)\)*\((\((a - b)\)*\((kd^2*\((x^2 - y^2)\) + 
                      2*δ*y^2)\) + 
                y^2*M^2*b)\)*\((
\ \ \ \ \ \ \ \ kd*M*
                  x*\((\((2*kd*x + M)\)^2 - 
                      2*y^2*\((2*δ + a*b - b^2)\) - a^2 + b^2)\) - 
                2*
\ \ \ \ \ \ \ \ \ \ \ \ y^2*\((4*δ*d - 
                      M^2*b^2)\))\) + \((4*\((kd^2*\((x^2 - 
                              1)\) + δ*\((1 - 
                              y^2)\))\)^2
			 + \((a - 
                      b)\)*\((\((a - b)\)*\((d - δ)\) - 
                      M^2*b)\)*\((1 - y^2)\)^2)\)*\((4*\((2*kd*M*b*x + 
                      2*M^2*b)\)*\((kd^2*\((x^2 - 1)\) + δ*\((1 - 
                            y^2)\))\)
				 + \((1 - 
                      y^2)\)*\((\((a - b)\)*\((M^2*b^2 - 
                            4*δ*d)\)
							 - \((4*kd*M*x + 
                            2*M^2)\)*\((\((a - b)\)*\((d - δ)\) - 
                            M^2*b)\))\))\);\)\)

Define functions f, e^2\gamma and w

f = Eme/Dme;
e2gamma = Eme/(16*kd^8*(x^2 - y^2)^4);
w = (y^2 - 1)*Fme/Eme;

Define metric components

gtt = -f;
gtphi = f*w;
grhorho = e2gamma/f;
gphiphi = (rho^2/f) - f*w^2;

Define x and y coordinates

x = (1/2/kd)*(rp + rm);
y = (1/2/kd)*(rp - rm);

Define r+ and r-

rp = Sqrt[rho^2 + (z + kd)^2];
rm = Sqrt[rho^2 + (z - kd)^2];

Define kd, δ and d

kd = Sqrt[d + δ];
d = (M^2 - (a - b)^2)/4;
δ = -(M^2*
                b^2)/(M^2 - (a - b)^2);

Define quadrupole moment Q

Q = -M*(a^2 - a*b + d - δ);

Define l=3  current moment

S3 = -M*(a^3 - 2*a^2*b + a*(b^2 + 2*d - 2*δ) + 
              b*(δ - d));

COMPUTE NUMERICAL EXAMPLE

Numerical examples: EOS FPS (choose one of the following models)

Slowly rotating model (no real analytic solution):

MNum = 1.8209*1.47536;
JNum = 2.5193;
QNum = -4.303;
aNum = N[JNum/MNum, 16];
S3Num = -6.549;
rpRNS = 10.487 + 2.961
rmRNS = 10.487 + 8.900

Rapidly rotating model:

MNum = 1.8637*1.47536;
JNum = 4.9508;
QNum = -19.795;
aNum = N[JNum/MNum, 16];
S3Num = -63.673;
rpRNS = 13.818 + 0.0
rmRNS = 13.818 + 9.492

Numerical value for quadrupole moment

Qb = N[Q  /.  {a -> aNum,  M -> MNum},  
        16]

Plot (analytical quadrupole - numerical quadrupole)

\(\(Plot[{0, Qb - \((QNum)\)}, {b, \(-6\), 6}, {Frame -> True, 
          Axes -> None, AspectRatio -> 1, 
          FrameLabel -> {"b", \*""Q-\!\(Q\_N\)""}}];\)\)

Find numerical value of b, for which Qb=numerical Q

bm = 
      FindRoot[Qb - (QNum) == 0,  {b,  {-0.557,  0.0}}]
bp = FindRoot[Qb - (QNum) == 0,  {b,  {0.3,  3}}]

Evaluate S3 for the two branches

S3bm  =  S3  /.  a -> aNum  /. M  ->  MNum /.  
        bm
S3bp  =  S3  /.  a -> aNum  /. M  ->  MNum /.  
        bp

Relative error (%) in S3 for the negative branch

100*(S3bm - S3Num)/S3Num

COMPUTATION OF ISCO

wz0  = 
        Simplify[w  /.  z -> 0 /.  M -> MNum  /.  
              a  -> aNum   /. bm];
fz0 = 
        Simplify[f  /.  z -> 0  /.  M -> MNum  /.  
              a  -> aNum   /. bm];
wr  =  
        Simplify[D[wz0, rho]  /.  M -> MNum  /.  
              a  -> aNum   /. bm];
wrr = 
        Simplify[D[wr, rho] /.  M -> MNum  /.  
              a  -> aNum   /. bm];
fr  =  
        Simplify[D[fz0, rho] /.  M -> MNum  /.  
              a  -> aNum   /. bm];
frr = 
        Simplify[D[fr, rho] /.  M -> MNum  /.  
              a  -> aNum   /. bm];
sqisco = (wr^2*fz0^4 + fr*rho*(2*fz0 - fr*rho))^(1/
              2);
iscopbm  =  ( 
            wr*wrr*fz0^5*rho*(2*fz0 - fr*rho)  +  
              wr^2*fz0^4*(2*fz0^2 + (-fr^2 + frr*fz0)*rho^2) + 
              wr*fz0^2*
                sqisco*(2*fz0^2 + 2*fr^2*rho^2 - 
                    fz0*rho*(4*fr + frr*rho))  +  
              rho*(2*fz0 - fr*rho)*(3*fr*fz0^2 - 4*fr^2*fz0*rho + 
                    fr^3*rho^2 + 
                    fz0^2*(frr*rho - wrr*fz0*sqisco)))/ (fz0^2*
              rho^2*(wr^2*fz0^4 + 3*fr*fz0*rho - fr^2*rho^2 - 
                  fz0^2*(2 + 
                        wr*sqisco)));
rhoiscopbm = 
        FindRoot[iscopbm, {rho, 8, 15}];
rplus = 
      Sqrt[gphiphi  /.  z -> 0 /.  M -> MNum  /.  
              a  -> aNum   /. bm  /. rhoiscopbm]

wz0  = 
        Simplify[w  /.  z -> 0 /.  M -> MNum  /.  
              a  -> -aNum   /. 
            bm];
fz0 = 
        Simplify[f  /.  z -> 0  /.  M -> MNum  /.  
              a  -> -aNum   /. 
            bm];
wr  =  
        Simplify[D[wz0, rho]  /.  M -> MNum  /.  
              a  -> -aNum   /. 
            bm];
wrr = 
        Simplify[D[wr, rho] /.  M -> MNum  /.  
              a  -> -aNum   /. 
            bm];
fr  =  
        Simplify[D[fz0, rho] /.  M -> MNum  /.  
              a  -> -aNum   /. 
            bm];
frr = 
        Simplify[D[fr, rho] /.  M -> MNum  /.  
              a  -> -aNum   /. 
            bm];
sqisco = (wr^2*fz0^4 + fr*rho*(2*fz0 - fr*rho))^(1/
              2);
iscombm  =  ( 
            wr*wrr*fz0^5*rho*(2*fz0 - fr*rho)  +  
              wr^2*fz0^4*(2*fz0^2 + (-fr^2 + frr*fz0)*rho^2) + 
              wr*fz0^2*
                sqisco*(2*fz0^2 + 2*fr^2*rho^2 - 
                    fz0*rho*(4*fr + frr*rho))  +  
              rho*(2*fz0 - fr*rho)*(3*fr*fz0^2 - 4*fr^2*fz0*rho + 
                    fr^3*rho^2 + 
                    fz0^2*(frr*rho - wrr*fz0*sqisco)))/ (fz0^2*
              rho^2*(wr^2*fz0^4 + 3*fr*fz0*rho - fr^2*rho^2 - 
                  fz0^2*(2 + 
                        wr*sqisco)));
rhoiscombm = 
        FindRoot[iscombm, {rho, 10, 30}];
rminus = 
      Sqrt[gphiphi  /.  z -> 0 /.  M -> MNum  /.  
              a  -> -aNum   /. bm  /. 
          rhoiscombm]

Relative errors (%)

100*(rplus - rpRNS)/rpRNS
100*(rminus - rmRNS)/rmRNS

COMPUTATION OF ISCO WITH THE SHIBATA-SASAKI FORMULA TO VARIOUS ORDERS

jp = aNum/MNum;
jm = -aNum/MNum;
Q2 = -QNum/MNum^3;
q3p = -S3Num/MNum^4;
q3m = S3Num/MNum^4;
Q4 = Q2^4;

    
rpSS4 = 6*
            MNum*(1 - 0.54433*j - 0.22619*j^2 + 0.17989*Q2 - 0.23002*j^3 + 
                0.26296*j*Q2 - 0.05317*q3 - 0.29693*j^4 + 0.44546*j^2*Q2 - 
                0.06249*Q2^2  + 0.01544*Q4 - 0.11310*j*q3 )  /. 
          j -> 
            jp  /. q3->q3p 
rpSS4noQ4 = 6*
            MNum*(1 - 0.54433*j - 0.22619*j^2 + 0.17989*Q2 - 0.23002*j^3 + 
                0.26296*j*Q2 - 0.05317*q3 - 0.29693*j^4 + 0.44546*j^2*Q2 - 
                0.06249*Q2^2  + -0.11310*j*q3 )  /. 
          j -> jp  /. q3 -> q3p
rpSS3 = 6*
            MNum*(1 - 0.54433*j - 0.22619*j^2 + 0.17989*Q2 - 0.23002*j^3 + 
                0.26296*j*Q2 - 0.05317*q3)  /. j -> jp  /. 
        q3 -> q3p
rmSS4 = 6*
            MNum*(1 - 0.54433*j - 0.22619*j^2 + 0.17989*Q2 - 0.23002*j^3 + 
                0.26296*j*Q2 - 0.05317*q3 - 0.29693*j^4 + 0.44546*j^2*Q2 - 
                0.06249*Q2^2  + 0.01544*Q4 - 0.11310*j*q3 )  /. 
          j -> jm  /. q3 -> q3m
rmSS4noQ4 = 6*
            MNum*(1 - 0.54433*j - 0.22619*j^2 + 0.17989*Q2 - 0.23002*j^3 + 
                0.26296*j*Q2 - 0.05317*q3 - 0.29693*j^4 + 0.44546*j^2*Q2 - 
                0.06249*Q2^2  - 0.11310*j*q3 )  /. 
          j -> 
            jm  /. q3->q3m 
rmSS3 = 6*
            MNum*(1 - 0.54433*j - 0.22619*j^2 + 0.17989*Q2 - 0.23002*j^3 + 
                0.26296*j*Q2 - 0.05317*q3 )  /. 
          j -> 
            jm  /. q3->q3m

Relative errors (%)

100*(rpSS4 - rpRNS)/rpRNS
100*(rpSS4noQ4 - rpRNS)/rpRNS
100*(rpSS3 - rpRNS)/rpRNS

100*(rmSS4 - rmRNS)/rmRNS
100*(rmSS4noQ4 - rmRNS)/rmRNS
100*(rmSS3 - rmRNS)/rmRNS

COMPARE WITH THE MANKO ET AL. ANSATZ

bmam = -Sqrt[
            a^2 + 2*a*M*Δ - 
              M^2];
bmap = 
        Sqrt[a^2 + 2*a*M*Δ - M^2];
Qma = -(1 + (1/
                      2)*Δ^2 - (a*Δ^2*(a*Δ^2 - M*Δ + bma))/(2*(M^2 - a^2 + 
                          a^2*Δ^2 - 
                          2*a*M*Δ)))*a^2*
          M;
Qmam = Qma   /. bma  -> bmam  /. a -> aNum  /. 
          M -> MNum;
Qmap = Qma   /. bma  -> bmap  /. a -> aNum  /. 
          M -> MNum;

Plot[
      Qmam - QNum, {Δ, 0.46, 0.52}]
Plot[Qmap - QNum, {Δ, 0.8969, 1.2}]

Deltam = 
      FindRoot[Qmam - QNum, {Δ, 0.48, 
          0.54}]
Deltap = FindRoot[Qmap - QNum, {Δ, 1, 10.2}]

bmamN = bmam  /. a -> aNum  /. M -> MNum  /.  
        Deltam
bmapN = bmap  /. a -> aNum  /. M -> MNum  /. 
        Deltap

KERR METRIC

Metric components

bKerr = Sqrt[a^2 - M^2];
gttK  =  gtt  /. b -> bKerr  /. a -> aNum  /.  
          M -> MNum;
gtphiK  =  gtphi  /. b ->   bKerr /. 
            a -> aNum  /.  M -> MNum;
grhorhoK  =  grhorho  /. b ->   bKerr /. 
            a -> aNum  /.  M -> MNum;
gphiphiK  =  gphiphi  /. b ->   bKerr /. 
            a -> aNum  /.  M -> MNum;

Q and S3

QKerr = Simplify[Q  /. b -> bKerr]
QKerrN = QKerr  /. a -> aNum  /. 
        M -> MNum
S3Kerr = Simplify[S3  /. b -> bKerr]
S3KerrN = S3Kerr  /. a -> aNum  /. M -> MNum

COMPARE METRIC COMPONENTS TO NUMERICAL METRIC AND KERR METRIC

Compute metric components for analytic metric

gttN = gtt /. bm  /. a -> aNum  /. M -> MNum;
gtphiN  = gtphi  /. bm  /. a -> aNum  /. 
          M -> MNum;
grhorhoN  = grhorho  /. bm  /. a -> aNum  /. 
          M -> MNum;
gphiphiN  = gphiphi  /. bm  /. a -> aNum  /. 
          M -> MNum;

Read numerical metric in equatorial plane and along symmetry axis
input is from numerical code RNS
re is the numerical equatorial coordinate radius
Le is a factor that converts omega to geometrized units

SetDirectory["C:\DOCUMENTS-XP\Research\Projects\ExteriorManko\N\
umericalModels"];
FileNames[]

numtable = 
        ReadList["out_FPS.dat",  Number,  
          RecordLists -> True];
re = Import["out_FPS_re.dat"];
re = re[[1, 1]];
SDIV = 261;
Le = 29.979;

rhoEqTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 2}];
gammaEqTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 3, 2}];
alphaEqTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 4, 3}];
omegaEqTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 5, 4}];
rhoAxTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 6, 5}];
gammaAxTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 7, 6}];
alphaAxTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 8, 7}];
omegaAxTab = 
        Table[numtable[[i, j]], {i, SDIV}, {j, 1, 9, 8}];

sList = 
        Table[rhoEqTab[[i, 1]], {i, SDIV}];
rList = 
        Table[re*sList[[i]]/(1.0 - sList[[i]]), {i, 
            SDIV}];
RhoList = 
        Table[Exp[gammaEqTab[[i, 2]]]*rList[[i]], {i, 
            SDIV}];
gttEqList = 
        Table[-Exp[
                gammaEqTab[[i, 2]] + 
                  rhoEqTab[[i, 2]]] + (omegaEqTab[[i, 2]]/Le)^2*
              Exp[(gammaEqTab[[i, 2]] - rhoEqTab[[i, 2]])]*
              rList[[i]]^2, {i, SDIV}];
gtphiEqList = 
        Table[-(omegaEqTab[[i, 2]]/Le)*
            Exp[(gammaEqTab[[i, 2]] - rhoEqTab[[i, 2]])]*
            rList[[i]]^2, {i, SDIV}];
gphiphiEqList = 
        Table[Exp[(gammaEqTab[[i, 2]] - rhoEqTab[[i, 2]])]*
            rList[[i]]^2, {i, SDIV}];

gttEq = 
        Table[{RhoList[[i]], gttEqList[[i]]}, {i, 1, 
            SDIV}];
gtphiEq = 
        Table[{RhoList[[i]], gtphiEqList[[i]]}, {i, 1, 
            SDIV}];
gphiphiEq = 
        Table[{RhoList[[i]], gphiphiEqList[[i]]}, {i, 1, 
            SDIV}];

Compute z-coordinate and g\(TraditionalForm\`\_tt\) along symmetry axis

BAxList = 
        Table[{rList[[i]], Exp[gammaAxTab[[i, 2]]]}, {i, 1, 
            SDIV - 4}];
zList = Table[Null, {SDIV - 4}];
zList[[1]] = 0.0;
Do[
        zList[[i + 1]] = 
          zList[[i]] + 
            0.5*(rList[[i + 1]] - 
                  rList[[i]])*(BAxList[[i, 2]] + 
                  BAxList[[i + 1, 2]]), {i, 1, 
          SDIV - 5}];

gttAxList = 
        Table[-Exp[gammaAxTab[[i, 2]] + rhoAxTab[[i, 2]]], {i, 1, 
            SDIV - 4}];
gttAx = 
        Table[{zList[[i]], gttAxList[[i]]}, {i, 1, 
            SDIV - 4}];

Find location of surface on symmetry axis

0.6*re
rList[[99]]
zList[[99]]

Compute relative errors between analytic and numerical and between Kerr and numerical metric components

gttN  /. z -> 0  /. 
      rho  ->  RhoList[[131]]
gttN /. rho -> 0  /.  
      z  ->  zList[[99]]

gttEqEr = 
        Table[{RhoList[[i]], 
            Abs[(gttN  - gttEqList[[i]])/gttEqList[[i]]  /. 
                  z -> 0  /. rho  ->  RhoList[[i]]]}, {i, 
            131, SDIV}];
gttEqErK = 
        Table[{RhoList[[i]], 
            Abs[(gttK  - gttEqList[[i]])/gttEqList[[i]]  /. 
                  z -> 0  /. rho  ->  RhoList[[i]]]}, {i, 
            131, SDIV}];
gtphiEqEr = 
        Table[{RhoList[[i]], 
            Abs[(gtphiN  - gtphiEqList[[i]])/
                    gtphiEqList[[i]]  /. z -> 0  /. 
                rho  ->  RhoList[[i]]]}, {i, 131, 
            SDIV}];
gtphiEqErK = 
        Table[{RhoList[[i]], 
            Abs[(gtphiK  - gtphiEqList[[i]])/
                    gtphiEqList[[i]]  /. z -> 0  /. 
                rho  ->  RhoList[[i]]]}, {i, 131, 
            SDIV}];
gphiphiEqEr = 
        Table[{RhoList[[i]], 
            Abs[(gphiphiN  - gphiphiEqList[[i]])/
                    gphiphiEqList[[i]]  /. z -> 0  /. 
                rho  ->  RhoList[[i]]]}, {i, 131, 
            SDIV}];
gphiphiEqErK = 
        Table[{RhoList[[i]], 
            Abs[(gphiphiK  - gphiphiEqList[[i]])/
                    gphiphiEqList[[i]]  /. z -> 0  /. 
                rho  ->  RhoList[[i]]]}, {i, 131, 
            SDIV}];
gttAxEr = 
        Table[{zList[[i]], 
            Abs[(gttN  - gttAxList[[i]])/gttAxList[[i]]  /. 
                  rho -> 0  /. z  ->  zList[[i]]]}, {i, 
            99, SDIV - 4}];
gttAxErK = 
        Table[{zList[[i]], 
            Abs[(gttK  - gttAxList[[i]])/gttAxList[[i]]  /. 
                  rho -> 0  /. z  ->  zList[[i]]]}, {i, 
            99, SDIV - 4}];

Show relative errors

e1 = 
        ListPlot[Log[10, gttEqEr],  PlotStyle  -> RGBColor[0, 0, 1], 
          DisplayFunction :> 
            Identity];
e2 = 
        ListPlot[Log[10, gttEqErK],  PlotStyle  -> RGBColor[0, 1, 0],  
          DisplayFunction :> 
            Identity];
Show[e1, e2, DisplayFunction :> $DisplayFunction,  
        plotoptions, 
        FrameLabel -> {"Log rho", "Log Error"}];

e1 = 
        ListPlot[Log[10, gphiphiEqEr],  PlotStyle  -> RGBColor[0, 0, 1], 
          DisplayFunction :> 
            Identity];
e2 = 
        ListPlot[Log[10, gphiphiEqErK],  PlotStyle  -> RGBColor[0, 1, 0], 
          DisplayFunction :> 
            Identity];
Show[e1, e2, DisplayFunction :> $DisplayFunction,  
        plotoptions, 
        FrameLabel -> {"Log rho", "Log Error"}];

e1 = 
        ListPlot[Log[10, gtphiEqEr],  PlotStyle  -> RGBColor[0, 0, 1], 
          DisplayFunction :> 
            Identity];
e2 = 
        ListPlot[Log[10, gtphiEqErK],  PlotStyle  -> RGBColor[0, 1, 0], 
          DisplayFunction :> 
            Identity];
Show[e1, e2, DisplayFunction :> $DisplayFunction,  
        plotoptions, 
        FrameLabel -> {"Log rho", "Log Error"}];

gttAxEr[[77]]
gttAxErK[[70]]

e1 = 
        ListPlot[Log[10, gttAxEr],  PlotStyle  -> RGBColor[0, 0, 1], 
          DisplayFunction ->  Identity];
e2 = 
        ListPlot[Log[10, gttAxErK], PlotStyle  -> RGBColor[0, 1, 0], 
          DisplayFunction ->  Identity];
Show[e1, e2, DisplayFunction ->  $DisplayFunction,  
        plotoptions, 
        FrameLabel -> {"Log z", "Log Error"}];