{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cosmological Constant Correction\n",
"The cosmological constant is the simplest model for dark energy, which provides a mechanism for the accelerating expansion of the universe. In the vicinity of a spherically symmetric point mass, $M$, the gravitational potential is approximately\n",
"\n",
"$$\\Phi = -\\frac{GM}{r} - \\frac{c^2 \\Lambda r^2}{6},$$\n",
"\n",
"where $G$ is Newton's gravitational constant and $\\Lambda$ is the cosmological constant. Consequently, objects straying too far from the central mass $M$ can be blown away by cosmic inflation if $\\Lambda > 0$.\n",
"\n",
"### Problem\n",
"1. Find the effective potential, $V_{eff}$, for a particle orbiting the mass $M$, explaining why angular momentum is conserved.\n",
"2. Sketch the effective potential for a range of angular momenta and cosmological constants, identifying where gravity becomes repulsive.\n",
"3. Determine a condition on this crossover radius. Test this by using `CosmologicalConstant.ipynb` to visualise the paths around $M$ for various starting positions, velocities and values of $\\Lambda$.\n",
"4. By considering the effective potential, explain the relationship between an orbit's energy and its precession.\n",
"5. Using the above crossover radius as an order of magnitude estimate for the attractive-repulsive crossover radius, discuss the implications for the solar system if the cosmological constant has the value $6 \\times 10^{-34}$m$^{-2}$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#NAME: Cosmological Constant\n",
"#DESCRIPTION: Investigating the cosmological constant correction term for the gravitational potential around a point mass.\n",
"\n",
"import numpy as np\n",
"from scipy.integrate import odeint\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import animation\n",
"#HTML is used to run the animation inline\n",
"from IPython.display import HTML\n",
"\n",
"#gravitational constant\n",
"G = 1.0\n",
"#mass of star\n",
"M = 100.0\n",
"#cosmological constant * c^2\n",
"C = 0.1\n",
"\n",
"#starting co-ordinates (x, y, v_x, v_y)\n",
"a0 = [5.0, 0.0, -3.0, 3.0]\n",
"\n",
"#time steps\n",
"dt = 0.05\n",
"Nsteps = 500\n",
"t = np.linspace(0, Nsteps*dt, Nsteps)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equation of motion for a particle in gravitational field $\\vec{g}$ is\n",
"$$\\frac{d^2 \\vec{r}}{dt^2} = \\vec{g}.$$\n",
"This can be separated into two sets of first-order differential equations which, setting \n",
"$$\\vec{a} = \\left( \\begin{matrix}\n",
" x\\\\\n",
" y\\\\\n",
" \\frac{dx}{dt}\\\\\n",
" \\frac{dy}{dt}\\\\\n",
" \\end{matrix} \\right)$$\n",
"can be written in 2D as\n",
"$$\\frac{d \\vec{a}}{dt} = \\left( \\begin{matrix}\n",
" a_2\\\\\n",
" a_3\\\\\n",
" g_0\\\\\n",
" g_1\\\\\n",
" \\end{matrix} \\right)$$"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#calculate gravitational field\n",
"def g(x,y):\n",
" #x and y are scalar floats\n",
" r = np.sqrt(x**2 + y**2)\n",
" x_component = (-G*M/(r**3) + C/3)*x\n",
" y_component = (-G*M/(r**3) + C/3)*y\n",
" return x_component, y_component\n",
"\n",
"#derivates of coordinates\n",
"def derivative(a, t):\n",
" gx, gy = g(a[0], a[1])\n",
" return [a[2], a[3], gx, gy]\n",
"\n",
"#solve the differential equation above\n",
"trajectory = odeint(derivative, a0, t)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#create the animation\n",
"fig = plt.figure(figsize = (8,6))\n",
"plt.xlim(-10,10)\n",
"plt.ylim(-10,10)\n",
"plt.title(\"Cosmological Constant Correction \\n to Newton's Law of Gravitation\")\n",
"plt.xlabel('x')\n",
"plt.ylabel('y')\n",
"#central mass\n",
"mass, = plt.plot([0.0],[0.0], 'ko')\n",
"#orbiting body\n",
"point, = plt.plot([trajectory[0][0]],[trajectory[0][1]], 'bo')\n",
"path, = plt.plot([trajectory[0][0]],[trajectory[0][1]], 'b-')\n",
"\n",
"def animator(i):\n",
" point.set_data([trajectory[i][0]], [trajectory[i][1]])\n",
" path.set_data(trajectory[:i,0], trajectory[:i,1])\n",
" return point, path\n",
"\n",
"anim = animation.FuncAnimation(fig, animator, Nsteps, interval = 20, repeat = False, blit = True)\n",
"html_video = anim.to_html5_video()\n",
"HTML(html_video)"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "Python 3"
},
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}