#!/usr/bin/env python
# coding: utf-8

# #  Lemaître-Tolman solutions
# 
# This Jupyter/SageMath notebook is relative to the lectures
# [Geometry and physics of black holes](https://relativite.obspm.fr/blackholes/).
# 
# The computations make use of tools developed through the [SageManifolds project](https://sagemanifolds.obspm.fr).

# In[1]:


version()


# First we set up the notebook to display mathematical objects using LaTeX rendering:

# In[2]:


get_ipython().run_line_magic('display', 'latex')


# ## Spacetime
# 
# We declare the spacetime manifold $M$:

# In[3]:


M = Manifold(4, 'M', structure='Lorentzian')
print(M)


# and declare the chart of **Lemaître synchronous coordinates** on it:

# In[4]:


X.<t,x,th,ph> = M.chart(r't:\tau x:(0,+oo):\chi th:(0,pi):\theta ph:(0,2*pi):\phi')
X


# The most general metric tensor, assuming spherical symmetry and synchronous coordinates:

# In[5]:


g = M.metric()
a = function('a')
r = function('r')
g[0,0] = -1
g[1,1] = a(t,x)^2
g[2,2] = r(t,x)^2
g[3,3] = (r(t,x)*sin(th))^2
g.display()


# ## Einstein equation

# The cosmological constant:

# In[6]:


Lamb = var('Lamb', latex_name=r'\Lambda')
Lamb


# The Ricci tensor:

# In[7]:


Ric = g.ricci()
print(Ric)


# In[8]:


Ric.display()


# The Einstein tensor:

# In[9]:


G = Ric - 1/2*g.ricci_scalar() * g
G.set_name('G')
print(G)


# In[10]:


G.display_comp()


# ### Dust matter model
# 
# Let us consider a pressureless fluid ("dust"). Moreover, we assume that the coordinates $(\tau,\chi,\theta,\phi)$ are **comoving**, i.e. that the fluid 4-velocity is equal to $\partial_\tau$:

# In[11]:


u = M.vector_field('u')
u[0] = 1
u.display()


# Since $(\tau,\chi,\theta,\chi)$ are synchronous, the above does define a unit timelike vector:

# In[12]:


g(u,u).display()


# Let us check that $u$ is a geodesic vector field:

# In[13]:


nabla = g.connection()
acc = u['^b']*nabla(u)['^a_b']
acc.display()


# The 1-form associated to the fluid 4-velocity by metric duality:

# In[14]:


u_form = u.down(g)
print(u_form)


# In[15]:


u_form.display()


# The pressureless energy-momentum tensor:

# In[16]:


rho = function('rho')
T = rho(t,x)*(u_form * u_form)
T.set_name('T')
print(T)


# In[17]:


T.display()


# The Einstein equation:

# In[18]:


E = G + Lamb*g - 8*pi*T 
E.set_name('E')
print(E)


# In[19]:


E.display_comp()


# ## Solving the Einstein equation

# ### $\tau\chi$ component

# Let us first consider the $01 = \tau\chi$ component of the Einstein equation:

# In[20]:


E[0,1]


# A slight rearrangement of the equation:

# In[21]:


eq = E[0,1]*r(t,x)/(-2*a(t,x)) 
eq


# We see that this equation is equivalent to 
# $$ \frac{\partial}{\partial\tau} \left( \frac{1}{a}\frac{\partial r}{\partial\chi} \right) = 0 $$
# since

# In[22]:


drdx = diff(r(t,x), x)
eq - diff(drdx/a(t,x), t)


# Hence there exists a function of $\chi$ only, $f(\chi)$ say, such that
# $\frac{1}{a}\frac{\partial r}{\partial\chi} = f(\chi)$.
# We disregard the case $f(\chi)=0$, which would imply $\frac{\partial r}{\partial\chi}=0$ 
# and would lead to the so-called *Datt model* (1938). Accordingly, we may write
# $$
#     a(\tau,\chi) = \frac{1}{f(\chi)}\frac{\partial r}{\partial\chi} 
# $$
# Let us call `af` this expression of $a$:

# In[23]:


f = function('f')
af(t,x) = drdx / f(x)
af(t,x)


# We check that if we substitute $a$ by `af` in the $\tau\chi$ component of the Einstein equation, we get identically zero:

# In[24]:


E[0,1].expr().substitute_function(a, af)


# *NB:* `expr()` returns a Sage symbolic expression from the coordinate function `E[0,1]`, so that we may apply `substitute_function`

# Hence the first Lemaitre-Tolman equation is

# In[25]:


LT1 = a(t,x) == af(t,x)
LT1


# ### $\chi\chi$ component
# 
# The $11 = \chi\chi$ component of Einstein equation is

# In[26]:


E[1,1]


# It is equivalent to

# In[27]:


eq = (- E[1,1] * r(t,x)^2).expr() == 0
eq


# Let us substitute for $a(\tau,\chi)$ the value found above when solving the $\tau\chi$ component:

# In[28]:


eq1 = eq.substitute_function(a, af)
eq1


# Some slight rearrangement and simplification:

# In[29]:


eq2 = (eq1 * f(x)^2 / diff(r(t,x), x)^2).simplify_full()
eq2


# In[30]:


eq3 = (eq2 * diff(r(t,x),t)).simplify_full()
eq3


# We notice that the left-hand side of this equation is nothing but the partial derivative w.r.t. $\tau$ of the following quantity:

# In[31]:


A = (diff(r(t,x),t)^2 + 1 - f(x)^2 - (Lamb/3)*r(t,x)^2) * r(t,x) 
A


# In[32]:


bool(eq3.lhs() == diff(A, t))


# Hence `eq3` tells that $A$ is independent of $\tau$, i.e. is a function of $\chi$ only, which we call $2 m(\chi)$:

# In[33]:


m = function('m')
eq4 = A - 2*m(x) == 0
eq4


# Let us solve extract $(\partial r/\partial\tau)^2$ from this equation:

# In[34]:


drdt2_sol = solve(eq4, diff(r(t,x),t)^2)
drdt2_sol


# We thus obtain the second Lemaitre-Tolman equation:

# In[35]:


LT2 = drdt2_sol[0].expand()
LT2


# In[36]:


drdt2 = LT2.rhs()
drdt2


# ### $\tau\tau$ component
# 
# The $00 = \tau\tau$ component of Einstein equation is

# In[37]:


E[0,0]


# It is equivalent to

# In[38]:


eq = (- E[0,0] * a(t,x)^3 * r(t,x)^2).expr() == 0
eq


# As above, we substitute for $a(\tau,\chi)$ the value found when solving the $\tau\chi$ component:

# In[39]:


eq1 = eq.substitute_function(a, af)
eq1


# In[40]:


eq2 = (eq1 * f(x)^3).simplify_full()
eq2


# Let us substitute for $\partial r/\partial \tau$ the positive square root of the value of $(\partial r/\partial \tau)^2$ found when solving the $\chi\chi$ component:

# In[41]:


drdt = sqrt(drdt2)
drdt


# In[42]:


eq3 = eq2.subs({diff(r(t,x),t): drdt, diff(r(t,x),t,x): diff(drdt, x)}).simplify_full()
eq3


# If we use the negative square root of $(\partial r/\partial \tau)^2$ instead, we get the same result:

# In[43]:


drdt = - sqrt(drdt2)
drdt


# In[44]:


eq3_minus = eq2.subs({diff(r(t,x),t): drdt, diff(r(t,x),t,x): diff(drdt, x)}).simplify_full()
eq3_minus


# In[45]:


eq3_minus == eq3


# Thus we continue with `eq3` and rearrange it to get the third Lemaitre-Tolman equation:

# In[46]:


eq4 = (eq3 / (2*diff(r(t,x),x)^2)).simplify_full()
eq4


# In[47]:


dmdx_sol = solve(eq4, diff(m(x),x))
dmdx_sol


# In[48]:


LT3 = dmdx_sol[0]
LT3


# ### $\theta\theta$ and $\phi\phi$ components
# 
# First we notice that the $\theta\theta$ and $\phi\phi$ components of the Einstein equation are equivalent:

# In[49]:


E[3,3] == E[2,2] * sin(th)^2


# Let us thus consider only the $22 = \theta\theta$ component:

# In[50]:


E[2,2]


# It is equivalent to

# In[51]:


eq = (- E[2,2] * a(t,x)^3).expr() == 0
eq


# We substitute for $a(\tau,\chi)$ the value found when solving the $\tau\chi$ component:

# In[52]:


eq1 = eq.substitute_function(a, af).simplify_full()
eq1


# In[53]:


eq2 = (eq1 * f(x)^3).simplify_full()
eq2


# Then we substitute for $\partial r/\partial\tau$ the value obtained when solving the $\tau\tau$ component:

# In[54]:


eq3 = eq2.subs({diff(r(t,x),t,t,x): diff(drdt,t,x), diff(r(t,x),t,t): diff(drdt,t)}).simplify_full()
eq4 = eq3.subs({diff(r(t,x),t): drdt, diff(r(t,x),t,x): diff(drdt,x)}).simplify_full()
eq4


# We conclude that the $\theta\theta$ component of Einstein equation does not add any independent equation. 

# ## Summary
# 
# Let us collect the three independent equations obtained from the Einstein equation, constituting the **Lemaître-Tolman system**:

# In[55]:


for eq in [LT1, LT2, LT3]:
    show(eq)


# The first equation is the unnumbered one just above Eq. (8.1) in Lemaître's article *L'univers en expansion*, Annales de la Société Scientifique de Bruxelles A **53**, 51 (1933), translated in English in  [Gen. Relativ. Gravit. **29**, 641](http://dx.doi.org/10.1023/A:1018855621348) (1997). The second equation is Eq. (8.2) in Lemaître's article, while the third one is Eq. (8.3).