The basic imports and the variables we'll be using:
In [39]:
from __future__ import division, print_function
import sympy
from sympy import *
from sympy import Rational as frac
import simpletensors
from simpletensors import Vector, xHat, yHat, zHat
from simpletensors import TensorProduct, SymmetricTensorProduct, Tensor
init_printing()
var('vartheta, varphi')
var('nu, m, delta, c, t')
# These are related scalar functions of time
var('v, gamma, r', cls=Function)
v = v(t)
x = v**2
Omega = v**3
#gamma = v**2*(1 + (1-nu/3)*v**2 + (1-65*nu/12)*v**4) # APPROXIMATELY!!! Change this later
#r = 1/gamma
gamma = gamma(t)
r = r(t)
r = 1/v**2 # APPROXIMATELY!!! Just for testing. Change this later
# These get redefined momentarily, but have to exist first
var('nHat, lambdaHat, ellHat', cls=Function)
# And now we define them as vector functions of time
nHat = Vector('nHat', r'\hat{n}', [cos(Omega*t),sin(Omega*t),0,])(t)
lambdaHat = Vector('lambdaHat', r'\hat{\lambda}', [-sin(Omega*t),cos(Omega*t),0,])(t)
ellHat = Vector('ellHat', r'\hat{\ell}', [0,0,1,])(t)
# These are the spin functions -- first, the individual components as regular sympy.Function objects; then the vectors themselves
var('S_n, S_lambda, S_ell', cls=Function)
var('Sigma_n, Sigma_lambda, Sigma_ell', cls=Function)
SigmaVec = Vector('SigmaVec', r'\vec{\Sigma}', [Sigma_n(t), Sigma_lambda(t), Sigma_ell(t)])(t)
SVec = Vector('S', r'\vec{S}', [S_n(t), S_lambda(t), S_ell(t)])(t)
Thorne (1980) gives a nice review, along with the following formula for $Y^{\ell,m}_L$. Also note the Eqs. (74) and the footnote on page 32 of Blanchet's Living Review (2013), which explains the necessary normalizations for getting the metric perturbation modes from the tensors.
In [60]:
from __future__ import division
import sympy
from sympy import *
from sympy import Rational as frac
import simpletensors
from simpletensors import Vector, xHat, yHat, zHat
from simpletensors import TensorProduct, SymmetricTensorProduct, Tensor
init_printing()
var('vartheta, varphi')
DefaultOrthogonalRightHandedBasis=[xHat(t), yHat(t), zHat(t)]
def C(ell,m):
return (-1)**abs(m) * sympy.sqrt( frac(2*ell+1,4) * frac(factorial(ell-m), factorial(ell+m)) / sympy.pi )
def a(ell,m,j):
return frac((-1)**j, 2**ell * factorial(j) * factorial(ell-j)) * frac(factorial(2*ell-2*j), factorial(ell-m-2*j))
def YlmTensor(ell, m, OrthogonalRightHandedBasis=DefaultOrthogonalRightHandedBasis):
if(ell<0 or abs(m)>ell):
raise ValueError("YlmTensor({0},{1}) is undefined.".format(ell,m))
from sympy import prod
xHat, yHat, zHat = OrthogonalRightHandedBasis
if(m<0):
mVec = Tensor(SymmetricTensorProduct(xHat), SymmetricTensorProduct(yHat,coefficient=-sympy.I))
#mVec = VectorFactory('mBarVec', [1,-sympy.I,0])(t)
else:
mVec = Tensor(SymmetricTensorProduct(xHat), SymmetricTensorProduct(yHat,coefficient=sympy.I))
#mVec = VectorFactory('mVec', [1,sympy.I,0])(t)
def TensorPart(ell,m,j):
return sympy.prod((mVec,)*m) * SymmetricTensorProduct(*((zHat,)*(ell-2*j-m))) \
* sympy.prod([sum([SymmetricTensorProduct(vHat, vHat) for vHat in OrthogonalRightHandedBasis]) for i in range(j)])
if(m<0):
Y = sum([TensorPart(ell,-m,j) * (C(ell,-m) * a(ell,-m,j))
for j in range(floor(frac(ell+m,2))+1) ]) * (-1)**(-m)
else:
Y = sum([TensorPart(ell,m,j) * (C(ell,m) * a(ell,m,j))
for j in range(floor(frac(ell-m,2))+1) ])
try:
Y.compress()
except AttributeError:
pass
return Y
def YlmTensorConjugate(ell, m, OrthogonalRightHandedBasis=DefaultOrthogonalRightHandedBasis):
return YlmTensor(ell, -m, OrthogonalRightHandedBasis) * (-1)**abs(m)
# This is Blanchet's version of the above
def alphalmTensor(ell, m, OrthogonalRightHandedBasis=DefaultOrthogonalRightHandedBasis):
return YlmTensor(ell, -m, OrthogonalRightHandedBasis) * ( (-1)**abs(m) * (4*pi*factorial(ell)) / factorial2(2*ell+1) )
NVec = Vector('NVec', r'\hat{N}', [sympy.sin(vartheta)*sympy.cos(varphi),
sympy.sin(vartheta)*sympy.sin(varphi),
sympy.cos(vartheta)])(t)
def NTensor(ell):
return SymmetricTensorProduct(*((NVec,)*ell))
# These give the SWSH modes components from the given radiative tensors
def Ulm(U_L, m):
ell = U_L.rank
return (alphalmTensor(ell, m) | U_L) * (frac(4,factorial(ell))*sqrt(frac((ell+1)*(ell+2), 2*ell*(ell-1))))
def Vlm(V_L, m):
ell = V_L.rank
return (alphalmTensor(ell, m) | V_L) * (frac(-8,factorial(ell))*sqrt(frac(ell*(ell+2), 2*(ell+1)*(ell-1))))
def hlm(U_L, V_L, m):
return ( -Ulm(U_L,m) + sympy.I * Vlm(V_L,m) ) / sympy.sqrt(2)
Let's take a look at those tensors:
In [3]:
for ell in range(2,5):
print('')
for m in range(-ell, ell+1):
print('(ell,m) = ({0},{1}):'.format(ell,m))
display( YlmTensor(ell,m) )
In [4]:
NTensor(2)
Out[4]:
So now, I can just contract the $Y^{\ell,m}_L$ tensors with the $N_L$ tensors to get the usual spherical harmonics:
In [5]:
YlmTensor(0,0)
Out[5]:
In [6]:
for ell in range(1,3):
print('')
for m in range(-ell, ell+1):
print('(ell,m) = ({0},{1}):'.format(ell,m))
display(exptrigsimp( YlmTensor(ell,m) | NTensor(ell) ))
These values match, e.g., the ones on Wikipedia's spherical-harmonics page.
We can also see that the traces are zero. This is nice because it means we don't have to explicitly remove the traces from either these or the tensors with which they will be contracted.
In [7]:
AllTracesZero = True
for ell in range(2,9):
print('')
for m in range(-ell, ell+1):
print('(ell,m) = ({0},{1}):'.format(ell,m))
tmp = YlmTensor(ell,m).trace()
if(tmp!=0): AllTracesZero=False
display( tmp )
if(AllTracesZero):
print("Indeed, all traces were explicitly zero.")
else:
print("Not all traces were explicitly zero! Maybe they just failed to simplify...", file=sys.stderr)
We can just add in a couple low-order terms for $I_{jk}$ and $J_{jk}$, and see how they come out. I expect to find some new terms proportional to $\Sigma_n + i\, \Sigma_\lambda$ in the $(2,\pm2)$ modes, and recover the usual leading-order terms in the $(2,\pm1)$ modes. Note, however, that I am not conjugating the YlmTensor objects, as I should, so the normalization and signs will be screwy.
In [49]:
I_jk = SymmetricTensorProduct(nHat, nHat, coefficient = (m*nu*r**2))
J_jk = SymmetricTensorProduct(SigmaVec, nHat, coefficient = (r*nu/c)*(-frac(3,2))) \
+ SymmetricTensorProduct(ellHat, nHat, coefficient = -nu*m*delta*r**2*v)
In [66]:
U_jk = diff(I_jk,t,2).subs(t,0)
V_jk = diff(J_jk,t,2).subs(t,0)
display(U_jk)
display(V_jk)
In [68]:
Ulm(U_jk, 2)
Out[68]:
In [86]:
Ulm(U_jk, 2).subs([(diff(v,t,2).subs(t,0), 0),#(9216*nu**2/25)*v.subs(t,0)**17),
((diff(v,t).subs(t,0))**2,0),
(diff(v,t).subs(t,0),(32*nu/5)*v.subs(t,0)**9)])
Out[86]:
In [63]:
expand(_ * -1/sqrt(2))
Out[63]:
In [70]:
expand(_69 * (-1 / (sqrt(2)*(2*m*nu*v.subs(t,0)**2)*sqrt(16*pi/5)) ) )
Out[70]:
In [ ]:
In [ ]:
dvdt = dOmega**1/3dt = (1/3)*Omega**(-2/3)*dOmegadt = (1/3)*v**-2*dOmegadt
In [ ]:
Omega = v**3
In [ ]:
dOmegadt = (96*nu/5)*v**11
In [ ]:
dvdt = (32*nu/5)*v**9
In [ ]:
d2vdt2 = (288*nu/5)*v**8 * dvdt = (288*nu/5)*v**8 * (32*nu/5)*v**9 = (9216*nu**2/25)*v**17
In [56]:
32*288
Out[56]:
In [58]:
dvdt = (32*nu/5)*v**9
diff(dvdt,v)
Out[58]:
In [ ]:
In [ ]: