{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Study of stable circular orbits in the equatorial plane\n", "\n", "This Jupyter/SageMath notebook is related to the article [Lamy et al, arXiv:1802.01635](https://arxiv.org/abs/1802.01635)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'SageMath version 8.1, Release Date: 2017-12-07'" ] }, "execution_count": 1, "metadata": { }, "output_type": "execute_result" } ], "source": [ "version()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "%display latex" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "from sage.manifolds.utilities import simplify_chain_real as simpl" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Spacetime" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "a0 = 0.7\n", "b0 = 1" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\theta}$ " ] }, "execution_count": 5, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('r a b')\n", "var('th', latex_name=r'\\theta')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}$ " ] }, "execution_count": 6, "metadata": { }, "output_type": "execute_result" } ], "source": [ "m(b,r) = (1-unit_step(r))* r^3/(r^3 - 2*b^2) + \\\n", " unit_step(r)* r^3/(r^3 + 2*b^2)\n", "m(b,r)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}1$ " ] }, "execution_count": 7, "metadata": { }, "output_type": "execute_result" } ], "source": [ "m(0,r) # Kerr limit" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "22e21cbf60931b001ab417e3b921ce78a3494004" }, "execution_count": 8, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot(m(1, r), (r, -6, 6), axes_labels=[r\"$r$\", r\"$m(r)$\"], \n", " title=r\"$b=1$\", title_pos=(0.9, 1.1))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "### Derivatives of the mass function" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{3 \\, r^{5} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} - \\frac{3 \\, r^{5} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} + \\frac{r^{3} \\delta\\left(r\\right)}{r^{3} + 2 \\, b^{2}} - \\frac{r^{3} \\delta\\left(r\\right)}{r^{3} - 2 \\, b^{2}} - \\frac{3 \\, r^{2} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{3 \\, r^{2} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}$ " ] }, "execution_count": 9, "metadata": { }, "output_type": "execute_result" } ], "source": [ "md(b, r) = diff(m(b, r), r)\n", "md(b, r)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Because of the step function, the derivative contains Dirac delta functions. We get rid of them by means of the following Python function:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "def remove_dirac(f):\n", " return f.subs({dirac_delta(r): 0}).simplify()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{3 \\, r^{5} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} - \\frac{3 \\, r^{5} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} - \\frac{3 \\, r^{2} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{3 \\, r^{2} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}$ " ] }, "execution_count": 11, "metadata": { }, "output_type": "execute_result" } ], "source": [ "md(b, r) = remove_dirac(md(b, r))\n", "md(b, r)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "bc8961c390dfd2791a231218e30453c08aa07f97" }, "execution_count": 12, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot(md(1, r), (r, -6, 6), axes_labels=[r\"$r$\", r\"$m'(r)$\"],\n", " title=r\"$b=1$\", title_pos=(0.9, 1))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{18 \\, r^{7} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{3}} + \\frac{18 \\, r^{7} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{3}} + \\frac{24 \\, r^{4} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} - \\frac{24 \\, r^{4} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} - \\frac{6 \\, r {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{6 \\, r \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}$ " ] }, "execution_count": 13, "metadata": { }, "output_type": "execute_result" } ], "source": [ "mdd(b, r) = remove_dirac(diff(md(b, r), r))\n", "mdd(b, r)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "623557bcfae37d342a147261d7efae89c22e0277" }, "execution_count": 14, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot(mdd(1, r), (r, -6, 6), axes_labels=[r\"$r$\", r\"$m''\\,(r)$\"],\n", " title=r\"$b=1$\", title_pos=(0.9, 1))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "### Horizons" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "8e3c72cd248c14e83ab9a58bdda4a7f413e7310e" }, "execution_count": 15, "metadata": { }, "output_type": "execute_result" } ], "source": [ "D=r^2+a^2-2*m(b,r)*r\n", "plot(D.subs({a:a0},{b:b0}),(r,-5,5))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "#find_root(D.subs({a:a0},{b:b0}),0,1)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Circular orbits in the equatorial plane (Condition1>0 & Condition2>0)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## $E$, $L$ and $\\Omega$ of a test particle in the equatorial plane" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "#t_dot(a,b,E,L,r,th) = ((r^2+a^2)*(E*(r^2+a^2)-a*L)/D - a*(a*E*sin(th)^2-L))/Sigma\n", "#phi_dot(a,b,E,L,r,th) = (a*(E*(r^2+a^2)-a*L)/D - (a*E-L/sin(th)^2))/Sigma\n", "#mu0(a,b,E,L,r,th) = -(g00*t_dot^2 + 2*g03*t_dot*phi_dot + g33*phi_dot^2) #normalization (m0=0 for a photon)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(m_{0}, b\\right)$ " ] }, "execution_count": 18, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('m0','b')" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "M = function('M')(r)\n", "#M=r^3/(r^3+2*b^2)\n", "Sigma = r^2+a^2*cos(th)^2\n", "Delta = r^2+a^2-2*M*r\n", "g00 = -(1-2*r*M/Sigma)\n", "g03 = -2*a*r*sin(th)^2*M/Sigma\n", "g33 = (r^2 + a^2 + 2*r*(a*sin(th))^2*M/Sigma)*sin(th)^2" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{3 \\, r^{5}}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} + \\frac{3 \\, r^{2}}{r^{3} + 2 \\, b^{2}}$ " ] }, "execution_count": 20, "metadata": { }, "output_type": "execute_result" } ], "source": [ "diff(r^3/(r^3+2*b^2),r)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "var('e l q');" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "dg00=diff(g00,r)\n", "dg03=diff(g03,r)\n", "dg33=diff(g33,r)\n", "Mdot=diff(M,r)\n", "omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 \n", "omega_neg=(-dg03-sqrt(dg03^2-dg00*dg33))/dg33" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{a r \\frac{\\partial}{\\partial r}M\\left(r\\right) + \\sqrt{-r \\frac{\\partial}{\\partial r}M\\left(r\\right) + M\\left(r\\right)} r^{\\frac{3}{2}} - a M\\left(r\\right)}{a^{2} r \\frac{\\partial}{\\partial r}M\\left(r\\right) + r^{3} - a^{2} M\\left(r\\right)}$ " ] }, "execution_count": 23, "metadata": { }, "output_type": "execute_result" } ], "source": [ "simpl(omega_pos.subs({b:0, th:pi/2}))" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "E_pos=-(g00+g03*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2)) \n", "E_neg=-(g00+g03*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "L_pos=(g03+g33*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2))\n", "L_neg=(g03+g33*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{a^{2} r \\frac{\\partial}{\\partial r}M\\left(r\\right) + r^{3} + 2 \\, \\sqrt{-r \\frac{\\partial}{\\partial r}M\\left(r\\right) + M\\left(r\\right)} a \\sqrt{r} M\\left(r\\right) - {\\left(a^{2} + 2 \\, r^{2}\\right)} M\\left(r\\right)}{\\sqrt{-a^{2} r^{3} \\frac{\\partial}{\\partial r}M\\left(r\\right)^{2} + r^{5} - 3 \\, a^{2} r M\\left(r\\right)^{2} + 2 \\, {\\left({\\left(a^{3} + 3 \\, a r^{2}\\right)} M\\left(r\\right) - {\\left(a^{3} r + a r^{3}\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right)\\right)} \\sqrt{-r \\frac{\\partial}{\\partial r}M\\left(r\\right) + M\\left(r\\right)} \\sqrt{r} - 3 \\, {\\left(a^{2} r^{2} + r^{4}\\right)} M\\left(r\\right) + {\\left(3 \\, a^{2} r^{3} + r^{5} + 4 \\, a^{2} r^{2} M\\left(r\\right)\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right)} \\sqrt{r}}$ " ] }, "execution_count": 26, "metadata": { }, "output_type": "execute_result" } ], "source": [ "simpl(E_pos.subs({th:pi/2,b:0})).simplify_full()" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{{\\left(a^{3} + 3 \\, a r^{2}\\right)} M\\left(r\\right) - {\\left(a^{3} r + a r^{3}\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right) - {\\left(a^{2} r^{2} + r^{4} + 2 \\, a^{2} r M\\left(r\\right)\\right)} \\sqrt{-\\frac{r \\frac{\\partial}{\\partial r}M\\left(r\\right) - M\\left(r\\right)}{r}}}{{\\left(a^{2} r \\frac{\\partial}{\\partial r}M\\left(r\\right) + r^{3} - a^{2} M\\left(r\\right)\\right)} \\sqrt{-\\frac{a^{2} r^{4} \\frac{\\partial}{\\partial r}M\\left(r\\right)^{2} - r^{6} + 3 \\, a^{2} r^{2} M\\left(r\\right)^{2} + 3 \\, {\\left(a^{2} r^{3} + r^{5}\\right)} M\\left(r\\right) - {\\left(3 \\, a^{2} r^{4} + r^{6} + 4 \\, a^{2} r^{3} M\\left(r\\right)\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right) - 2 \\, {\\left({\\left(a^{3} r^{2} + 3 \\, a r^{4}\\right)} M\\left(r\\right) - {\\left(a^{3} r^{3} + a r^{5}\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right)\\right)} \\sqrt{-\\frac{r \\frac{\\partial}{\\partial r}M\\left(r\\right) - M\\left(r\\right)}{r}}}{a^{4} r^{2} \\frac{\\partial}{\\partial r}M\\left(r\\right)^{2} + r^{6} - 2 \\, a^{2} r^{3} M\\left(r\\right) + a^{4} M\\left(r\\right)^{2} + 2 \\, {\\left(a^{2} r^{4} - a^{4} r M\\left(r\\right)\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right)}}}$ " ] }, "execution_count": 27, "metadata": { }, "output_type": "execute_result" } ], "source": [ "L_pos.subs({th:pi/2,b:0}).simplify_full()" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{{\\left(a^{3} + 3 \\, a r^{2}\\right)} M\\left(r\\right) - {\\left(a^{3} r + a r^{3}\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right) + {\\left(a^{2} r^{2} + r^{4} + 2 \\, a^{2} r M\\left(r\\right)\\right)} \\sqrt{-\\frac{r \\frac{\\partial}{\\partial r}M\\left(r\\right) - M\\left(r\\right)}{r}}}{{\\left(a^{2} r \\frac{\\partial}{\\partial r}M\\left(r\\right) + r^{3} - a^{2} M\\left(r\\right)\\right)} \\sqrt{-\\frac{a^{2} r^{4} \\frac{\\partial}{\\partial r}M\\left(r\\right)^{2} - r^{6} + 3 \\, a^{2} r^{2} M\\left(r\\right)^{2} + 3 \\, {\\left(a^{2} r^{3} + r^{5}\\right)} M\\left(r\\right) - {\\left(3 \\, a^{2} r^{4} + r^{6} + 4 \\, a^{2} r^{3} M\\left(r\\right)\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right) + 2 \\, {\\left({\\left(a^{3} r^{2} + 3 \\, a r^{4}\\right)} M\\left(r\\right) - {\\left(a^{3} r^{3} + a r^{5}\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right)\\right)} \\sqrt{-\\frac{r \\frac{\\partial}{\\partial r}M\\left(r\\right) - M\\left(r\\right)}{r}}}{a^{4} r^{2} \\frac{\\partial}{\\partial r}M\\left(r\\right)^{2} + r^{6} - 2 \\, a^{2} r^{3} M\\left(r\\right) + a^{4} M\\left(r\\right)^{2} + 2 \\, {\\left(a^{2} r^{4} - a^{4} r M\\left(r\\right)\\right)} \\frac{\\partial}{\\partial r}M\\left(r\\right)}}}$ " ] }, "execution_count": 28, "metadata": { }, "output_type": "execute_result" } ], "source": [ "L_neg.subs({th:pi/2,b:0}).simplify_full()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Conditions 1 & 2" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "reset('M')" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "M=m(b,r)\n", "Md = md(b ,r)\n", "Mdd = mdd(b, r)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{3 \\, r^{5} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} + \\frac{3 \\, r^{5} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} + \\frac{3 \\, r^{2} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{3 \\, r^{2} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}} - \\frac{\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}}{r}$ " ] }, "execution_count": 31, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition1 = M/r - Md\n", "Condition1" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{1}{r}$ " ] }, "execution_count": 32, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition1.subs({b:0}) # Kerr limit" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "Condition2_pos=-(a^2*r^2*Md^2 - r^4 + 3*a^2*M^2 + 3*(a^2*r + r^3)*M - (3*a^2*r^2 + r^4 + 4*a^2*r*M)*Md - 2*((a^3 + 3*a*r^2)*M - (a^3*r + a*r^3)*Md*sqrt(-(Md - M/r))))" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "Condition2_neg=-(a^2*r^2*Md^2 - r^4 + 3*a^2*M^2 + 3*(a^2*r + r^3)*M - (3*a^2*r^2 + r^4 + 4*a^2*r*M)*Md + 2*((a^3 + 3*a*r^2)*M - (a^3*r + a*r^3)*Md*sqrt(-(Md - M/r))))" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r^{4} + 2 \\, a^{3} - 3 \\, a^{2} r + 6 \\, a r^{2} - 3 \\, r^{3} - 3 \\, a^{2}$ " ] }, "execution_count": 35, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition2_pos.subs({b:0}) # Kerr limit" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r^{4} - 3 \\, r^{3}$ " ] }, "execution_count": 36, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition2_pos.subs({a:0},{b:0}) # Schwarzschild limit" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "p1 = plot(20, -1.587,0, fill=True, fillalpha=0.3)\n", "p2 = plot(20, 1.587, 5, fill=True, fillalpha=0.3)\n", "p3 = plot(20, 1.31, 1.34, fill=True, fillcolor='black', fillalpha=1)\n", "p4 = plot(20, 0.89, 0.92, fill=True, fillcolor='black', fillalpha=1)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 76 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 75 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "data": { "image/png": "59d7106437053c07c991eb6bb82d07c5c4e53c8d" }, "execution_count": 38, "metadata": { }, "output_type": "execute_result" } ], "source": [ "p = plot((Condition1.subs({b:b0}),\n", " Condition2_pos.subs({a:a0},{b:b0}),\n", " Condition2_neg.subs({a:a0},{b:b0})), (r,-3,5), ymax=12, ymin=-1.5,\n", " legend_label=['$A(r)$','$B_+(r)$','$B_-(r)$'], color=('blue','red','green'), \n", " axes_labels=[r'$r/m$',''])\n", "graph = p+p1+p2#+p3+p4\n", "show(graph)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "graph.save('Circular_orbits_a=07_b=1.pdf')" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "#find_root(Condition2_pos.subs({a:a0},{b:b0}),-1,-0.3)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}8.52502568341426 - 0.0409208183571186i$ " ] }, "execution_count": 41, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition2_pos.subs({a:a0},{b:b0},{r:1.587})" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Stable circular orbits in the equatorial plane (Condition1>0, Condition2>0 & Condition 3<0), ISCO" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Defining $R$" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(E, L, Q, m_{0}\\right)$ " ] }, "execution_count": 42, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('E L Q m0')" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "reset('M')\n", "M=m(b,r)\n", "#Md = md(b ,r)\n", "#Mdd = mdd(b, r)\n", "Sigma = r^2+a^2*cos(th)^2\n", "Delta = r^2+a^2-2*M*r\n", "g00 = -(1-2*r*M/Sigma)\n", "g03 = -2*a*r*sin(th)^2*M/Sigma\n", "g33 = (r^2 + a^2 + 2*r*(a*sin(th))^2*M/Sigma)*sin(th)^2\n", "dg00=diff(g00,r)\n", "dg03=diff(g03,r)\n", "dg33=diff(g33,r)\n", "omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 \n", "omega_neg=(-dg03-sqrt(dg03^2-dg00*dg33))/dg33\n", "E_pos=-(g00+g03*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2)) \n", "E_neg=-(g00+g03*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))\n", "L_pos=(g03+g33*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2))\n", "L_neg=(g03+g33*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\left({\\left(a^{2} + r^{2}\\right)} E - L a\\right)}^{2} - {\\left(m_{0}^{2} r^{2} + {\\left(E a - L\\right)}^{2} + Q\\right)} {\\left(a^{2} + 2 \\, {\\left(\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}\\right)} r + r^{2}\\right)}$ " ] }, "execution_count": 44, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R = ((r^2+a^2)*E - a*L)^2 - Delta*((a*E-L)^2 + m0^2*r^2 + Q)\n", "R" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}2 \\, m_{0}^{2} r^{3} + {\\left(E^{2} - m_{0}^{2}\\right)} r^{4} - Q a^{2} + {\\left(E^{2} a^{2} - a^{2} m_{0}^{2} - L^{2} - Q\\right)} r^{2} + 2 \\, {\\left(E^{2} a^{2} - 2 \\, E L a + L^{2} + Q\\right)} r$ " ] }, "execution_count": 45, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R.subs({b:0}).simplify_full() # Kerr limit" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}2 \\, m_{0}^{2} r^{3} + {\\left(E^{2} - m_{0}^{2}\\right)} r^{4} - {\\left(L^{2} + Q\\right)} r^{2} + 2 \\, {\\left(L^{2} + Q\\right)} r$ " ] }, "execution_count": 46, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R.subs({a:0},{b:0}).simplify_full() # Schwarzschild limit" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-2 \\, {\\left(a^{2} + 2 \\, {\\left(\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}\\right)} r + r^{2}\\right)} m_{0}^{2} r + 4 \\, {\\left({\\left(a^{2} + r^{2}\\right)} E - L a\\right)} E r - 2 \\, {\\left(m_{0}^{2} r^{2} + {\\left(E a - L\\right)}^{2} + Q\\right)} {\\left(\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}} - 3 \\, {\\left(\\frac{r^{5} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} - \\frac{r^{5} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} - \\frac{r^{2} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{r^{2} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}\\right)} r + r\\right)}$ " ] }, "execution_count": 47, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rd = remove_dirac(diff(R, r))\n", "Rd" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}2 \\, E^{2} a^{2} + 6 \\, m_{0}^{2} r^{2} + 4 \\, {\\left(E^{2} - m_{0}^{2}\\right)} r^{3} - 4 \\, E L a + 2 \\, L^{2} + 2 \\, {\\left(E^{2} a^{2} - a^{2} m_{0}^{2} - L^{2} - Q\\right)} r + 2 \\, Q$ " ] }, "execution_count": 48, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rd.subs({b:0}).simplify_full() # Kerr limit" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}6 \\, m_{0}^{2} r^{2} + 4 \\, {\\left(E^{2} - m_{0}^{2}\\right)} r^{3} + 2 \\, L^{2} - 2 \\, {\\left(L^{2} + Q\\right)} r + 2 \\, Q$ " ] }, "execution_count": 49, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rd.subs({a:0},{b:0}).simplify_full() # Schwarzschild limit" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-8 \\, {\\left(\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}} - 3 \\, {\\left(\\frac{r^{5} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} - \\frac{r^{5} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} - \\frac{r^{2} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{r^{2} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}\\right)} r + r\\right)} m_{0}^{2} r + 8 \\, E^{2} r^{2} - 2 \\, {\\left(a^{2} + 2 \\, {\\left(\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}\\right)} r + r^{2}\\right)} m_{0}^{2} + 2 \\, {\\left(m_{0}^{2} r^{2} + {\\left(E a - L\\right)}^{2} + Q\\right)} {\\left(\\frac{6 \\, r^{5} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} - \\frac{6 \\, r^{5} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} - 6 \\, {\\left(\\frac{3 \\, r^{7} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{3}} - \\frac{3 \\, r^{7} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{3}} - \\frac{4 \\, r^{4} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{{\\left(r^{3} - 2 \\, b^{2}\\right)}^{2}} + \\frac{4 \\, r^{4} \\mathrm{u}\\left(r\\right)}{{\\left(r^{3} + 2 \\, b^{2}\\right)}^{2}} + \\frac{r {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} - \\frac{r \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}}\\right)} r - \\frac{6 \\, r^{2} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2 \\, b^{2}} + \\frac{6 \\, r^{2} \\mathrm{u}\\left(r\\right)}{r^{3} + 2 \\, b^{2}} - 1\\right)} + 4 \\, {\\left({\\left(a^{2} + r^{2}\\right)} E - L a\\right)} E$ " ] }, "execution_count": 50, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rdd = remove_dirac(diff(Rd, r))\n", "Rdd" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}2 \\, E^{2} a^{2} - 2 \\, a^{2} m_{0}^{2} + 12 \\, m_{0}^{2} r + 12 \\, {\\left(E^{2} - m_{0}^{2}\\right)} r^{2} - 2 \\, L^{2} - 2 \\, Q$ " ] }, "execution_count": 51, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rdd.subs({b:0}).simplify_full() # Kerr limit" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}12 \\, m_{0}^{2} r + 12 \\, {\\left(E^{2} - m_{0}^{2}\\right)} r^{2} - 2 \\, L^{2} - 2 \\, Q$ " ] }, "execution_count": 52, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rdd.subs({a:0},{b:0}).simplify_full() # Schwarzschild limit" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "Rd2=Rdd.subs({E:e},{L:l},{Q:q})\n", "Rd2_normalized=Rd2.subs({m0:1}).simplify() #normalization by m0" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}2 \\, a^{2} e^{2} + 12 \\, {\\left(e^{2} - 1\\right)} r^{2} - 2 \\, a^{2} - 2 \\, l^{2} - 2 \\, q + 12 \\, r$ " ] }, "execution_count": 54, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rd2_normalized.subs({b:0}).simplify_full() # Kerr limit" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}12 \\, {\\left(e^{2} - 1\\right)} r^{2} - 2 \\, l^{2} - 2 \\, q + 12 \\, r$ " ] }, "execution_count": 55, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Rd2_normalized.subs({a:0},{b:0}).simplify_full() # Schwarzschild limit" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Conditions 1, 2 & 3" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "q1=0 #value of q=Q/m0^2 corresponding to a motion in the equatorial plane \n", "e1_pos=E_pos.subs({a:a0},{b:b0}) #value of e=E/m0 corresponding to a circular orbit\n", "e1_neg=E_neg.subs({a:a0},{b:b0})\n", "l1_pos=L_pos.subs({a:a0},{b:b0}) #value of l=L/m0 corresponding to a circular orbit\n", "l1_neg=L_neg.subs({a:a0},{b:b0})" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "Condition3_pos = Rd2_normalized.subs({e:e1_pos},{l:l1_pos},{q:q1}).simplify()\n", "Condition3_neg = Rd2_normalized.subs({e:e1_neg},{l:l1_neg},{q:q1}).simplify() " ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}NaN$ " ] }, "execution_count": 58, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition3_neg.subs({a:a0},{b:b0},{th:pi/2},{r:1})" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 67 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: 'unable to convert -488.7212108566773 + 1.1209327348532101e-13*I to float; use abs() or real_part() as desired'\n" ] }, { "data": { "image/png": "7ddbfbdd85e91f60169c4561b659ee35ba679d56" }, "execution_count": 59, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot(Condition3_neg.subs({a:a0},{b:b0},{th:pi/2}),(r,0,10))\n", "#pc3.save('ISCO_a0_b0.pdf')" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 47 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 72 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "data": { "image/png": "23179e68c3408106f33b912c6070d81f442190c5" }, "execution_count": 60, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot((Condition1.subs({b:b0}), \n", " Condition2_neg.subs({a:a0},{b:b0}),\n", " Condition3_neg.subs({a:a0},{b:b0},{th:pi/2})), (r,-3,10), ymax=10, ymin=-10,\n", " legend_label=['A(r)','B(r)','C(r)'], color=('blue','red','green'), \n", " title='Counter-rotating', axes_labels=[r'$r$',''])" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 46 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 57 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 46 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 73 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''\n" ] }, { "data": { "image/png": "508721ba993913ad76dac432dd25f77aad19ef49" }, "execution_count": 61, "metadata": { }, "output_type": "execute_result" }, { "data": { "image/png": "154b5abb1e7810fa23fe5cda84334857d935f2e2" }, "execution_count": 61, "metadata": { }, "output_type": "execute_result" } ], "source": [ "p_pos = plot((Condition1.subs({b:b0}), \n", " Condition2_pos.subs({a:a0},{b:b0}),\n", " Condition3_pos.subs({a:a0},{b:b0},{th:pi/2})), (r,-3,10), ymax=10, ymin=-10,\n", " legend_label=['A(r)','B(r)','C(r)'], color=('blue','red','green'), \n", " title='Co-rotating',axes_labels=[r'$r$',''])\n", "p_neg = plot((Condition1.subs({b:b0}), \n", " Condition2_neg.subs({a:a0},{b:b0}),\n", " Condition3_neg.subs({a:a0},{b:b0},{th:pi/2})), (r,-3,10), ymax=10, ymin=-10,\n", " legend_label=['A(r)','B(r)','C(r)'], color=('blue','red','green'), \n", " title='Counter-rotating', axes_labels=[r'$r$',''])\n", "show(p_pos)\n", "show(p_neg)" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}0.0$ " ] }, "execution_count": 62, "metadata": { }, "output_type": "execute_result" } ], "source": [ "find_root(Condition2_pos.subs({a:a0},{b:b0},{th:pi/2}),1.001,2)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 157 points.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "verbose 0 (3749: plot.py, generate_plot_points) Last error message: 'unable to convert -8.32736020802459 + 2.683660004211231*I to float; use abs() or real_part() as desired'\n" ] }, { "data": { "image/png": "c25851821cb8abff2c1cbbc5a9eebeafaf6654a3" }, "execution_count": 63, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot(Condition3_pos.subs({a:a0},{b:b0},{th:pi/2}),(r,0,2))" ] }, { "cell_type": "code", "execution_count": 64, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}1.62683109949525 + 1.36936138033619i$ " ] }, "execution_count": 64, "metadata": { }, "output_type": "execute_result" } ], "source": [ "Condition3_pos.subs({a:a0},{b:b0},{th:pi/2},{r:1.04})" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Frequency of the orbits at the ISCO" ] }, { "cell_type": "code", "execution_count": 65, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "%display latex" ] }, { "cell_type": "code", "execution_count": 66, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\theta}$ " ] }, "execution_count": 66, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('r a b M0 m')\n", "var('th', latex_name=r'\\theta')" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "%display latex\n", "var('r a b M0 m')\n", "var('th', latex_name=r'\\theta')\n", "\n", "m0=4.31*10^6*2*10^30\n", "M(b,r) = m0*r^3/(r^3 + 2*m0*b^2)\n", "#M=function('M',r)\n", "Sigma = r^2 + a^2*cos(th)^2\n", "Delta = r^2 - 2*M*r + a^2\n", "g00(a,b,r) = -(1 - 2*r*M/Sigma)\n", "g03(a,b,r) = -2*a*r*sin(th)^2*M/Sigma\n", "g33(a,b,r) = (r^2 + a^2 + 2*r*(a*sin(th))^2*M/Sigma)*sin(th)^2" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "dg00=diff(g00,r)\n", "dg03=diff(g03,r)\n", "dg33=diff(g33,r)\n", "omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 \n", "omega_neg=(-dg03-sqrt(dg03^2-dg00*dg33))/dg33" ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "a1=m0*0.9\n", "b1=m0*1\n", "r1_pos=m0*1.59\n", "r1_neg=m0*8.43" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( a, b, r \\right) \\ {\\mapsto} \\ 0.00107187558536188$ " ] }, "execution_count": 70, "metadata": { }, "output_type": "execute_result" } ], "source": [ "om_pos=omega_pos.subs({a:a1},{b:b1},{r:r1_pos},{th:pi/2})\n", "(3*10^8)^3/(6.67*10^(-11))*om_pos" ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( a, b, r \\right) \\ {\\mapsto} \\ -0.00197809424213708$ " ] }, "execution_count": 71, "metadata": { }, "output_type": "execute_result" } ], "source": [ "om_neg=omega_neg.subs({a:a1},{b:b1},{r:r1_neg},{th:pi/2})\n", "(3*10^8)^3/(6.67*10^(-11))*om_neg" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# General case: only $R(r,E,L,Q) \\geq 0$ allowed" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "reset('m0')" ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}m_{0}$ " ] }, "execution_count": 73, "metadata": { }, "output_type": "execute_result" } ], "source": [ "var('m0')" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\left(E^{2} - 1\\right)} r^{4} - Q a^{2} + {\\left({\\left(E^{2} - 1\\right)} a^{2} - L^{2} - Q\\right)} r^{2} + 2 \\, r^{3} + 2 \\, {\\left(E^{2} a^{2} - 2 \\, E L a + L^{2} + Q\\right)} r$ " ] }, "execution_count": 74, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R.subs({b:0},{m0:1}).simplify_full() # Kerr limit" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}E^{2} r^{4} - {\\left(L^{2} + Q\\right)} r^{2} + 2 \\, {\\left(L^{2} + Q\\right)} r$ " ] }, "execution_count": 75, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R.subs({a:0},{b:0},{m0:0}).simplify_full() # Schwarzschild limit" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "630ead0275296434a1a2f7b8df3f48541993a1bd" }, "execution_count": 76, "metadata": { }, "output_type": "execute_result" } ], "source": [ "a2 = 0.9\n", "b2 = 1\n", "e2 = 1\n", "l2 = 2\n", "q2 = -1\n", "plot0=plot(R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2}), (r,-2,2), axes_labels = [r'$r/m$', r'$R(r,E_1,L_1,Q_1)$'])\n", "show(plot0)" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\left(r^{2} - 0.990000000000000\\right)}^{2} - 0.420000000000000 \\, {\\left(\\frac{r^{3} {\\left(\\mathrm{u}\\left(r\\right) - 1\\right)}}{r^{3} - 2} - \\frac{r^{3} \\mathrm{u}\\left(r\\right)}{r^{3} + 2}\\right)} r - 0.210000000000000 \\, r^{2} - 0.170100000000000$ " ] }, "execution_count": 77, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2})" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "#find_root(R_null(a2,b2,e2,l2,q2,r),-2,2)" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-0.240000000000000$ " ] }, "execution_count": 79, "metadata": { }, "output_type": "execute_result" } ], "source": [ "R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2},{r:1})" ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "003292f45f0e362068aa437788f6ee062c4c6aac" }, "execution_count": 80, "metadata": { }, "output_type": "execute_result" } ], "source": [ "#p1 = plot(11, -0.665666,0.713731,ymax=8, fill=True, fillalpha=0.3)\n", "p2 = plot(25, -0.054222, 2, ymax=22, ymin=-10, fill=True, fillalpha=0.3)\n", "p3 = plot(11,-2, -1.446841, fill=True, fillalpha=0.3)\n", "p4 = plot(25, 1.425, 1.44, fill=True, fillcolor='black', fillalpha=1)\n", "p5 = plot(25, 0.545, 0.56, fill=True, fillcolor='black', fillalpha=1)\n", "plot2=plot0+p2+p4+p5\n", "show(plot2)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "#plot2.save('R_photons_dark_zone_b=0.pdf')" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath (stable)", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.15" } }, "nbformat": 4, "nbformat_minor": 0 }