{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simulation of a Kerr-Black-Hole\n",
"by **Florian Hollants**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This Jupyter Notebook simulates a rotating Black Hole in SageMath using the Kerr Metric. The Project is based on the work of Florentin Jaffredo. The parameters for mass and angular momentum can be altered in Cell 5 by changing \"m_val\" and \"a_val\", while the \"Angle\" can be changed in Cell 18. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# n_cpu = 4 # 4 Go Ram minimum\n",
"# n_geod = 100\n",
"# nx, ny = 180, 90"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"n_cpu = 8 # 8 Go Ram minimum\n",
"n_geod = 1000\n",
"nx, ny = 720, 360"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# n_cpu = 36 # 144 Go Ram minimum\n",
"# n_geod = 30000\n",
"# nx, ny = 4000, 2000"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"%display latex"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}t :\\ \\left( -\\infty, +\\infty \\right) ;\\quad r :\\ \\left( 4.20000000000000 , +\\infty \\right) ;\\quad {\\theta} :\\ \\left( 0 , \\pi \\right) ;\\quad {\\phi} :\\ \\left( -\\infty, +\\infty \\right)\n",
"\\end{math}"
],
"text/plain": [
"t: (-oo, +oo); r: (4.20000000000000, +oo); th: (0, pi); ph: (-oo, +oo)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = Manifold(4, 'M', structure='Lorentzian')\n",
"\n",
"m = var('m')\n",
"a = var('a')\n",
"\n",
"m_val = 2\n",
"a_val = 1\n",
"\n",
"C. = M.chart(r't r:(4.2,+oo) th:(0,pi):\\theta ph:\\phi') #check that minimum radius is bigger than singularity in g[1,1]\n",
"C.coord_range()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{2 \\, m r}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} - 1 & 0 & 0 & -\\frac{2 \\, a m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} \\\\\n",
"0 & \\frac{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}}{a^{2} - 2 \\, m r + r^{2}} & 0 & 0 \\\\\n",
"0 & 0 & a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2} & 0 \\\\\n",
"-\\frac{2 \\, a m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} & 0 & 0 & {\\left(\\frac{2 \\, a^{2} m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} + a^{2} + r^{2}\\right)} \\sin\\left({\\theta}\\right)^{2}\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 2*m*r/(a^2*cos(th)^2 + r^2) - 1 0 0 -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2)]\n",
"[ 0 (a^2*cos(th)^2 + r^2)/(a^2 - 2*m*r + r^2) 0 0]\n",
"[ 0 0 a^2*cos(th)^2 + r^2 0]\n",
"[ -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) 0 0 (2*a^2*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) + a^2 + r^2)*sin(th)^2]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = M.metric()\n",
"g[0,0] = (2*m*r)/(r^2+(a*cos(th))^2)-1\n",
"g[0,3] = -2*m*r*a*sin(th)^2/(r^2+(a*cos(th))^2)\n",
"g[1,1] = (r^2+(a*cos(th))^2)/(r^2-2*m*r+a^2)\n",
"g[2,2] = r^2+(a*cos(th))^2\n",
"g[3,3] = (r^2+a^2+(2*m*r*a^2*sin(th)^2)/(r^2+(a*cos(th))^2))*sin(th)^2\n",
"g[:]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{E}^{3} \\\\ & \\left(t, r, {\\theta}, {\\phi}\\right) & \\longmapsto & \\left(x, y, z\\right) = \\left(r \\cos\\left({\\phi}\\right) \\sin\\left({\\theta}\\right), r \\sin\\left({\\phi}\\right) \\sin\\left({\\theta}\\right), r \\cos\\left({\\theta}\\right)\\right) \\end{array}\n",
"\\end{math}"
],
"text/plain": [
"M --> E^3\n",
" (t, r, th, ph) |--> (x, y, z) = (r*cos(ph)*sin(th), r*sin(ph)*sin(th), r*cos(th))"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E. = EuclideanSpace()\n",
"phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])\n",
"phi.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## One Geodesic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try the following for circular Orbit"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = M((0, 14.98, pi/2, 0))\n",
"Tp = M.tangent_space(p)\n",
"v = Tp((2, 0, 0.005, 0.05))\n",
"v = v / sqrt(-g.at(p)(v, v))\n",
"\n",
"tau = var('tau')\n",
"\n",
"curve = M.integrated_geodesic(g, (tau, 0, 5000), v)\n",
"sol = curve.solve(step = 1, method=\"ode_int\", parameters_values={m: m_val, a: a_val})\n",
"\n",
"interp = curve.interpolate()\n",
"\n",
"P = curve.plot_integrated(mapping=phi, color=\"red\", thickness=2, plot_points=5000)\n",
"P += sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try the following for light at different Values for a, aimed directly at the Black Hole"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# p = M((0, 20, pi/2, 0))\n",
"# Tp = M.tangent_space(p)\n",
"# v = Tp((2, -1, 0.000, 0.00))\n",
"# v = v / sqrt(-g.at(p)(v, v))\n",
"\n",
"# A=[0,1,2,4] #list of values for a\n",
"# Color=[[\"red\"], [\"blue\"], [\"green\"],[\"orange\"], [\"black\"]] #color for different Geodesics\n",
"\n",
"# P = sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"\n",
"# tau = var('tau')\n",
"# for i in range(4):\n",
"# curve = M.integrated_geodesic(g, (tau, 0, 200), initial_tangent_vector=v, across_charts=True)\n",
"# curve.solve_across_charts(step=0.2, parameters_values={m:m_val, a:A[i]})\n",
"# interp = curve.interpolate()\n",
"# P += curve.plot_integrated(mapping=phi, color=Color[i], thickness=2, plot_points=10000, across_charts=True)\n",
"\n",
"# P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multiple Geodesics"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import multiprocessing\n",
"from ipywidgets import FloatProgress\n",
"from IPython.display import display"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def chunks(l, n):\n",
" \"\"\"Yield successive n-sized chunks from l.\"\"\"\n",
" for i in range(0, len(l), n):\n",
" yield l[i:i + n]\n",
"\n",
"n_batches_per_cpu = 3"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"curve = M.integrated_geodesic(g, (tau, 0, 200), v, across_charts=True)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"args = []\n",
"start_index = 0\n",
"\n",
"for chunk in chunks(range(n_geod), n_geod//(n_batches_per_cpu*n_cpu)):\n",
" args += [(loads(curve.dumps()), start_index, len(chunk))] #\n",
" start_index += len(chunk)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\mathit{dt} = \\frac{2 \\, {\\left(a^{3} m - 2 \\, a m^{2} r_{0} + a m r_{0}^{2}\\right)} y - \\sqrt{-2 \\, a^{2} m r_{0} - 4 \\, m r_{0}^{3} + r_{0}^{4} + {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{2} + {\\left(a^{6} - 6 \\, a^{4} m r_{0} - 6 \\, m r_{0}^{5} + r_{0}^{6} + 3 \\, {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{4} - 4 \\, {\\left(3 \\, a^{2} m + 2 \\, m^{3}\\right)} r_{0}^{3} + 3 \\, {\\left(a^{4} + 4 \\, a^{2} m^{2}\\right)} r_{0}^{2}\\right)} y^{2}} r_{0}}{2 \\, a^{2} m + 4 \\, m r_{0}^{2} - r_{0}^{3} - {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}}, \\mathit{dt} = \\frac{2 \\, {\\left(a^{3} m - 2 \\, a m^{2} r_{0} + a m r_{0}^{2}\\right)} y + \\sqrt{-2 \\, a^{2} m r_{0} - 4 \\, m r_{0}^{3} + r_{0}^{4} + {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{2} + {\\left(a^{6} - 6 \\, a^{4} m r_{0} - 6 \\, m r_{0}^{5} + r_{0}^{6} + 3 \\, {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{4} - 4 \\, {\\left(3 \\, a^{2} m + 2 \\, m^{3}\\right)} r_{0}^{3} + 3 \\, {\\left(a^{4} + 4 \\, a^{2} m^{2}\\right)} r_{0}^{2}\\right)} y^{2}} r_{0}}{2 \\, a^{2} m + 4 \\, m r_{0}^{2} - r_{0}^{3} - {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}}\\right]\n",
"\\end{math}"
],
"text/plain": [
"[dt == (2*(a^3*m - 2*a*m^2*r0 + a*m*r0^2)*y - sqrt(-2*a^2*m*r0 - 4*m*r0^3 + r0^4 + (a^2 + 4*m^2)*r0^2 + (a^6 - 6*a^4*m*r0 - 6*m*r0^5 + r0^6 + 3*(a^2 + 4*m^2)*r0^4 - 4*(3*a^2*m + 2*m^3)*r0^3 + 3*(a^4 + 4*a^2*m^2)*r0^2)*y^2)*r0)/(2*a^2*m + 4*m*r0^2 - r0^3 - (a^2 + 4*m^2)*r0), dt == (2*(a^3*m - 2*a*m^2*r0 + a*m*r0^2)*y + sqrt(-2*a^2*m*r0 - 4*m*r0^3 + r0^4 + (a^2 + 4*m^2)*r0^2 + (a^6 - 6*a^4*m*r0 - 6*m*r0^5 + r0^6 + 3*(a^2 + 4*m^2)*r0^4 - 4*(3*a^2*m + 2*m^3)*r0^3 + 3*(a^4 + 4*a^2*m^2)*r0^2)*y^2)*r0)/(2*a^2*m + 4*m*r0^2 - r0^3 - (a^2 + 4*m^2)*r0)]"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dt, y, r0 = var('dt, y, r0')\n",
"\n",
"p = M((0, r0, pi/2, 0))\n",
"Tp = M.tangent_space(p)\n",
"v = Tp((dt, -1, 0, y))\n",
"\n",
"sol = g.at(p)(v, v).solve(dt)\n",
"sol"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def calc_some_geodesics(args):\n",
" \"\"\"\n",
" Compute nb geodesics starting at index n0\n",
" \"\"\"\n",
" curve, n0, nb = args\n",
" res = {}\n",
" r = 100\n",
" posi = [0, r, pi/2, 0]\n",
" p = M(posi)\n",
" Tp = M.tangent_space(p)\n",
" for i in range(n0, n0+nb):\n",
" dy = i*0.006/n_geod\n",
" v = Tp([sol[0].rhs()(r0=r, y=dy, m=m_val, a=a_val).n(), -1, 0, dy])\n",
" curve._initial_tangent_vector = v\n",
" curve.solve_across_charts(step=0.2, parameters_values={m:m_val, a:a_val})\n",
" res[i] = (p.coord(), curve._solutions.copy())\n",
" curve._solutions.clear()\n",
" return res"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fef6a25ffd99454da6f5a0bada40f10b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"FloatProgress(value=0.0, max=1000.0)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"geo = {}\n",
"pool = multiprocessing.Pool(n_cpu)\n",
"\n",
"%display plain\n",
"f = FloatProgress(min=0, max=n_geod)\n",
"display(f)\n",
"\n",
"for i, some_res in enumerate(pool.map(calc_some_geodesics, args)):\n",
" f.value += len(some_res)\n",
" geo.update(some_res)\n",
"\n",
"pool.close()\n",
"pool.join()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"\n",
"for i in range(0, n_geod, 5*n_geod/100): \n",
" curve._solutions = geo[i][1]\n",
" interp = curve.interpolate()\n",
" P += curve.plot_integrated(mapping=phi, color=[\"red\"], thickness=2, plot_points=150, \n",
" label_axes=False, across_charts=True)#\n",
"\n",
"P"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8.46600505906165\n"
]
}
],
"source": [
"Angle = pi/20 #angle used in the actual calculation\n",
"\n",
"Xi=a_val/m_val\n",
"Z_1=1+(1-Xi^2)^(1/3)*((1+Xi)^(1/3)+(1-Xi)^(1/3))\n",
"Z_2=(3*Xi^2+Z_1^2)^(1/2)\n",
"#calculating ISCO\n",
"\n",
"disk_min = m_val*(3+Z_2-((3-Z_1)*(3+Z_1+2*Z_2))^(1/2))\n",
"print(numerical_approx(disk_min))\n",
"if disk_min<2*m_val:\n",
" disk_min=2*m_val\n",
"\n",
"disk_max = (50^2+disk_min^2-12^2)^(1/2)\n",
"#disk_max =25\n",
"#keeping area of accretion Disk the same as in original work of Florentin Jaffredo\n",
"\n",
"#disk_min = 4\n",
"#disk_max = 50\n",
"alpha = -pi/20\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"