Energy and angular momentum fluxes

In this example we will compute the fluxes of energy and angular momentum for a point mass on a circular orbit in Schwarzschild spacetime.

Schwarzschild spacetime

We want to consider point-particle perturbations of a Schwarzschild spacetime. Let’s first load the Schwarzschild spacetime using GeneralRelativityTensors

<<GeneralRelativityTensors`

g=ToMetric["Schwarzschild"];
MatrixForm[TensorValues[g]]

We will later also need the square root of the negative of the determinant of the metric

SqrtDetg=Simplify[(Sqrt[-Det[TensorValues[g]]]),Assumptions->r>0&&0≤θ≤π]

Spin raising and lowering operators

The spin-raising and lowering operators

ð=(m^μ)(∂_μ)-(b+s)β+(b-s)(OverBar[β])'
ð'=((OverBar[m])^μ)(∂_μ)+(b+s)β'-(b-s)(OverBar[β])

take us from one spin-weighted spherical harmonic to another. Here, β and β'=-α are spin coefficients and an overbar denotes complex conjugation. If we work with the Kinnersley tetrad in Schwarzschild spacetime, this gives the identity

ð (Subscript[, s])(Subscript[Y, ℓ𝓂])(θ,ϕ)=-(1/((Sqrt[2])r))(Sqrt[ℓ(ℓ+1)-s(s+1)])(Subscript[, s+1])(Subscript[Y, ℓ𝓂])(θ,ϕ)
ð' (Subscript[, s])(Subscript[Y, ℓ𝓂])(θ,ϕ)=(1/((Sqrt[2])r))(Sqrt[ℓ(ℓ+1)-s(s-1)])(Subscript[, s-1])(Subscript[Y, ℓ𝓂])(θ,ϕ)

Let’s check this explicitly using the SpinWeightedSphericalHarmonicY function provided by the SpinWeightedSpheroidalHarmonics package. First, we define our differential operators acting on a scalar of spin-weight s and boost-weight b.

ð[scalar_,{s_,b_}]:=KinnersleyDerivative["Schwarzschild","delta"][scalar]-((b+s) SpinCoefficient["Schwarzschild","beta"]+(b-s) SpinCoefficient["Schwarzschild","alpha",True])scalar;
ðp[scalar_,{s_,b_}]:=KinnersleyDerivative["Schwarzschild","deltaStar"][scalar]-((b+s) SpinCoefficient["Schwarzschild","alpha"]+(b-s) SpinCoefficient["Schwarzschild","beta",True])scalar;

For an arbitrary function f(θ,ϕ) these operators in the Schwarzschild case are given explicitly by

Simplify[ð[f[θ,ϕ],{s,b}]]

Simplify[ðp[f[θ,ϕ],{s,b}]]

Now, let’s check the above identities for a range of values of s, ℓ, 𝓂 in the spin-weighted spherical harmonics (for arbitrary θ, ϕ):

<<SpinWeightedSpheroidalHarmonics`

Table[Simplify[ð[SpinWeightedSphericalHarmonicY[s,l,m,θ,ϕ],{s,0}]==(1/((Sqrt[2])r))(-(Sqrt[l(l+1)-s(s+1)]) )SpinWeightedSphericalHarmonicY[s+1,l,m,θ,ϕ]],{s,-2,2},{l,Abs[s],4},{m,-l,l}]

Table[Simplify[ðp[SpinWeightedSphericalHarmonicY[s,l,m,θ,ϕ],{s,0}]==(1/((Sqrt[2])r))((Sqrt[l(l+1)-s(s-1)]) )SpinWeightedSphericalHarmonicY[s-1,l,m,θ,ϕ]],{s,-2,2},{l,Abs[s],4},{m,-l,l}]

Boost raising and lowering operators

Similar to spin-raising and lowering operators, the boost-raising and lowering operators

þ=(l^μ)(∂_μ)-(b+s)ϵ-(b-s)(OverBar[ϵ])
þ'=(n^μ)(∂_μ)+(b+s)ϵ'+(b-s)(OverBar[ϵ])'

take us from one boost-weighted function to another. Here, ϵ and ϵ'=-γ are spin coefficients and an overbar denotes complex conjugation. Let’s now define our differential operators acting on a scalar of spin-weight s and boost-weight b.

þ[scalar_,{s_,b_}]:=KinnersleyDerivative["Schwarzschild","D"][scalar]-((b+s) SpinCoefficient["Schwarzschild","epsilon"]+(b-s) SpinCoefficient["Schwarzschild","epsilon",True])scalar;
þp[scalar_,{s_,b_}]:=KinnersleyDerivative["Schwarzschild","Delta"][scalar]-((b+s) SpinCoefficient["Schwarzschild","gamma"]+(b-s) SpinCoefficient["Schwarzschild","gamma",True])scalar;

If we work with the Kinnersley tetrad in Schwarzschild spacetime, for an arbitrary function f(t,r) these operators are given explicitly by

Collect[þ[f[t,r],{s,b}],{Derivative[__][f][t,r],f[t,r]},Simplify]

Collect[þp[f[t,r],{s,b}],{Derivative[__][f][t,r],f[t,r]},Simplify]

Point particle stress-energy

Tensor components of the stress-energy tensor

The stress-energy tensor of a point particle is given by 

(T^(μν))=-(1/(Sqrt[-g]))(((u^μ)(u^ν))/(u^t))δ[r-(Subscript[r, 0])(t)]δ[θ-(Subscript[θ, 0])(t)]δ[ϕ-(Subscript[ϕ, 0])(t)]

We can use the completeness relation of the spin-weighted spherical harmonics,

(Subscript[Σ, ℓ𝓂])(Subscript[[, s])(Subscript[(OverBar[Y]), ℓ𝓂])(θ',ϕ')(Subscript[, s])(Subscript[Y, ℓ𝓂])(θ,ϕ)]=δ[cosθ'-cosθ]δ[ϕ'-ϕ]

to rewrite this as

(T^(μν))=(Subscript[Σ, ℓ𝓂])(-(1/(r^2))(((u^μ)(u^ν))/(u^t))δ[r-(Subscript[r, 0])(t)](Subscript[, s])(Subscript[(OverBar[Y]), ℓ𝓂])((Subscript[θ, 0]),(Subscript[ϕ, 0])))(Subscript[, s])(Subscript[Y, ℓ𝓂])(θ,ϕ)

Let’s consider a single ℓ,𝓂 mode (i.e. a single term in the sum) but note that we have not yet specified the spin-weight s (that will come later). The four-velocity can be written in terms of the constants of motion, which for the case of a point particle on a circular orbit with radius (Subscript[r, 0]) are given by

<<KerrGeodesics`

M=1;
$Assumptions=r0>3;

Let’s now define the four-velocity vector, setting its covariant components

u=ToTensor["u",g,{-"ℰ",0,0,"ℒ"}/.Simplify[KerrGeoConstantsOfMotion[0,r0,0,1]],{-α}];
u[-α]==TensorValues[u[-α]]

We can then define the stress-energy tensor

T=ToTensor["T",g,Simplify[-δ[r-r0]((1/((r^2)u[t]))TensorValues[MergeTensors[u[-α]u[-β]]](OverBar[Y])[π/2,Ω t]/.r->r0)], {-α,-β}];
TensorValues[T]//MatrixForm

Kinnersley tetrad projection of the stress-energy tensor

The source for the s=-2 Teukolsky equation is given in terms of the tetrad components of the stress-energy tensor. We thus start by projecting the stress-energy tensor onto the Kinnersley tetrad

{l,n,m,(OverBar[m])}=KinnersleyNullTetrad["Schwarzschild"];

Tnn=Simplify[TensorValues[MergeTensors[T[-μ,-ν]n[μ]n[ν]]]/.x_ δ[r-r0]:>(x/.r->r0)δ[r-r0]]

Tnmb=Simplify[TensorValues[MergeTensors[T[-μ,-ν]n[μ](OverBar[m])[ν]]]/.x_ δ[r-r0]:>(x/.r->r0)δ[r-r0]/.θ->π/2]

Tmbmb=Simplify[TensorValues[MergeTensors[T[-μ,-ν](OverBar[m])[μ](OverBar[m])[ν]]]/.x_ δ[r-r0]:>(x/.r->r0)δ[r-r0]/.θ->π/2]

Now, we use the fact that the tetrad components have spin and boost weights as follows:

(Subscript[T, nn]): s=0, b=-2
(Subscript[T, n(OverBar[m])]): s=-1, b=-1
(Subscript[T, (OverBar[m])(OverBar[m])]): s=-2, b=0

We will later choose a harmonic with spin-weight s appropriate for the component.

Source for the Teukolsky equation

The source for the s=-2 Teukolsky equation is then given by

𝒯=(ð'-(OverBar[τ])-4τ')((þ'-2(OverBar[ρ])')(Subscript[T, n(OverBar[m])])-(ð'-(OverBar[τ])')(Subscript[T, nn]))+(þ'-4ρ'-(OverBar[ρ])')((ð'-2(OverBar[τ]))(Subscript[T, n(OverBar[m])])-(þ'-(OverBar[ρ])')(Subscript[T, (OverBar[m])(OverBar[m])]))

where τ, τ'=-ϖ, ρ and ρ'=-μ are spin coefficients and (Subscript[T, nn]), (Subscript[T, n(OverBar[m])]) and (Subscript[T, (OverBar[m])(OverBar[m])]) are tetrad components of the stress-energy tensor.

Now, acting on these with the appropriate differential operators above and making use of the spin-lowering property of ð' along with the fact that τ=0=τ' for Schwarzschild we get three terms, all of spin-weight s=-2

SpinCoefficient["Schwarzschild","tau"]

SpinCoefficient["Schwarzschild","pi"]

𝒯0=-(1/((Sqrt[2])r))(1/((Sqrt[2])r))(Sqrt[With[{s=0},Factor[ℓ(ℓ+1)-s(s-1)]]With[{s=-1},Factor[ℓ(ℓ+1)-s(s-1)]]])Simplify[(Tnn/. Derivative[0,1][(OverBar[Y])][(π/2),t Ω]->-I 𝓂 (OverBar[Y])[(π/2),t Ω] /.(OverBar[Y])[(π/2),t Ω]->(OverBar[Y])[0,(π/2),t Ω])]

ρp=-SpinCoefficient["Schwarzschild","mu"];
ρpb=-SpinCoefficient["Schwarzschild","mu",True];

𝒯1=(Sqrt[With[{s=-1},Factor[ℓ(ℓ+1)-s(s-1)]]])Simplify[(((1/((Sqrt[2])r))(þp[#,{-1,-1}]-2ρpb #)&[Tnmb])+((þp[#,{-2,-1}]-4ρp#-ρpb #)&[(1/((Sqrt[2])r))Tnmb])/. Derivative[0,1][(OverBar[Y])][(π/2),t Ω]->-I 𝓂 (OverBar[Y])[(π/2),t Ω] /.(OverBar[Y])[(π/2),t Ω]->(OverBar[Y])[-1,(π/2),t Ω])]

𝒯2=-Simplify[((þp[#,{-2,-1}]-4ρp#-ρpb #)&[(þp[#,{-2,0}]-ρpb #)&[Tmbmb]]/.{Derivative[0,1][(OverBar[Y])][(π/2),t Ω]->-I 𝓂 (OverBar[Y])[(π/2),t Ω],Derivative[0,2][(OverBar[Y])][(π/2),t Ω]->(((-I 𝓂))^2) (OverBar[Y])[-2,(π/2),t Ω] }/.(OverBar[Y])[(π/2),t Ω]->(OverBar[Y])[-2,(π/2),t Ω])]

Variation of parameters weighting coefficients

The inhomogeneous solutions to the Teukolsky equation can be computed using variation of parameters and integrating homogeneous solutions against the source,

ψ(r)=(C^in)(r) (R^in)(r)+(C^up)(r)(R^up)(r)

where

(C^in)(r) = (∫_r^(Subscript[r, ℐ]))(((R^up)(r'))/(W(r')Δ))𝒯(r')dr' and (C^up)(r) = (∫_(Subscript[r, ℋ])^r)(((R^in)(r'))/(W(r')Δ))𝒯(r')dr'.

We can do these integrals explicitly in this case as we have delta-distributions. Leaving out the constant invariant Wronskian (Δ^(s+1))W(r)=const where W(r)=(R^in)(∂_r)(R^up)-(R^up)(∂_r)(R^in) we get contributions from each of the three terms.

Δ=r(r-2);
Σ=(r^2);

C0=Collect[Integrate[-4(Δ^(-2))R[r](r^4)𝒯0 Σ/.δ->DiracDelta,{r,2,∞},Assumptions->2<r0<∞],{(OverBar[Y])[_,(π/2),t Ω],Derivative[_][R][_],R[_]},Simplify]

C1=Collect[Integrate[-4(Δ^(-2))R[r](r^4)𝒯1 (r^2)/.δ->DiracDelta,{r,2,∞},Assumptions->2<r0<∞],{(OverBar[Y])[_,(π/2),t Ω],Derivative[_][R][_],R[_]},Simplify]

C2=Collect[Integrate[-4(Δ^(-2))R[r](r^4)𝒯2 (r^2)/.δ->DiracDelta,{r,2,∞},Assumptions->2<r0<∞],{(OverBar[Y])[_,(π/2),t Ω],Derivative[_][R][_],R[_]},Simplify]/.(Ω^2)->M/(r0^3)

Numerical calculation

We will now use the results from the previous sections to numerically compute the energy flux

Inhomogenous mode amplitudes

We can use the Teukolsky package to numerically compute the homogeneous solutions to the Teukolsky equation

<<Teukolsky`

Let’s consider the (2,2) mode of (Subscript[ψ, 4]) for a circular orbit of radius (Subscript[r, 0])=10 about a Schwarzschild black hole. We use KerrGeodesics to compute the orbital frequency

𝓈=-2;ℓ=2;𝓂=2;

a=0;
r0=10.0`32;
Ω=KerrGeoFrequencies[0,r0,0,1]["(Subscript[Ω, ϕ])"];
ω0=𝓂 Ω;

Δ0=r0(r0-2);

Now compute the Teukolsky homogeneous solutions for this case

{RIn,RUp}={"In","Up"}/.TeukolskyRadial[𝓈,ℓ,𝓂,a,ω0]

Next, compute the invariant Wronskian

Winv=(RIn[r0]RUp'[r0]-RUp[r0]RIn'[r0])((Δ0)^(𝓈+1))

And finally compute the weighting coefficients

CUp=((C0+C1+C2)/Winv)/.{(OverBar[Y])[s_,(π/2),t Ω]:>SpinWeightedSphericalHarmonicY[s,ℓ,𝓂,π/2,0],R->RIn}

CIn=((C0+C1+C2)/Winv)/.{(OverBar[Y])[s_,(π/2),t Ω]:>SpinWeightedSphericalHarmonicY[s,ℓ,𝓂,π/2,0],R->RUp}

In fact, the Teukolsky package has all of this built-in, so we can directly compute it using TeukolskyPointParticleMode

ψ4=TeukolskyPointParticleMode[𝓈,ℓ,𝓂,0,0,KerrGeoOrbit[a,r0,0,1]]

As expected, our computed values for the weighting coefficients agree perfectly with the asymptotic amplitudes

((Values[ψ4["Amplitudes"]])/(2π{CUp,CIn}))

Fluxes

We the asymptotic amplitudes at hand, we can now compute the flux of energy through future null infinity

ℱℐ=Abs[((((2π CUp))^2)/(4π ((ω0)^2)))]

and through the horizon

αℓ𝓂ω=With[{rh=2M},((256(((2M rh))^5)(((ω0)^2)+4 ((((1/(4rh))))^2))(((ω0)^2)+16 ((((1/(4rh))))^2))((ω0)^4))/((Abs[(ℓ-1)ℓ(ℓ+1)(ℓ+2)+(((-1))^(ℓ+𝓂))12 I M ω0])^2))];
ℱℋ=αℓ𝓂ω Abs[((((2π CIn))^2)/(4π ((ω0)^2)))]

Again, these are also built in to TeukolskyMode, which produces the same results

ψ4["EnergyFlux"]

The package can also compute the angular momentum flux

ψ4["AngularMomentumFlux"]

Compare with Regge-Wheeler calculation

We can compare the previous calculation to the same thing done within the Regge-Wheeler formalism (we can ignore the warning which is just telling us that there are two packages providing the same RenormalizedAngularMomentum function)

<<ReggeWheeler`

ψRW=ReggeWheelerPointParticleMode[𝓈,ℓ,𝓂,0,KerrGeoOrbit[a,r0,0,1]]

ψRW["EnergyFlux"]

Compare with post-Newtonian approximation

Finally, we can also compare against a post-Newtonian approximation using the series expansions included in the PostNewtonianSelfForce package

<<PostNewtonianSelfForce`

ℱℐPN=PostNewtonianExpansion[{"Schwarzschild","Circular","Flux","Infinity-l2m2","!Spin"}];
ℱℐPN["Series"]+((O[y])^8)

Block[{$MaxExtraPrecision=1000},((Normal[ℱℐPN["Series"]])/2)/.y->1/10.`64]

Exercise

[Note: Before starting this exercise, make sure to clear previously defined values.]

M=.;
r0=.;
Ω=.;

The stress-energy for a point particle with spin to linear order in the spin magnitude σ is given by

(T^(μν))=(T^((pp)μν))+σ (T^((σ)μν))

where (T^((pp)μν)) is the previous point-particle stress-energy tensor (but with an extra contribution to the four-velocity coming from the spin) and

(T^((σ)μν))=-(1/((u^t)(Sqrt[-g])))(+(∂_ρ){((OverTilde[S])^(ρ(μ))(u^(ν)))δ[r-(Subscript[r, 0])(t)]δ[θ-(Subscript[θ, 0])(t)]δ[ϕ-(Subscript[ϕ, 0])(t)]}+((OverTilde[S])^(ρ(μ))(Subscript[Γ, ρσ]^(ν)))(u^σ)δ[r-(Subscript[r, 0])(t)]δ[θ-(Subscript[θ, 0])(t)]δ[ϕ-(Subscript[ϕ, 0])(t)])
where a tilde on a tensor denotes that it has been divided by the spin magnitude σ

The spin vector for the aligned-spin, circular orbit case only has one non-zero component, ((OverBar[S])^θ)=-M/(Subscript[r, 0]). This is related to the spin tensor by ((OverTilde[S])^(μν))=-(ϵ^(μ ν α β))(Subscript[(OverTilde[S]), α])(Subscript[u, β]). Let’s now compute the components of ((OverTilde[S])^(μν)). First we define the Levi-Civita tensor

ϵT=ToTensor["ϵ",g,SqrtDetg Normal[LeviCivitaTensor[4]],{-μ,-α,-β,-γ}]

The four-velocity will be given in terms of arbitrary covariant components (Subscript[u, t]) and (Subscript[u, ϕ]) for now

u=ToTensor["u",g,{ut,0,0,Ω ut},{α}];

These are given explicitly to linear order in spin by

(u^t)=(1/(Sqrt[1-3M/(Subscript[r, 0])]))-σ ((3(M^(5/2)))/(2(((Subscript[r, 0])((Subscript[r, 0])-3M))^(3/2))))
Ω=(Sqrt[(M/(Subscript[r, 0]^3))])-σ ((3(M^2))/(2(Subscript[r, 0]^3)))

We define a spin vector with the components given above

Svec=ToTensor["(OverTilde[S])",g,{0,0,-M/r0,0},{μ}];

S=ToTensor["(OverTilde[S])",g,Simplify[TensorValues[MergeTensors[-ϵT[μ,ν,α,β]Svec[-α]u[-β]]]/.{θ->π/2,r->r0}],{μ,ν}];
TensorValues[S[μ,ν]]//MatrixForm

Now we compute (Sqrt[-g])(T^((σ)μν))

Γ=ChristoffelSymbol[g];

δTens=ToTensor["δ",g,δ[r]δ[θ]δ[ϕ]];

derivTerm=Simplify[TensorValues[MergeTensors[((-SqrtDetg)/(SqrtDetg u[t]))(1/2)(S[ρ,μ]u[ν]+S[ρ,ν]u[μ])D[δTens,-ρ]]]/.{δ[r]->δ[r],δ[θ]->δ[θ],δ'[ϕ]->-I m δ[ϕ],δ[ϕ]->δ[ϕ],θ->π/2}]

ΓTerm=Factor[TensorValues[MergeTensors[-(SqrtDetg/(SqrtDetg u[t]))(1/2)(S[ρ,μ]Γ[ν,-ρ,-σ]+S[ρ,ν]Γ[μ,-ρ,-σ])u[σ]δTens]]/.{δ[r]->δ[r],δ[θ]->δ[θ],δ'[ϕ]->-I m δ[ϕ],δ[ϕ]->δ[ϕ],θ->π/2,r->r0}]

Now that we have our two contributions, we define a tensor with those components

Tσ=ToTensor["(T^((σ)))",g,Collect[derivTerm+ΓTerm,{δ[θ],δ[ϕ],δ'[r]},Simplify],{μ,ν}];

Coefficient[TensorValues[Tσ[μ,ν]],δ[r]]/.{ut->(Sqrt[(r0/(r0-3))]),Ω->(r0^(-3/2))}//Simplify//Factor//MatrixForm

Coefficient[TensorValues[Tσ[μ,ν]],δ'[r]]/.{ut->(Sqrt[(r0/(r0-3))]),Ω->(r0^(-3/2))}//Simplify//Factor//MatrixForm

Next, we proceed exactly as we did in the point-particle case. Project onto the Kinnersley tetrad, act on the tetrad components with differential operators, then derive an expression for the weighting coefficients.