The basic imports and the variables we'll be using:


```python
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)
```

# Defining the spherical-harmonic tensors

[Thorne (1980)](http://link.aps.org/doi/10.1103/RevModPhys.52.299) 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)](http://arxiv.org/abs/1310.1528), which explains the necessary normalizations for getting the metric perturbation modes from the tensors.


```python
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:


```python
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) )
```


```python
NTensor(2)
```

# Spherical-harmonic functions

So now, I can just contract the $Y^{\ell,m}_L$ tensors with the $N_L$ tensors to get the usual spherical harmonics:


```python
YlmTensor(0,0)
```


```python
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](https://en.wikipedia.org/wiki/Spherical_harmonics#List_of_spherical_harmonics).

## Traces of the Spherical-harmonic tensors

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.


```python
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)
```

# Testing contractions with the multipole moments

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.


```python
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)
```


```python
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)
```


```python
Ulm(U_jk, 2)
```


```python
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)])
```


```python
expand(_ * -1/sqrt(2))
```


```python
expand(_69 * (-1 / (sqrt(2)*(2*m*nu*v.subs(t,0)**2)*sqrt(16*pi/5)) ) )
```


```python

```


```python
dvdt = dOmega**1/3dt = (1/3)*Omega**(-2/3)*dOmegadt = (1/3)*v**-2*dOmegadt
```


```python
Omega = v**3
```


```python
dOmegadt = (96*nu/5)*v**11
```


```python
dvdt = (32*nu/5)*v**9
```


```python
d2vdt2 = (288*nu/5)*v**8 * dvdt = (288*nu/5)*v**8 * (32*nu/5)*v**9 = (9216*nu**2/25)*v**17
```


```python
32*288
```


```python
dvdt = (32*nu/5)*v**9
diff(dvdt,v)
```


```python

```


```python

```

7fb0a8ba4d2da5ea4c00a051356fc7bd4da74b56
ipynb
Jupyter Notebook
Waveforms/SphericalHarmonicTensors.ipynb
keefemitman/PostNewtonian
853d6577cb0002da5eebe1cb55f0c28fbc114324
MIT
2015-03-26T01:04:36.000Z
2022-02-01T19:26:21.000Z
Waveforms/SphericalHarmonicTensors.ipynb
keefemitman/PostNewtonian
853d6577cb0002da5eebe1cb55f0c28fbc114324
MIT
2015-01-08T23:46:29.000Z
2017-09-20T19:13:51.000Z
Waveforms/SphericalHarmonicTensors.ipynb
keefemitman/PostNewtonian
853d6577cb0002da5eebe1cb55f0c28fbc114324
MIT
2016-05-13T02:36:14.000Z
2021-11-23T21:36:32.000Z
Qwen/Qwen-72B
1. YES
2. YES