Notebook

Radial Oscillation Solver¶

Authors: Phil Chang¶

Can work for any TOV star.

Table of Contents¶

$\label{toc}$

This notebook is organized as follows:

Step 1: Initialize core Python/NRPy+ modules
Step 2: Make a Model TOV Star using TOV.TOV_Solver NRPy+ module
Step 3: Radial Oscillations
Step 4: Output this notebook to $\LaTeX$ -formatted PDF file

Step 1: Initialize core Python/NRPy+ modules [Back to top]¶

$\label{initializenrpy}$

In [1]:

# Step 1: Import needed Python/NRPy+ modules
import numpy as np                  # NumPy: A numerical methods module for Python
import scipy.integrate as si        # SciPy: Python module for mathematics, science, and engineering applications
import math, sys                    # Standard Python modules for math; multiplatform OS-level functions
import TOV.Polytropic_EOSs as ppeos # NRPy+: Piecewise polytrope equation of state support

Step 2: Make a Model TOV Star using `TOV.TOV_Solver` NRPy+ module [Back to top]¶

$\label{initializetov}$

In [2]:

import TOV.TOV_Solver as TOV

############################
# Single polytrope example #
############################
# Set neos = 1 (single polytrope)
neos = 1

# Set rho_poly_tab (not needed for a single polytrope)
rho_poly_tab = []

# Set Gamma_poly_tab
Gamma_poly_tab = [2.0]

# Set K_poly_tab0
K_poly_tab0 = 1. # ZACH NOTES: CHANGED FROM 100.

# Set the eos quantities
eos = ppeos.set_up_EOS_parameters__complete_set_of_input_variables(neos,rho_poly_tab,Gamma_poly_tab,K_poly_tab0)

# Set initial condition (Pressure computed from central density)
rho_baryon_central  = 0.129285

# output file
tov_file = "outputTOVpolytrope-oscillations.txt"

TOV.TOV_Solver(eos,
               outfile=tov_file,
               rho_baryon_central=rho_baryon_central,
               verbose = True,
               accuracy="medium",
               integrator_type="default",
               no_output_File = False)

1256 1256 1256 1256 1256 1256
Just generated a TOV star with
* M        = 1.405030336771405e-01 ,
* R_Schw   = 9.566044579232513e-01 ,
* R_iso    = 8.100085557410308e-01 ,
* M/R_Schw = 1.468768334847266e-01

Step 3: Radial Oscillations [Back to top]¶

$\label{oscillations}$

According to Bardeen, Thorne and Meltzer (1966; hereafter BTM), the line element is written as

$ds^2 = e^\nu dt^2 - e^{\lambda} dr^2 - r^2 d\Omega^2,$ The metric that we defined previously is (for

$G=c=1$ and not different sign convention)

$ds^2 = -e^\nu dt^2 + \left(1 - \frac{2m}{r^2}\right)^{-1} dr^2 + r^2 d\Omega^2,$ So the identification is

$e^{\lambda} = \left(1 - \frac{2m}{r^2}\right)^{-1}$ BTM states that the eigenfunction

$u_n$ for the displacement of the form

$\delta r = e^{\lambda/2}r^{-2} u_n e^{i\omega_n t}$ satisfies

$\frac{d}{dr}P\frac{du_n}{dr} + \left(Q + \omega_n^2 W\right)u_n = 0$ where

$P = e^{(\lambda + 3\nu)/2} r^{-2}\gamma p \\ Q = -4e^{(\lambda + 3\nu)/2} r^{-3}\frac{dp}{dr} - 8\pi e^{3(\nu+\lambda)/2}r^{-2}p(p+\rho) + e^{(\lambda + 3\nu)/2} r^{-2}(p + \rho)^{-1}\frac{dp}{dr}\\ W = e^{(3\lambda + \nu)/2}r^{-2} (p + \rho)$ where

$\gamma$ is the adiabatic index

The solution to the Sturm-Liouville problem must satisfy a BC at $r=0$ and at $r=R$ . In particular at $r\rightarrow 0$

$u_n = r^3$ At

$r=R$ , we have the pressure perturbation is zero. This gives a singular point at

$r=R$ and the blowup is of the form

$u_n \propto (R-r)^{-s} \quad\textrm{and}\quad s = \left[\frac{d\ln p/dr}{d\ln (p/\rho)/dr}\right]_{r=R}$

So to find the eigenvalues, $\omega_n^2$ , integrate from $r=0$ outward $r=R$ or until it start blowing up.

For convenience we will define $\Phi = e^{(\lambda + 3\nu)/2} r^{-2}$

$P = \Phi\gamma p \\ Q = -4\Phi r^{-1}\frac{dp}{dr} - 8\pi e^{\lambda}\Phi p(p+\rho) + \Phi (p + \rho)^{-1}\left(\frac{dp}{dr}\right)^2\\ W = \Phi e^{\lambda-\nu}(p + \rho)$

The formal equations that we are integrating is then

$\frac{d^2u_n}{dr^2} + P^{-1}\frac{dP}{dr}\frac{du_n}{dr} + P^{-1}(Q + \omega_n^2 W)u_n = 0$

The alternative way is to write

$\frac{d\psi}{dr} = -\frac 1 r \left(3\psi + \frac{\Delta p}{\Gamma p}\right) - \frac {dp}{dr}\frac{\psi}{p+\rho} \\ \frac{d\Delta p}{dr} = \psi\left(\omega^2 e^{\lambda-\nu} (p+\rho)r - 4 \frac {dp}{dr}\right) + \psi\left(\left(\frac{dp}{dr}\right)^2\frac r {p+\rho} - 8\pi e^\lambda(p+\rho)p r\right) + \Delta p\left(\frac{dp}{dr}\frac 1 {p+\rho} - 4\pi (p+\rho)r e^\lambda\right)$ The BC of this equations are

$\psi(0) = 1$ ,

$\Delta p(0) = -3\psi\Gamma p(0)$ , and

$\Delta p(R) = 0$

The one thing that is confusing is the relativistic adiabatic index

$\Gamma = \frac{\rho}{P}\frac{dP}{d\rho}$ Now this is in terms of

$\rho_b$

$\Gamma = \frac{\rho}{P}\frac{dP}{d\rho_b}\frac{d\rho_b}{d\rho}$

$\Gamma = \frac{\rho_b + \rho_b^{\Gamma_b}}{\rho_b^{\Gamma_b}}\Gamma_b\rho_b^{\Gamma_b -1} \frac{1}{1 + \Gamma_b\rho_b^{\Gamma_b - 1}}\\ = \Gamma_b \frac{1 + \rho_b^{\Gamma_b - 1}}{1 + \Gamma_b\rho_b^{\Gamma_b - 1}}$

In [3]:

import scipy.integrate as si

def deriv( r, y, rArr_np, P_np, dPdr_np, Q_np, W_np, omega2):
    un = y[0]
    dundr = y[1]

    P = np.interp(r, rArr_np, P_np)
    Q = np.interp(r, rArr_np, Q_np)
    W = np.interp(r, rArr_np, W_np)
    dPdr = np.interp(r, rArr_np, dPdr_np)

    d2undr2 = -(dundr*dPdr/P + (Q + omega2*W)*un/P)
    #print((Q + omega2*W)/Q)
    return [dundr, d2undr2]

def gamma_rel( r, rArr_np, rho_baryon_np, gammab=2) :
    rho_b = np.interp(r, rArr_np, rho_baryon_np)
    P = rho_b**gammab
    return gammab*(1+rho_b**(gammab-1))/(1+gammab*rho_b**(gammab-1))
    return (rho_b + P)/P * gammab*rho_b**(gammab - 1)/(1+gammab*rho_b**(gammab-1))

def alt_deriv(r, y, rArr_np, eLambda_np, eNu_np, p_np, rho_np, dpdr_np, rho_baryon_np, m_np, omega2) :
    psi = y[0]
    deltap = y[1]
    p = np.interp(r, rArr_np, p_np)
    rho = np.interp(r, rArr_np, rho_np)
    dpdr = np.interp(r, rArr_np, dpdr_np)
    eLambda = np.interp(r, rArr_np, eLambda_np)
    eNu = np.interp(r, rArr_np, eNu_np)
    m = np.interp(r, rArr_np, m_np)
    rhoPlusP = p + rho
    gamma = gamma_rel( r, rArr_np, rho_baryon_np)
    dPdrSchw = -(rho + p)*(m + 4.*math.pi*r**3*p)/(r*r*(1.-2.*m/r))
    dpdr = dPdrSchw

    dpsidr = -(3*psi + deltap/(gamma*p))/r - dpdr*psi/rhoPlusP
    #eNu = 1
    #eLambda = 1
    ddeltapdr = psi*(omega2*eLambda/eNu*rhoPlusP*r - 4*dpdr) + \
        psi*(dpdr**2*r/rhoPlusP - 8*np.pi*eLambda*rhoPlusP*p*r) + \
        deltap*(dpdr/rhoPlusP - 4*np.pi*rhoPlusP*r*eLambda)
    return [dpsidr, ddeltapdr]

# read in the initial data
r_SchwArr_np,rhoArr_np,rho_baryonArr_np,PArr_np,mArr_np,exp2phiArr_np,confFactor_exp4phi_np,rbarArr_np = np.loadtxt(tov_file, unpack=True)

boolArray = rhoArr_np > 0
expnu = exp2phiArr_np[boolArray]
explambda = 1./(1- 2*mArr_np/r_SchwArr_np)[boolArray]

r_SchwArr_np = r_SchwArr_np[boolArray]
rhoArr_np = rhoArr_np[boolArray]
PArr_np = PArr_np[boolArray]
rho_baryonArr_np = rho_baryonArr_np[boolArray]
mArr_np = mArr_np[boolArray]
dpdr_np = np.zeros(PArr_np.size)
dpdr_np[:-1] = (PArr_np[1:]-PArr_np[:-1])/(r_SchwArr_np[1:] - r_SchwArr_np[:-1])
dpdr_np[-1] = dpdr_np[-2] # copy over the last one

min_dDeltaP = 1
min_omega2 = 0
for n in np.arange(2.0, 15.0, 0.2) :


    omega2 = n*rho_baryonArr_np[0]

    r = si.ode(alt_deriv).set_integrator('dopri5')

    # integrate out a little bit
    r0 = 1e-3

    psi = 1
    p = np.interp(r0, r_SchwArr_np, PArr_np)
    gamma = gamma_rel( r0, r_SchwArr_np, rho_baryonArr_np)
    #print(gamma)
    y0 = [psi, -3*gamma*psi*p]

    #r.set_initial_value(y0, r0).set_f_params(r_SchwArr_np, P_np, dPdr_np, Q_np, W_np, omega2)
    r.set_initial_value(y0, r0).set_f_params(r_SchwArr_np, explambda, expnu, PArr_np, rhoArr_np, dpdr_np,rho_baryonArr_np,mArr_np, omega2)

    R = r_SchwArr_np[-1]
    dr = 1e-2

    r_arr = []
    psi_arr = []
    Deltap_arr = []

    while r.successful() and r.t < R:
        dr = min(max(r.t*1e-2, 1e-5), R-r.t)
        r_arr.append(r.t+dr)
        psi, deltap = r.integrate(r.t+dr)
        psi_arr.append(psi)
        Deltap_arr.append(deltap)

    r_arr = np.array(r_arr)
    Deltap_arr = np.array( Deltap_arr)
    print( math.sqrt(omega2)/(2*math.pi), Deltap_arr[-1])

    if( min_dDeltaP > np.abs(Deltap_arr[-1])) :
        min_omega2 = omega2
        min_dDeltaP = np.abs(Deltap_arr[-1])

    #pl.plot(r_arr, Deltap_arr)

print( (min_omega2)**0.5/(2*math.pi))
#pl.show()

0.08092993644066569 0.0009266974129402323
0.08488003342081815 0.0010534195609376677
0.0886543035320239 0.0011630368590904265
0.09227432468448224 0.0012567190261121282
0.09575759216483186 0.0013355600114076866
0.0991185245977526 0.0014005803879500699
0.10236917201806559 0.0014528638599218846
0.10551972725937707 0.0014933026512813181
0.10857890357763954 0.0015228410503980794
0.11155421894016991 0.0015422657550877418
0.114452213716382 0.0015524689447448277
0.11727861990068524 0.0015541667704478434
0.12003849443840264 0.0015481439141351672
0.12273632554491537 0.0015350109880955482
0.12537611841772922 0.0015154807035301656
0.12796146502258202 0.0014900987656813593
0.13049560142761327 0.0014595027654243577
0.1329814552980359 0.0014242760481190818
0.13542168553969686 0.0013848889245836665
0.1378187156217963 0.0013418016015926024
0.14017476176855295 0.0012954864707387785
0.1424918569536505 0.0012463500494504487
0.14477187143685272 0.00119483628570832
0.14701653043300142 0.0011413168211486969
0.14922742938812106 0.0010860953142048473
0.15140604724718124 0.0010295503375526342
0.15355375802709842 0.0009719603842770745
0.15567184095228564 0.0009136211405781264
0.1577614893651259 0.0008547500890532673
0.15982381858763026 0.0007956885433670953
0.16185987288133144 0.0007365663415519544
0.16387063162870216 0.0006776646590828148
0.16585701483994822 0.0006192068001031311
0.16781988807304046 0.0005612653013117582
0.1697600668416364 0.0005039817206930438
0.17167832057457305 0.000447645800614521
0.17357537618145866 0.00039228152877463613
0.17545192127122436 0.00033804054909160694
0.17730860706404789 0.0002850475631361193
0.17914605103161482 0.00023337927428010373
0.1809648392960659 0.0001831634977604156
0.18276552881405175 0.00013443055187650286
0.1845486493689645 8.72758764257258e-05
0.18631470539154452 4.1739138479586896e-05
0.18806417762659386 -2.1260732318257964e-06
0.18979752466140346 -4.425022931862742e-05
0.19151518432966375 -8.463706436586268e-05
0.19321757500303668 -0.00012324790006858427
0.1949050967811851 -0.0001600492128691676
0.19657813258984871 -0.00019501596936231962
0.19823704919550525 -0.00022817829593481192
0.1998821981442324 -0.00025950631348657156
0.20151391663157958 -0.0002890229613744426
0.20313252830954528 -0.00031671358971690306
0.20473834403613125 -0.0003425904826522005
0.20633166257238916 -0.0003666913404417286
0.2079127712313875 -0.0003890137421460794
0.20948194648309026 -0.00040957203178132283
0.21103945451875417 -0.0004284048570553927
0.21258555177810728 -0.000445536048601437
0.2141204854422653 -0.00046098211509420654
0.21564449389506882 -0.00047481904952238014
0.21715780715527913 -0.0004870474551833342
0.2186606472818518 -0.0004977307540138093
0.22015322875430973 -0.00050688913949875
0.18806417762659386