{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Radial Oscillation Solver\n",
"\n",
"## Authors: Phil Chang\n",
"\n",
"Can work for any TOV star.\n",
"\n",
"[comment]: <> (Introduction: TODO)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This notebook is organized as follows:\n",
"\n",
"1. [Step 1](#initializenrpy): Initialize core Python/NRPy+ modules\n",
"1. [Step 2](#initializetov): Make a Model TOV Star using `TOV.TOV_Solver` NRPy+ module\n",
"1. [Step 3](#oscillations): Radial Oscillations\n",
"1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
"$$\\label{initializenrpy}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:14.029446Z",
"iopub.status.busy": "2021-03-07T17:15:14.028658Z",
"iopub.status.idle": "2021-03-07T17:15:14.267154Z",
"shell.execute_reply": "2021-03-07T17:15:14.266452Z"
}
},
"outputs": [],
"source": [
"# Step 1: Import needed Python/NRPy+ modules\n",
"import numpy as np # NumPy: A numerical methods module for Python\n",
"import scipy.integrate as si # SciPy: Python module for mathematics, science, and engineering applications\n",
"import math, sys # Standard Python modules for math; multiplatform OS-level functions\n",
"import TOV.Polytropic_EOSs as ppeos # NRPy+: Piecewise polytrope equation of state support"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: Make a Model TOV Star using `TOV.TOV_Solver` NRPy+ module \\[Back to [top](#toc)\\]\n",
"$$\\label{initializetov}$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:14.273702Z",
"iopub.status.busy": "2021-03-07T17:15:14.272920Z",
"iopub.status.idle": "2021-03-07T17:15:14.429269Z",
"shell.execute_reply": "2021-03-07T17:15:14.429740Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1256 1256 1256 1256 1256 1256\n",
"Just generated a TOV star with\n",
"* M = 1.405030336771405e-01 ,\n",
"* R_Schw = 9.566044579232513e-01 ,\n",
"* R_iso = 8.100085557410308e-01 ,\n",
"* M/R_Schw = 1.468768334847266e-01 \n",
"\n"
]
}
],
"source": [
"import TOV.TOV_Solver as TOV\n",
"\n",
"############################\n",
"# Single polytrope example #\n",
"############################\n",
"# Set neos = 1 (single polytrope)\n",
"neos = 1\n",
"\n",
"# Set rho_poly_tab (not needed for a single polytrope)\n",
"rho_poly_tab = []\n",
"\n",
"# Set Gamma_poly_tab\n",
"Gamma_poly_tab = [2.0]\n",
"\n",
"# Set K_poly_tab0\n",
"K_poly_tab0 = 1. # ZACH NOTES: CHANGED FROM 100.\n",
"\n",
"# Set the eos quantities\n",
"eos = ppeos.set_up_EOS_parameters__complete_set_of_input_variables(neos,rho_poly_tab,Gamma_poly_tab,K_poly_tab0)\n",
"\n",
"# Set initial condition (Pressure computed from central density)\n",
"rho_baryon_central = 0.129285\n",
"\n",
"# output file\n",
"tov_file = \"outputTOVpolytrope-oscillations.txt\"\n",
"\n",
"TOV.TOV_Solver(eos,\n",
" outfile=tov_file,\n",
" rho_baryon_central=rho_baryon_central,\n",
" verbose = True,\n",
" accuracy=\"medium\",\n",
" integrator_type=\"default\",\n",
" no_output_File = False)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 3: Radial Oscillations \\[Back to [top](#toc)\\]\n",
"$$\\label{oscillations}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"According to [Bardeen, Thorne and Meltzer (1966; hereafter BTM)](https://ui.adsabs.harvard.edu/abs/1966ApJ...145..505B/abstract), the line element is written as\n",
"$$\n",
"ds^2 = e^\\nu dt^2 - e^{\\lambda} dr^2 - r^2 d\\Omega^2,\n",
"$$\n",
"The metric that we defined previously is (for $G=c=1$ and not different sign convention)\n",
"$$\n",
"ds^2 = -e^\\nu dt^2 + \\left(1 - \\frac{2m}{r^2}\\right)^{-1} dr^2 + r^2 d\\Omega^2,\n",
"$$\n",
"So the identification is \n",
"$$\n",
"e^{\\lambda} = \\left(1 - \\frac{2m}{r^2}\\right)^{-1}\n",
"$$\n",
"BTM states that the eigenfunction $u_n$ for the displacement of the form\n",
"$$\n",
"\\delta r = e^{\\lambda/2}r^{-2} u_n e^{i\\omega_n t}\n",
"$$\n",
"satisfies \n",
"$$\n",
"\\frac{d}{dr}P\\frac{du_n}{dr} + \\left(Q + \\omega_n^2 W\\right)u_n = 0\n",
"$$\n",
"where\n",
"$$\n",
"P = e^{(\\lambda + 3\\nu)/2} r^{-2}\\gamma p \\\\\n",
"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}\\\\\n",
"W = e^{(3\\lambda + \\nu)/2}r^{-2} (p + \\rho)\n",
"$$\n",
"where $\\gamma$ is the adiabatic index\n",
"\n",
"The solution to the Sturm-Liouville problem must satisfy a BC at $r=0$ and at $r=R$. In particular at $r\\rightarrow 0$\n",
"$$\n",
"u_n = r^3\n",
"$$\n",
"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\n",
"$$\n",
"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} \n",
"$$\n",
"\n",
"So to find the eigenvalues, $\\omega_n^2$, integrate from $r=0$ outward $r=R$ or until it start blowing up.\n",
"\n",
"For convenience we will define $\\Phi = e^{(\\lambda + 3\\nu)/2} r^{-2}$\n",
"$$\n",
"P = \\Phi\\gamma p \\\\\n",
"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\\\\\n",
"W = \\Phi e^{\\lambda-\\nu}(p + \\rho)\n",
"$$\n",
"\n",
"The formal equations that we are integrating is then\n",
"$$\n",
"\\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\n",
"$$\n",
"\n",
"The alternative way is to write\n",
"$$\n",
"\\frac{d\\psi}{dr} = -\\frac 1 r \\left(3\\psi + \\frac{\\Delta p}{\\Gamma p}\\right) - \\frac {dp}{dr}\\frac{\\psi}{p+\\rho} \\\\\n",
"\\frac{d\\Delta p}{dr} = \\psi\\left(\\omega^2 e^{\\lambda-\\nu} (p+\\rho)r - 4 \\frac {dp}{dr}\\right) +\n",
"\\psi\\left(\\left(\\frac{dp}{dr}\\right)^2\\frac r {p+\\rho} - 8\\pi e^\\lambda(p+\\rho)p r\\right) + \n",
"\\Delta p\\left(\\frac{dp}{dr}\\frac 1 {p+\\rho} - 4\\pi (p+\\rho)r e^\\lambda\\right)\n",
"$$\n",
"The BC of this equations are $\\psi(0) = 1$, $\\Delta p(0) = -3\\psi\\Gamma p(0)$, and $\\Delta p(R) = 0$\n",
"\n",
"The one thing that is confusing is the relativistic adiabatic index \n",
"$$\n",
"\\Gamma = \\frac{\\rho}{P}\\frac{dP}{d\\rho}\n",
"$$\n",
"Now this is in terms of $\\rho_b$\n",
"$$\n",
"\\Gamma = \\frac{\\rho}{P}\\frac{dP}{d\\rho_b}\\frac{d\\rho_b}{d\\rho} \n",
"$$\n",
"$$\n",
"\\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}}\\\\\n",
"= \\Gamma_b \\frac{1 + \\rho_b^{\\Gamma_b - 1}}{1 + \\Gamma_b\\rho_b^{\\Gamma_b - 1}}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:14.469959Z",
"iopub.status.busy": "2021-03-07T17:15:14.469211Z",
"iopub.status.idle": "2021-03-07T17:15:33.448860Z",
"shell.execute_reply": "2021-03-07T17:15:33.448157Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.08092993644066569 0.0009266974129402323\n",
"0.08488003342081815 0.0010534195609376677\n",
"0.0886543035320239 0.0011630368590904265\n",
"0.09227432468448224 0.0012567190261121282\n",
"0.09575759216483186 0.0013355600114076866\n",
"0.0991185245977526 0.0014005803879500699\n",
"0.10236917201806559 0.0014528638599218846\n",
"0.10551972725937707 0.0014933026512813181\n",
"0.10857890357763954 0.0015228410503980794\n",
"0.11155421894016991 0.0015422657550877418\n",
"0.114452213716382 0.0015524689447448277\n",
"0.11727861990068524 0.0015541667704478434\n",
"0.12003849443840264 0.0015481439141351672\n",
"0.12273632554491537 0.0015350109880955482\n",
"0.12537611841772922 0.0015154807035301656\n",
"0.12796146502258202 0.0014900987656813593\n",
"0.13049560142761327 0.0014595027654243577\n",
"0.1329814552980359 0.0014242760481190818\n",
"0.13542168553969686 0.0013848889245836665\n",
"0.1378187156217963 0.0013418016015926024\n",
"0.14017476176855295 0.0012954864707387785\n",
"0.1424918569536505 0.0012463500494504487\n",
"0.14477187143685272 0.00119483628570832\n",
"0.14701653043300142 0.0011413168211486969\n",
"0.14922742938812106 0.0010860953142048473\n",
"0.15140604724718124 0.0010295503375526342\n",
"0.15355375802709842 0.0009719603842770745\n",
"0.15567184095228564 0.0009136211405781264\n",
"0.1577614893651259 0.0008547500890532673\n",
"0.15982381858763026 0.0007956885433670953\n",
"0.16185987288133144 0.0007365663415519544\n",
"0.16387063162870216 0.0006776646590828148\n",
"0.16585701483994822 0.0006192068001031311\n",
"0.16781988807304046 0.0005612653013117582\n",
"0.1697600668416364 0.0005039817206930438\n",
"0.17167832057457305 0.000447645800614521\n",
"0.17357537618145866 0.00039228152877463613\n",
"0.17545192127122436 0.00033804054909160694\n",
"0.17730860706404789 0.0002850475631361193\n",
"0.17914605103161482 0.00023337927428010373\n",
"0.1809648392960659 0.0001831634977604156\n",
"0.18276552881405175 0.00013443055187650286\n",
"0.1845486493689645 8.72758764257258e-05\n",
"0.18631470539154452 4.1739138479586896e-05\n",
"0.18806417762659386 -2.1260732318257964e-06\n",
"0.18979752466140346 -4.425022931862742e-05\n",
"0.19151518432966375 -8.463706436586268e-05\n",
"0.19321757500303668 -0.00012324790006858427\n",
"0.1949050967811851 -0.0001600492128691676\n",
"0.19657813258984871 -0.00019501596936231962\n",
"0.19823704919550525 -0.00022817829593481192\n",
"0.1998821981442324 -0.00025950631348657156\n",
"0.20151391663157958 -0.0002890229613744426\n",
"0.20313252830954528 -0.00031671358971690306\n",
"0.20473834403613125 -0.0003425904826522005\n",
"0.20633166257238916 -0.0003666913404417286\n",
"0.2079127712313875 -0.0003890137421460794\n",
"0.20948194648309026 -0.00040957203178132283\n",
"0.21103945451875417 -0.0004284048570553927\n",
"0.21258555177810728 -0.000445536048601437\n",
"0.2141204854422653 -0.00046098211509420654\n",
"0.21564449389506882 -0.00047481904952238014\n",
"0.21715780715527913 -0.0004870474551833342\n",
"0.2186606472818518 -0.0004977307540138093\n",
"0.22015322875430973 -0.00050688913949875\n",
"0.18806417762659386\n"
]
}
],
"source": [
"\n",
"import scipy.integrate as si\n",
"\n",
"def deriv( r, y, rArr_np, P_np, dPdr_np, Q_np, W_np, omega2):\n",
" un = y[0]\n",
" dundr = y[1]\n",
"\n",
" P = np.interp(r, rArr_np, P_np)\n",
" Q = np.interp(r, rArr_np, Q_np)\n",
" W = np.interp(r, rArr_np, W_np)\n",
" dPdr = np.interp(r, rArr_np, dPdr_np)\n",
"\n",
" d2undr2 = -(dundr*dPdr/P + (Q + omega2*W)*un/P)\n",
" #print((Q + omega2*W)/Q)\n",
" return [dundr, d2undr2]\n",
"\n",
"def gamma_rel( r, rArr_np, rho_baryon_np, gammab=2) :\n",
" rho_b = np.interp(r, rArr_np, rho_baryon_np)\n",
" P = rho_b**gammab\n",
" return gammab*(1+rho_b**(gammab-1))/(1+gammab*rho_b**(gammab-1))\n",
" return (rho_b + P)/P * gammab*rho_b**(gammab - 1)/(1+gammab*rho_b**(gammab-1))\n",
"\n",
"def alt_deriv(r, y, rArr_np, eLambda_np, eNu_np, p_np, rho_np, dpdr_np, rho_baryon_np, m_np, omega2) :\n",
" psi = y[0]\n",
" deltap = y[1]\n",
" p = np.interp(r, rArr_np, p_np)\n",
" rho = np.interp(r, rArr_np, rho_np)\n",
" dpdr = np.interp(r, rArr_np, dpdr_np)\n",
" eLambda = np.interp(r, rArr_np, eLambda_np)\n",
" eNu = np.interp(r, rArr_np, eNu_np)\n",
" m = np.interp(r, rArr_np, m_np)\n",
" rhoPlusP = p + rho\n",
" gamma = gamma_rel( r, rArr_np, rho_baryon_np)\n",
" dPdrSchw = -(rho + p)*(m + 4.*math.pi*r**3*p)/(r*r*(1.-2.*m/r))\n",
" dpdr = dPdrSchw\n",
"\n",
" dpsidr = -(3*psi + deltap/(gamma*p))/r - dpdr*psi/rhoPlusP\n",
" #eNu = 1\n",
" #eLambda = 1\n",
" ddeltapdr = psi*(omega2*eLambda/eNu*rhoPlusP*r - 4*dpdr) + \\\n",
" psi*(dpdr**2*r/rhoPlusP - 8*np.pi*eLambda*rhoPlusP*p*r) + \\\n",
" deltap*(dpdr/rhoPlusP - 4*np.pi*rhoPlusP*r*eLambda)\n",
" return [dpsidr, ddeltapdr]\n",
"\n",
"# read in the initial data\n",
"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)\n",
"\n",
"boolArray = rhoArr_np > 0\n",
"expnu = exp2phiArr_np[boolArray]\n",
"explambda = 1./(1- 2*mArr_np/r_SchwArr_np)[boolArray]\n",
"\n",
"r_SchwArr_np = r_SchwArr_np[boolArray]\n",
"rhoArr_np = rhoArr_np[boolArray]\n",
"PArr_np = PArr_np[boolArray]\n",
"rho_baryonArr_np = rho_baryonArr_np[boolArray]\n",
"mArr_np = mArr_np[boolArray]\n",
"dpdr_np = np.zeros(PArr_np.size)\n",
"dpdr_np[:-1] = (PArr_np[1:]-PArr_np[:-1])/(r_SchwArr_np[1:] - r_SchwArr_np[:-1])\n",
"dpdr_np[-1] = dpdr_np[-2] # copy over the last one\n",
"\n",
"min_dDeltaP = 1\n",
"min_omega2 = 0\n",
"for n in np.arange(2.0, 15.0, 0.2) :\n",
"\n",
"\n",
" omega2 = n*rho_baryonArr_np[0]\n",
"\n",
" r = si.ode(alt_deriv).set_integrator('dopri5')\n",
"\n",
" # integrate out a little bit\n",
" r0 = 1e-3\n",
"\n",
" psi = 1\n",
" p = np.interp(r0, r_SchwArr_np, PArr_np)\n",
" gamma = gamma_rel( r0, r_SchwArr_np, rho_baryonArr_np)\n",
" #print(gamma)\n",
" y0 = [psi, -3*gamma*psi*p]\n",
"\n",
" #r.set_initial_value(y0, r0).set_f_params(r_SchwArr_np, P_np, dPdr_np, Q_np, W_np, omega2)\n",
" 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)\n",
"\n",
" R = r_SchwArr_np[-1]\n",
" dr = 1e-2\n",
"\n",
" r_arr = []\n",
" psi_arr = []\n",
" Deltap_arr = []\n",
"\n",
" while r.successful() and r.t < R:\n",
" dr = min(max(r.t*1e-2, 1e-5), R-r.t)\n",
" r_arr.append(r.t+dr)\n",
" psi, deltap = r.integrate(r.t+dr)\n",
" psi_arr.append(psi)\n",
" Deltap_arr.append(deltap)\n",
"\n",
" r_arr = np.array(r_arr)\n",
" Deltap_arr = np.array( Deltap_arr)\n",
" print( math.sqrt(omega2)/(2*math.pi), Deltap_arr[-1])\n",
"\n",
" if( min_dDeltaP > np.abs(Deltap_arr[-1])) :\n",
" min_omega2 = omega2\n",
" min_dDeltaP = np.abs(Deltap_arr[-1])\n",
"\n",
" #pl.plot(r_arr, Deltap_arr)\n",
"\n",
"print( (min_omega2)**0.5/(2*math.pi))\n",
"#pl.show()\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}