Fundamental constants

Physical constants used throughout the book

(*First, we clear all previous definitions to ensure a clean start. *)

ClearAll["Global`*"];

cSI=Quantity[299792458,"Meters"/"Seconds"]; (*Speed of light*)

cCGS=UnitConvert[cSI,"Centimeters"/"Seconds"];

ergPerJoule=Quantity[10^7,"Ergs"/"Joules"];

(*Planck constant and related quantities*)
hSI=Quantity[6.62606957*10^-34,"Joule"*"Seconds"];

hbarSI=hSI/(2*Pi);

hCGS=UnitConvert[hSI,"Ergs"*"Seconds"];

ElectronVoltsPerJoule=Quantity[6.24150934*10^18,"Electronvolts"/"Joule"];

heV=hSI*ElectronVoltsPerJoule;

Quantity[4.135667510864479`*^-15,"Electronvolts" "Seconds"]

ElectronVoltsPerJoule hSI

(*Boltzmann constant*)
kBSI=Quantity[1.3806488*10^-23,"Joules"/"Kelvins"];

kBCGS=UnitConvert[kBSI,"Ergs"/"Kelvins"];

kBeVperK=UnitConvert[kBSI,"Electronvolts"/"Kelvins"];

WienDisplacementConstantHzperK=Quantity[5.878925757646824946*10^10,"Hertz"/"Kelvins"];

(*Newtonian gravitational constant*)
GSI=Quantity[6.673848*10^-11,"Meters"^3/"Kilograms"/"Seconds"^2];

mPerKg = GSI/cSI^2; (* Conversion to geometrized units *)

(* Vacuum permittivity (epsilon0SI) this name has changed*)
epsilon0SI = 1/(4*Pi*Quantity[8.988*10^9, "Newtons"*"Meters"^2/"Coulombs"^2];);

(*Radiation density constant*)
RadiationDensityConstantSI=(8*Pi^5*kBSI^4)/(15*hSI^3*cSI^3)//UnitSimplify;

RadiationDensityConstantCGS=(8*Pi^5*kBCGS^4)/(15*hCGS^3*cCGS^3)//UnitSimplify;

(*Cosmological constants*)
secPerDay=Quantity[86400,"Seconds"/"Days"];

secPerYear=Quantity[365.2425,"Days"/"Years"]*secPerDay;

1.0006643835616438`

parsecSI=Quantity[3.08567758149137*10^16,"Meters"];

lightyearSI=Quantity[365.2425,"Days"]*secPerDay*cSI;

(*Hubble parameter*)
h=0.7; (*Example dimensionless Hubble parameter*)

hubbleSI=Quantity[(100*h*10^3)/(3.08567758149137*10^22),"Seconds"^-1];

Quantity[2.2685455026110534`*^-18,(1/("Seconds"))]

(*Electron*)
ChargeElectronSI=Quantity[1.602*10^-19,"Coulombs"];

ElectronMassSI=Quantity[9.109*10^-31,"Kilograms"];

(*Proton and Neutron*)
MassProtonSI=Quantity[1.672621778*10^-27,"Kilograms"];

MassNeutronSI=Quantity[1.67492735*10^-27,"Kilograms"];

MassProtoneV=UnitConvert[MassProtonSI*cSI^2,"Electronvolts"];

MassProtonGU=MassProtonSI*GSI/cSI^2;

(*Earth*)
earthMassSI=Quantity[5.9736*10^24,"Kilograms"];

earthRadiusSI=Quantity[6.371*10^6,"Meters"];

earthEquatorialRadiusSI=Quantity[6378137,"Meters"];

earthOmegaSI=2*Pi/(Quantity[23*3600+56*60+4.091,"Seconds"]);

earthEquatorSpeedSI=earthEquatorialRadiusSI*earthOmegaSI;

(*Moon*)
MoonDistSI=Quantity[384399000,"Meters"];

MoonMassSI=Quantity[7.3477*10^22,"Kilograms"];

(*Sun*)
MassSunSI=Quantity[1.98855*10^30,"Kilograms"];

RadiusSunSI=Quantity[6.958*10^8,"Meters"];

SunLuminositySI=Quantity[3.828*10^26,"Watts"];

Chapter 1: Special relativity

Exercise 1.1(a,c,e,g)

Convert the following to units in which c=1, expressing everything in terms of [m] and [kg]

(*1.1(a) 10[J] in natural units*)
AnswerCh1Ex1a=N[10/cSI^2,3]

(*1.1(c) hbar [J s] in natural units*)

AnswerCh1Ex1c=N[hbarSI/cSI]

(*1.1(e) 30 x 10^3[kg m/s] in natural units*)

AnswerCh1Ex1e=N[30000/cSI]

(*1.1(g) Density of water, 1000[kg/m^3] in natural units*)

Here we do not need to convert the units because there are no seconds to convert to meters! Later, in Chapter 8, we will see how to convert kilograms to meters.

Exercise 1.2 (a,c,e):

Convert the following from natural units (c=1) to SI units:

(*1.2(a) v=10^-2*)
AnswerCh1Ex2a=0.01*cSI;

0.01` cSI

(*1.2(c) Time, t=10^18*)
AnswerCh1Ex2c=N[Quantity[10^18,"Meters"]/cSI,4];

Quantity[3.3356409519815204957557671447`4.*^9,"Seconds"]

(*1.2(e) Acceleration a=10 m^-1*)

AnswerCh1Ex2e=N[Quantity[10,"Meters"^-1]*cSI^2,4];

Quantity[8.9875517873681764`4.*^17,(("Meters")/(("Seconds")^2))]

Exercise 1.3 (c,g,h,i):

(* Exer. 1.3(a) The worldline of a particle, e.g. a clock, stationary in this coordinate system. *)Graphics[{Black,Thick,Line[{{1,-1},{1,3}}]},PlotRange->{{-1,3},{-1,3}},Axes->True,AxesLabel->{"x","t"}]

(* Exer. 1.3(b) The worldline of a particle moving at a tenth of the speed of light along the x axis *)
Plot[10*(x-0.5),{x,-1,3},PlotRange->{-1,3},PlotStyle->Black,AxesLabel->{"x","t"}]

(* Exer 1.3(c): The axes of an observer moving at half the speed of light along the x axis *)(* xbarLabel=OverBar["x"];
tbarLabel=OverBar["t"]; *)
xbarLabel=OverBar["x"];
tbarLabel=OverBar["t"];

AnswerEx1o3c=Plot[{2*x,1/2*x},{x,-1,4},PlotRange->{{-1,3.5},{-1,4}},PlotStyle->Directive[Black,AbsoluteThickness[0.5],Solid],AxesLabel->{"x","t"},(*Keep regular axes labels*)PlotLabels->{Placed[tbarLabel,Top],Placed[xbarLabel,Right]},(*These would now label the curves*)LabelStyle->{FontFamily->"Times",FontSize->18,FontSlant->"Italic"},AspectRatio->1]

(OverBar["x"])

AnswerEx1o3cbis=Plot[{2*x,1/2*x},{x,-1,4},PlotRange->{{-1,4},{-1,4}},(*Adjusted to show the full graph and labels*)PlotStyle->Directive[Black,AbsoluteThickness[0.5],Solid],AxesLabel->{"x","t"},
BaseStyle->Directive[Black,FontFamily->"Times",FontSlant->"Italic",FontSize->18],
TicksStyle->Directive[Black,FontFamily->"Helvetica",FontSlant->"Plain",FontSize->16],(*Applies to the numbers on the ticks (0,1,2,...,6) and (-1.0,-0.5,...)*)
Epilog->{Text[Style[xbarLabel,FontFamily->"Times",FontSlant->"Italic",FontSize->18],{3.6,1.5}](*Absolute coordinates*),Text[Style[tbarLabel,FontFamily->"Times",FontSlant->"Italic",FontSize->18],{1.6,3.6}]},AspectRatio->1]

Plot[{2 x,(x/2)},{x,-1,4},PlotRange->{{-1,4.5},{-1,4.5}},PlotStyle->Black,AxesLabel->{"x","t"},PlotLabels->{Placed[tbarLabel,Top],Placed[xbarLabel,Right]},AspectRatio->1]

(* Exer 1.3(d): Interval Delta s^2 = -1. These events are time-like separated from the origin. *)Plot[{Sqrt[x^2+1],-Sqrt[x^2+1]},{x,-3,3},PlotStyle->Black,AxesLabel->{"x","t"},AspectRatio->1]

(* Exer 1.3(e): Interval Delta s^2 = 1. These events are space-like separated from the origin *)Plot[{Sqrt[x^2-1],-Sqrt[x^2-1]},{x,-3,3},PlotStyle->Black,AxesLabel->{"x","t"},AspectRatio->1]

(* Exer.  1.3(f) Calibration ticks at 1m intervals along the x-bar and t-bar axes of (c). Below, The thin dotted lines are the invariant hyperbolae. The 6 hyperbolae crossing the t-axes have squared interval from the origin of {1,4,9} m^2,and thus mark of {-3,-2,-1,+1,+2,+3}m along both time axes.The 6 hyperbolae that cross the x-axes have negative squared interval {-1,-4,-9} m^2 and thus mark off {-3,-2,-1,1,2,3}m. *)

Plot[{2*x,(1/2)*x,Sqrt[x^2+1],-Sqrt[x^2+1],Sqrt[2^2+x^2],-Sqrt[2^2+x^2],Sqrt[3^2+x^2],-Sqrt[3^2+x^2],Sqrt[x^2-1],-Sqrt[x^2-1],Sqrt[x^2-2^2],-Sqrt[x^2-2^2],Sqrt[x^2-3^2],-Sqrt[x^2-3^2]},{x,-3.75,3.75},PlotRange->{{-3.75,3.75},{-3.75,3.75}},PlotStyle->{Directive[Black,AbsoluteThickness[2],Solid],(*2*x,using Dashed for visual distinctness as solid black is default for the first line*)Directive[Black,AbsoluteThickness[2],Solid],(*(1/2)*x*)Directive[Black,AbsoluteThickness[1],Dotted],(*Sqrt[x^2+1]*)Directive[Black,AbsoluteThickness[1],Dotted],(*-Sqrt[x^2+1]*)Directive[Black,AbsoluteThickness[1],Dotted],(*Sqrt[2^2+x^2]*)Directive[Black,AbsoluteThickness[1],Dotted],(*-Sqrt[2^2+x^2]*)Directive[Black,AbsoluteThickness[1],Dotted],(*Sqrt[3^2+x^2]*)Directive[Black,AbsoluteThickness[1],Dotted],(*-Sqrt[3^2+x^2]*)Directive[Black,AbsoluteThickness[1],Dotted],(*Sqrt[x^2-1]*)Directive[Black,AbsoluteThickness[1],Dotted],(*-Sqrt[x^2-1]*)Directive[Black,AbsoluteThickness[1],Dotted],(*Sqrt[x^2-2^2]*)Directive[Black,AbsoluteThickness[1],Dotted],(*-Sqrt[x^2-2^2]*)Directive[Black,AbsoluteThickness[1],Dotted],(*Sqrt[x^2-3^2]*)Directive[Black,AbsoluteThickness[1],Dotted] (*-Sqrt[x^2-3^2]*)},AxesLabel->{"x","t"},AspectRatio->1]

(* Exer. 1.3(g): Locus of points with interval from the origin, delta s^2 = 0. *)

AnswerEx1o3g=Plot[{x,-x},{x,-1,4},PlotStyle->{Directive[Black,AbsoluteThickness[0.5],Dashed] ,Directive[Black,AbsoluteThickness[0.5],Dashed]},AxesLabel->{"x","t"},BaseStyle->{FontFamily->"Helvetica",FontSize->14},(*For general plot text like numbers on axes*)LabelStyle->{FontFamily->"Times",FontSize->18,FontSlant->"Italic"},(*For AxesLabel text*)PlotRange->{{-1,4},{-2,4}},AspectRatio->Automatic (*To ensure correct scaling,as mentioned in previous answers*)]

AnswerEx1o3h=Plot[2,{x,-1,4},PlotStyle->Directive[Black,AbsoluteThickness[0.75],Solid] ,AxesLabel->{"x","t"},BaseStyle->{FontFamily->"Helvetica",FontSize->14},(*For general plot text like numbers on axes*)LabelStyle->{FontFamily->"Times",FontSize->18,FontSlant->"Italic"},(*For AxesLabel text*)PlotRange->{{-1,4},{-1,4}},AspectRatio->Automatic ]

AnswerEx1o3i=Plot[4/Sqrt[3]+(x-2/Sqrt[3])*(1/2),{x,-1,4},PlotStyle->Directive[Black,AbsoluteThickness[1.5],Solid] ,AxesLabel->{"x","t"},BaseStyle->{FontFamily->"Helvetica",FontSize->14},(*For general plot text like numbers on axes*)LabelStyle->{FontFamily->"Times",FontSize->18,FontSlant->"Italic"},(*For AxesLabel text*)PlotRange->{{-1,4},{-1,4}},AspectRatio->Automatic ]

AnswerEx1o3hyperbole=Plot[Sqrt[2^2+x^2],{x,-1,4},PlotStyle->Directive[Black,AbsoluteThickness[0.5],Dotted] ,AxesLabel->{"x","t"},BaseStyle->{FontFamily->"Helvetica",FontSize->14},(*For general plot text like numbers on axes*)LabelStyle->{FontFamily->"Times",FontSize->18,FontSlant->"Italic"},(*For AxesLabel text*)PlotRange->{{-1,4},{-1,4}},AspectRatio->Automatic ]

Show[AnswerEx1o3cbis,AnswerEx1o3g,AnswerEx1o3h,AnswerEx1o3i,AnswerEx1o3hyperbole]

Exercise 1.14

Suppose that the velocity v of referential Obar relative to O is small, |v| << 1. Show that the time dilation, Lorentz contraction and velocity-addition formulae can be approximated by, respectively: (see textbook)

(a) Time dilation (including relative error (d)):

(* Define the true function for the Lorentz factor. Note naming convention:
capital letters for the first letter is reserved for Mathematica build-in 
functions. *)
lfTrue[v_] := 1/Sqrt[1 - v^2]

(* Series expansion of Ftrue around v = 0 up to v^4 *)
Series[lfTrue[v], {v, 0, 5}]

(* We can understand the series expansion above by compting the 
first few terms of a Taylor series about v=0. E.g. the v^2, v^3, and v^4 
 *)
 (* The first term is simply the function evaluated at v=0*)
term1 = lfTrue[0]
(* The second (non-zero) term contains the second-order derivative of the function 
evaluated at v=0 *)
d2 = D[lfTrue[v], {v, 2}] // Simplify;
term2 = (d2 /. v -> 0)/Factorial[2]
(* The third order term contains the third-order derivative of the function 
. When evaluated at v=0 this term (and all odd-order terms) vanish !*)
d3 = D[lfTrue[v], {v, 3}] // Simplify;
d3 /. v -> v
term3 = (d3 /. v -> 0)/Factorial[3]
(* The fourth order term contains the fourth-order derivative of the function 
. When evaluated at v=0 this term (and all even-order terms) are non-zero *)
d4 = D[lfTrue[v], {v, 4}] // Simplify;
d4 /. v -> v
term4 = (d4 /. v -> 0)/Factorial[4]

(* Alternative: Taylor expansion of the Lorentz factor using x = v^2 . 
This leads to the same result but more directly, avoiding the odd power terms
that evaluated to zero above. *)
(* Here is the function definition for the Lorentz factor *)
lf[x_] := 1/Sqrt[1 - x]
(* Below we define the first derivative with respect to x:*) 
d1x = D[lf[x], {x}] // Simplify;
(* which we evaluate with x = v^2 to show its functional dependence *)
d1x /. x -> v^2
(* and then we evaluate it at x=0 to get the constant for that term in the Taylor series *) 
term1x = d1x /. x->0
(* Below we repeat the above steps for the second derivative, including the
2! in the denominator *)
d2x = D[lf[x], {x, 2}] // Simplify;
d2x /. x -> v^2
term2x = (d2x /. x -> 0)/Factorial[2]

(* Approximation to order v^2 *)
lfApprox[v_] := 1 + (1/2)*v^2

(* Finally, the result we are asked for: Relative error at v = 0.1 *)
RelativeError = N[(lfTrue[0.1] - lfApprox[0.1])/lfTrue[0.1]]

(b) Length contraction (including relative error (d))

(* The (exact) Length contraction factor, the factor to multiply the proper length to get
the contracted length observed when the rod is moving.  *)
lcfTrue[v_] := Sqrt[1 - v^2];
lcfTrue[x] /. x-> v

(* Find the Taylor series expansion for the above function to order 5 *)
Series[lcfTrue[v], {v, 0, 5}]

(* Length contraction factor: approximation to order v^2, by inspection of the
above result *)
lcfApprox[v_] := 1 - (1/2)*v^2;
lcfApprox[x] /. x-> v

(* The result we are asked for: Relative error at v = 0.1  *)
RelativeError = N[(lcfTrue[0.1] - lcfApprox[0.1])/lcfTrue[0.1]]

(c) Velocity addition (including relative error (d))

(* The (exact) velocity addition formula  *)
wprimeTrue[v_, w_] := (v + w)/(1 + v*w)
wprimeTrue[v,w]

(* Define the approximate expression *)
wprimeApprox[v_, w_] := (v + w)*(1 - v*w)
wprimeApprox[v, w]

(* Compute the relative error at v = 0.1, w = 0.1 *)
RelativeError = N[(wprimeTrue[0.1, 0.1] - wprimeApprox[0.1, 0.1])/wprimeTrue[0.1, 0.1]]

Chapter 12: Gravitational wave astronomy

Exercise 12.1

Use Eq. 12.2 to find the final gravitation wave frequencies at merger of pairs of equal-mass black holes of masses 10, million and billion solar masses, respectively.

(* The problem provides the key formula, which is also derived and discussed in
 Section 12.2 of the textbook.
The formula is given as: 
ffinal = 4.0 / Mchirp
where ffinal is the final frequency in kHz and Mchirp is the "chirp mass" of 
the system in units of solar masses. *)

(*First, define a Mathematica function based on Equation (12.2) from the 
textbook to easily calculate the final frequency for any given chirp mass.*)

ffinal[Mchirp_] := 4.0 / Mchirp

(* The formula for the chirp mass is given in Chapter 9, Equation (9.151):*)
(*Mchirp = (m1*m2)^(3/5) / (m1+m2)^(1/5)
where m1 and m2 are the component masses of the binary system, m1 >= m2 by 
convention. *)

(* Second, define a Mathematica function to calculate the chirp mass *)
chirpMass[m1_, m2_] := (m1 * m2)^(3/5) / (m1 + m2)^(1/5)

(* Calculate the chirp mass for the 3 binary systems: *)  

(* S1: m1 = m2 = 10 solar masses *)

mChirpS1 = chirpMass[10, 10]
N[mChirpS1]

(* Finally frequency in kHz *)
freqS1 = ffinal[mChirpS1]

(* S2: m1 = m2 = 10^6 solar masses *)
mChirpS2 = chirpMass[10^6, 10^6]
(* Finally frequency in kHz *)
freqS2 = ffinal[mChirpS2]

(* S3: m1 = m2 = 10^9 solar masses *)
mChirpS3 = chirpMass[10^9, 10^9]
(* Finally frequency in kHz *)
freqS3 = ffinal[mChirpS3]

Exercise 12.7 (f)

Here we explore solving Eq.(12.26) for the inclination angle iota.

(f) This result shows that iota needs to get fairly large before epsilon becomes measurable, or equivalently the ratio r becomes measurably different from 1. Suppose that r is measured to be 1, with an error of 1%. Show that iota might be as large as 0.53 radians or 30°.

(*Exercise 12.7(f) asks us to explore the uncertainty in measuring the 
inclination angle, ι, of a nearly face-on binary system. 
Specifically, it states:
"Suppose that r is measured to be 1, with an error of 1%. Show that ι might 
be as large as 0.53 radians, or 30°."*)
(*Here, `r` is the ratio of the measured polarization amplitudes, |h+|/|h×|, 
and ι is the inclination angle of the binary's orbit (where ι = 0 is a 
face-on orbit).*)
(*The problem states that `r` is measured to be 1 with an error of 1%. We have constructed
r >= 1 so this means the likely measured value could be up to 1.01 *)
(* We could simply apply the arccos function to Scott eqn(12.37):
z = r - sqrt(r^2 - 1) = cos(ι)
ι = arccos(z)
*)
inclinationAngle[r_] :=ArcCos[r-Sqrt[r^2-1]]
iotaDirect=inclinationAngle[1.01]

The above is the inclination in radians. Below we convert to degrees.

(* The above result is the maximum likely inclination angle in radians. In degrees: *)
iotaDirect*180/Pi

The above result is approximately 30°, confirming the problem statement.

The above is a legitimate solution, but does not follow the spirit of the exercise. The Exercise 12.7(f) builds directly on the results from the previous parts of Exercise 12.7, in which explore the relationship between the measured amplitude ratio `r` and the inclination `ι` (Greek letter iota). A key relationship for this part comes from the approximation derived in Exercise 12.7(e) for nearly face-on systems (where `r` is close to 1 and `ι` is small). In Exercise 12.7(d) we defined variable, ϵ, as the deviation of the measured 
amplitude ratio from 1:
ϵ = r - 1
Then in Exercise 12.7(e) we showed that for small ι, this is related to the inclination by:
ϵ ≈ ι^4 / 8
(See tip in Exercise 12.11(d) for entering Greek letters.)
As stated above, the problem gives that `r` is measured to be 1 with an error of 1%.
 We have constructed r >= 1 so this means the likely measured value could be up to 
 1.01 and the maximum  likely value for ϵ is 0.01.

Clear[ι]   (* This is how we clear a previously defined variable *)

(*
We can now set up an equation in Mathematica to solve for the inclination 
angle ι that corresponds to this maximum error using the result of Exer 12.7(e).*)

(* The maximum deviation from r=1 is 1%, or 0.01 *)
epsilon = 0.01;

(* The relationship from Exercise 12.7(e) *)
equation = epsilon == ι^4 / 8;

(* Of course we could invert this relation to find ι , but as a demonstration
of Mathematica, let's use its Solve function. This remindes us that we must be
explicit about wanting a real (and not complex) answer. Since ι must be a real, 
positive angle, we will select the appropriate solution as indicated below *)

(* Solve the equation for ι over the real numbers *)
solutions = Solve[equation, ι, Reals]

(*The positive solution is the physically meaningful one for the inclination angle.*)

(* Here is how we could extract the positive solution for ι in radians *)
iotaInRadians = ι /. solutions[[2]]

(*This matches the value of 0.53 radians given in the problem statement. Now, let's convert this to degrees to confirm the second part of the question.*)

The above is the inclination in radians. Below we convert to degrees.

(* Convert the result from radians to degrees *)
iotaInDegrees = iotaInRadians * (180 / Pi)

The above result is approximately 30°, again confirming the problem statement.

Exercise 12.7(g)

Show that the percentage change of 1 + cos^2(ι) in this range is about 7%. This will be the uncertainty in the inferred distance for a face-on signal, even though the amplitudes themselves are measured to 1% accuracy.

This part investigates how the uncertainty in the inclination angle (found in part f) affects the inferred amplitude of the gravitational wave, and therefore the inferred distance.

From Chapter 9, Schutz Eq. (9.100), the amplitude of the h+ polarization is proportional to the factor (1 + Cos[ι]^2). We need to calculate how much this factor changes as inclination varies from iota = 0 (perfectly face-on) to the maximum value of iota = 0.531829 radians that we found in Exer. 12.7(f).

(* Define a function for the amplitude factor *)
amplitudeFactor[ι_] := 1 + Cos[ι]^2

Now we apply this function at the two extreme values of inclination:

(* Calculate the factor for a perfectly face-on system (iota = 0) *)
factorFaceOn = amplitudeFactor[0]

(* Calculate the factor at the maximum possible inclination from part (f) *)
factorMaxIota = amplitudeFactor[0.5318295896944989]
(* Calculate the factor at the midway inclination from part (f) *)
factorMidwayIota = amplitudeFactor[0.5318295896944989/2]

Now we calculate the percentage change between these two values. The percentage change is the absolute difference divided by the mid-point value (midway between face-on and maximum likely inclination), multiplied by 100 to convert to percent:

percentageChange = Abs[factorMaxIota - factorFaceOn] / factorMidwayIota * 100

The result is approximately 13.3%, which is almost twice the value of “about 7%” stated in the problem. A likely explanation is that the 7% uncertainty is about the midway value, so we can justify dividing by two:

percentageChange/2

Rounding off, we confirm an uncertainty of about 7% about the midway inclination.

Discussion 

This pair of exercises provides a crucial insight into the challenges of gravitational wave astronomy. As discussed in Section 12.3 of the textbook, the inclination angle is a key parameter needed to accurately determine the distance to a standard siren.

Our calculations show that for nearly face-on systems (which produce the strongest and most easily detectable signals), even a small 1% uncertainty in the measured polarization amplitudes leads to a very large uncertainty in the inclination angle (up to 30°).

Furthermore, this large uncertainty in iota propagates into a significant uncertainty (about 7%) in the inferred amplitude of the wave. Since the inferred distance is inversely proportional to the measured amplitude, this translates directly into a ~7% uncertainty in the distance measurement. This makes it difficult to get a precise measurement of cosmological parameters like the Hubble-LemaΙtre constant from a single, nearly face-on event.

Exercise 12.9

Here we examine how an optimum all-sky search might correct for signal modulation and pulsar spin-down, and we use that to do a rough estimation of the number of arithmetic operations that would be required to search for a gravitational wave pulsar in six months of data. Each arithmetic operation is called, in the language of computing, a floating-point operation. For the purposes of this calculation we assume that the search goes up to 2 kHz in gravitational wave frequency, because theoretically it would be possible for a neutron star to spin at 1 kHz. We assume that the data are sampled at 4 kHz, which is the minimum sampling rate to allow a 2 kHz bandwidth search. The method will be to demodulate the recorded data for a particular patch on the sky, as described in the text, then to compensate for a particular spin-down value by the method described below, and then to do a Fourier transform for each combination of sky patch and spin-down value, in order to find a signal. Since this is a rough estimate, we regularly round our results to one significant figure. It should be emphasized that there can be other approaches to this problem that are more efficient, but not by huge factors. The approach here is straightforward to describe and shows rather readily that an optimum search is far out of reach.

(a) Show that if the duration of the observation is s, the number of data samples is , and the frequency resolution of the observation is

(*Part (a): Number of Data Samples*)

(*The problem asks for Ns, the number of samples in a data stream of 
duration T = 0.15e8 seconds (about one half year) with a sampling rate of 
4000 Hz.*)

(* Duration of the observation in seconds, Tobs *)
Tobs = 0.15*10^8;
(* Sampling rate in Hz *)
samplingRate = 4000;
(* Number of samples Ns *)
Ns = Tobs * samplingRate

(* For the sake of curiousity, here is a much more precise calculation *)
secperyear = 365.2425 * 24 * 3600;
Nsprecise = 0.4*^4 * secperyear * (1/2)

The frequency resolution is calculated from the reciprocal of the longest period that is resolved, which here is 6 months. Schutz is more conservative and  assumes half this resolution.

(* Frequency resolution, Delta f in Hz*)
Deltaf = 1/Tobs

(b) ... Assuming that the interpolation requires of order ten operations per data point, show that each demodulation requires  operations.There is, of course, one demodulation for each resolvable patch on the sky.

(* Number of operations to demodulate the Doppler shift for data from a given patch
of the sky, Ndemod *)
Ndemod = 10*Ns

(c) To correct for the pulsar’s changing signal frequency, we will assume that the only parameter we need is the first time-derivative of the frequency, i.e. that the frequency changes linearly in time over the observation period. We shall handle this in the same way as the demodulation, by shifting to a different time measure. We assume that the source’s wavefronts are emitted at a uniform rate with respect to a clock that is slowing down. Once we map that clock’s time onto barycenter time, we restore the pulsar to a single-frequency signal and we can use a Fourier transform to search for it. Letting  represent the start of the observation time in the barycenter, we can give concrete expression to this by relating the pulsar’s `clock’  to barycentric time t as follows:

where  is a parameter known as the slowdown time, a measure of how long it takes the pulsar to change its frequency by a significant amount. Since realistic slowdown times are larger than  years, we always have , so that the fractional changes in clock rates are small.

(*Part (c): Computational Cost Estimation*)

(*This part is broken into several steps to build up the final estimate of 
the total computational cost.*)

(* Define constants from the problem *)
Tobs = 0.15*^8; (* Observation time for this part (6 months) *)
f0 = 2000; (* Maximum frequency to search in Hz *)
Deltaf = 1.3*^-7; (* Frequency resolution from a 6-month observation *)
numSkyPatches = 0.3*^13; (* Number of independent sky patches to search *)

(*Sub-part (iii): Maximum Spin-Down*)
(*We need to find the maximum inverse spin-down time, (1/τ)_max, for a 
pulsar that slows from 2000 Hz to 1000 Hz in 1000 years.*)

fInitial = 2000; (* Hz *)
fFinal = 1000; (* Hz *)
timeToSlow = 1000 * secperyear; (* seconds *)

(* The rate of change of frequency is df/dt = (fFinal - fInitial)/timeToSlow *)
dfdt = (fFinal - fInitial) / timeToSlow;
(* From the text, df/dt = -2*f0/τ. We solve for 1/τ. *)
invTauMax = Abs[dfdt / (2*f0)]
(* Compare with Scott eqn.(12.70)*)

(*Sub-part (iv): Spin-Down Resolution*)
(*Next, we find the resolution with which we can measure the spin-down 
parameter, Δ(1/τ).*)

deltaInvTau = Deltaf / (2 * Tobs * f0)

(* Compare with Scott eqn.(12.72) for the conservative resolution of Schutz. 
Compare with Scott eqn.(12.73) for the standard resolution.*)

(*Sub-part (vi): Number of Independent Searches*)
(*The total number of searches is the number of sky patches multiplied by the 
number of spin-down templates  we need to check. 
The number of spin-down templates, n in Scott eqn(12.75), is the total range 
of spin-down values divided by the resolution.*)

numSpinDownTemplates = invTauMax / deltaInvTau  (* compare with Scott eqn(12.75)*)
numberSearches = numSkyPatches * numSpinDownTemplates (* compare with Scott eqn(12.76)*)

(*Sub-part (vii): Flops per Search*)
(*The exercise provides a formula for the number of floating-point operations 
per search.*)
(* We need the natural lorgarithm of the number of the number of samples Ns *)
N[Log[Nsprecise],6]
(* We use the Nsprecise value from part (a) *)
flopsPerSearch = (Log[Nsprecise] + 10 + 10) * Nsprecise (* Compare with Scott eqn.(12.77) *)

(*Sub-part (viii): Total Computing Resources*)
(*Finally, we calculate the total flops required for the entire search over six
 months and determine how many 400 petaflop supercomputers would be needed to 
 complete it in time.*)

(* Total flops needed for the search *)
flopsEverySixMonths = flopsPerSearch * numberSearches (* Scott eqn(12.78)*)

(* Convert this to a rate in flops per second *)
secondsInSixMonths = 0.5 * secperyear;
flopsPreciseRate = flopsEverySixMonths / secondsInSixMonths
(* Scott eqn.(12.79)*)

(* Define the speed of one supercomputer in flops/sec *)
petaflop = 10^15;
supercomputerSpeed = 400 * petaflop;

(* Calculate the number of machines required *)
numMachines = flopsPreciseRate / supercomputerSpeed

Conclusion for Exercise 12.9
Even correcting the factor of ten thousand exaggeration of the question, the final result is still staggering. To perform an optimal, all-sky search for a single unknown pulsar in six months of data would require nearly 5 million (and not 50 billion) of the world’s fastest supercomputers.
This calculation powerfully illustrates why such “blind” searches are computationally impossible to perform optimally. Instead, researchers must use hierarchical methods and other clever algorithms to conduct sub-optimal searches, which trade some sensitivity for a computationally feasible problem. It also highlights the importance of targeted searches, where the sky position and frequency are already known from radio pulsar observations, dramatically reducing the number of required templates.

Exercise 12.11(d)

(d) In Schutz’s Table 12.1, look at the masses given for GW170817. The collaboration published uncertainty ranges for these values as well (abbott2019properties). The uncertainty in the chirp mass was one part in , so for the purposes of this exercise, we take it as exactly known. The 90% uncertainty range of the total mass  was asymmetrical, given as . Calculate the corresponding 90% bounds on the component masses  and . (These are only approximate bounds; a full analysis would take into account the uncertainty in the chirp mass and its correlation with the uncertainty in the total mass.).

(*From Table 12.1 in the textbook, we have the best-fit masses for the two 
neutron stars in GW170817.*)

m1 = 1.44; (* Mass of the primary component in solar masses *)
m2 = 1.28; (* Mass of the secondary component in solar masses *)

(*From these, we can calculate the total mass (M) and the reduced mass (μ).*)

M = m1 + m2            (* Scott eqn.(12.100) *)
μ = (m1 * m2) / M      (* Scott eqn.(12.101) *)

Tip: 
To get the Greek letter mu above we have several options
1) use the syntax :
	\[x] 
and replace “x” with the name of the Greek letter in capital letters, e.g. replace “x” with “Mu” gives a lowercase mu, μ.
2) use keyboard shortcuts:
	esc m esc
gives μ.
For ψ type:
	esc psi esc
3) Using palettes. Select Insert from the Default toolbar at the top of the page. (Make sure the Default toolbar appears at the top with “Evaluation” at the top left, and Insert at the top right. If not, from the main menu, Window, Toolbar, select Default.)

(*The formulas from part (b) of the exercise give the rate of change of the 
component masses with respect to a change in the total mass.*)

dm1dM = 1/2 + (1/2 - μ/(3*M)) / Sqrt[1 - 4*μ/M];
dm2dM = 1/2 - (1/2 - μ/(3*M)) / Sqrt[1 - 4*μ/M];


(*Calculating the Uncertainty Bounds*)

(*The problem gives the asymmetric 90% confidence interval for the total mass.*)

deltaMUpper = 0.04; (* Upper uncertainty bound for M *)
deltaMLower = -0.01; (* Lower uncertainty bound for M *)

(*We can now calculate the corresponding uncertainties for m1 and m2 by 
multiplying the derivatives by these δM values.*)

(* Uncertainty for m1 *)
deltaM1Upper = dm1dM * deltaMUpper
deltaM1Lower = dm1dM * deltaMLower
(* Compare with result just after Scott eqn.(12.103)*)

(* Uncertainty for m2 *)
deltaM2Upper = dm2dM * deltaMUpper
deltaM2Lower = dm2dM * deltaMLower
(* Compare with result just after Scott eqn.(12.104)*)

(*Finally, we apply these uncertainties to the original mass values to find the 
90% confidence bounds.*)
(* For the primary mass *)
m1LowerBound = m1 + deltaM1Lower
m1UpperBound = m1 + deltaM1Upper
(* Compare with result just after Scott eqn.(12.103)*)

(* For the secondary mass *)
m2LowerBound = m2 + deltaM2Lower
m2UpperBound = m2 + deltaM2Upper
(* Compare with result just after Scott eqn.(12.104)*)

(*7. Conclusion for Exercise 12.11*)
π
ψ

(* ::Text:: *)
(*The final 90% confidence intervals are approximately:*)
(*m1: [1.34, 1.84] Mö*)
(*m2: [0.92, 1.37] Mö*)
(*This calculation demonstrates an important concept in parameter estimation 
for binary systems. Even though the total mass M is determined with a relatively
 small uncertainty (in this case, +0.04/-0.01 M_O), the uncertainties in the 
 individual component masses can be much larger.*)
(*The reason, as shown by the formulas for dm1/dM and dm2/dM, is that for nearly equal-mass systems (where μ/M is close to its maximum value of 1/4), the denominator Sqrt[1 - 4μ/M] becomes very small. This amplifies the effect of any uncertainty in M, leading to larger error bars for m1 and m2. This is a common feature in gravitational wave data analysis: some parameters (like the chirp mass) are measured very precisely, while others that depend on them can have significantly larger uncertainties.*)

Supplementary Problem 12.3

Here we consider the effect of the movement of a gravitational wave detector during a detection event and compare this to the wavelength of the gravitational waves. The duration of the detection event and the wavelength of the gravitational waves depends upon the source. To be concrete, consider the merger of neutron stars with masses . We saw in Exercise 12.2 that the detection event duration in this case is about 150 s.

(a) Calculate the minimum wavelength of the gravitational waves.

(* Define Necessary Constants*)

(*First, we define the physical and astronomical constants used in the 
calculations. These values are taken from earlier sections of the textbook for 
consistency. The `N[]` function is used to get numerical values where 
appropriate.*)

cSI = 299792458; (* Speed of light in m/s *)
secperday = 3600*24;
secperyear = 365.2425*secperday;
auSI = 149597870700; (* Astronomical Unit in meters *)
lightyearSI = 365.2425*secperday*cSI;
earthOmegaSI = 2*Pi/(23*3600 + 56*60 + 4.091); (* Earth's rotational angular velocity in rad/s *)
earthEquatorialRadiusSI = 6378137; (* Earth's equatorial radius in meters *)

(* Chirp Mass of a Neutron Star Binary*)

(*Maple uses `:=` for assignment and `subs` for substitution. 
In Mathematica, we define a function for the chirp mass using `:=` and then 
call that function with the specific mass values.*)

(* Chirp Mass Function (from Schutz Eq. 9.151) *)
mChirp[m1_, m2_] := (m1*m2)^(3/5) / (m1 + m2)^(1/5)

(* Calculate the chirp mass for two 1.4 solar mass neutron stars *)
mChirpValue = mChirp[1.4, 1.4]

(* Final Frequency and Wavelength*)

(*We now calculate the final frequency (`ffinal`) in kHz and the corresponding 
wavelength (`lambdafinal`) in meters.*)

(* The final frequency in kHz (from Schutz Eq. 12.2) *)
ffinal = 4.0 / mChirpValue;  (* suppress output with ";" *)

(* The final frequency in Hz *)
ffinalHz = ffinal * 1000   (* compare to Scott eqn(12.107) *)

(* The corresponding wavelength in meters *)
lambdafinalSI = cSI / ffinalHz (* compare to Scott eqn(12.108) *)

(b) Calculate how much the LIGO Livingston detector has moved relative to the Sun during this detection event; express the result as a fraction of the minimum wavelength. Is the rotation of the Earth about its axis important?

(* Velocity and Displacement Calculations*)

(*This section calculates various velocities and compares the resulting 
displacements over a 150-second interval to the final gravitational wavelength.*)

(*Earth's Velocity in Orbit*)

(*Maple's `evalf` is equivalent to Mathematica's `N[]` for numerical evaluation.*)

earthOrbitVelocitySI = N[2*Pi*auSI] / secperyear (* Scott eqn.(12.109) *)

earthDisplacementSI = 150 * earthOrbitVelocitySI
earthDisplacementSI / lambdafinalSI

(* Earth's equatorial speed from rotation *)
earthEquatorialVelocitySI=N[earthOmegaSI*earthEquatorialRadiusSI]
(* compare with Scott eqn.(12.110) *)

equatorialDisplacementSI = 150 * earthEquatorialVelocitySI
equatorialDisplacementSI / lambdafinalSI

(* Livingston detector at 30 deg North*)
livingstonDisplacementSI = Cos[30*Pi/180]equatorialDisplacementSI 

livingstonDisplacementSI/ lambdafinalSI

(c) How much has the Sun moved relative to the Galactic centre during this neutron star merger detection event?

(*Sun's Orbital Velocity*)
sunOrbitVelocitySI = N[(2*Pi*28000)*lightyearSI] / (0.240*^9*secperyear)
sunDisplacementSI = 150*sunOrbitVelocitySI
sunDisplacementSI/lambdafinalSI

(d) Finally, what about the motion of the Milky Way galaxy? How much has the Milky Way moved during this neutron star merger detection event, and with respect to what reference frame?

(*Milky Way's Velocity*)

milkyWayVelocitySI = 600000;
milkyWayDisplacementSI = 150*milkyWayVelocitySI
milkyWayDisplacementSI/lambdafinalSI

Chapter 13: Cosmology

Exercise 13.6

Prove the statement leading to 
   , 
 			Schutz Eq. (13.9)

that we can deduce  of our three-spaces by setting  in the expression for the Einstein tensor (Schutz Eqs 10.15-17) for the spacetime with line element for a spherically symmetric star (Schutz Eq. (10.7)).

As an exercise in digging into the basics with Mathematica, below we use Mathematica’s core tools to calculate the Christoffel symbols important tensors of GR from the metric tensor. (Of course, there are more streamlined and powerful tools in Mathematica to compute these and much more. But below, we’ll see some great demonstrations of the power and simplicity of Mathematica.)

(*1. Setting up the Metric in Mathematica*)

(*First, we clear all previous definitions to ensure a clean start. Then, we define 
our coordinates and the 3D metric tensor, g. We will use a general function, λ[r], 
for the g_rr component. This corresponds to the spatial 
line element in Schutz's Eq.(13.8) above:*)
(*dl^2 = Exp[2λ[r]] dr^2 + r^2 (dθ^2 + Sin[θ]^2 dϕ^2)*)

ClearAll["Global`*"];

(* Define the coordinates for our 3D space *)
coords = {r, θ, ϕ};

(* Define the metric tensor g_ij as a 3x3 matrix *)
g = {{Exp[2*λ[r]], 0, 0}, 
     {0, r^2, 0}, 
     {0, 0, r^2 * Sin[θ]^2}};

MatrixForm[g]

(*2. Calculating the Riemann Tensor from First Principles*)

(*To calculate the Riemann tensor, R^ρσμν, we need to follow a series of steps:*)
(*1. Calculate the inverse metric, g^ij.*)
(*2. Calculate the Christoffel symbols of the second kind, Γ^ρμν.*)
(*3. Calculate the Riemann tensor using the Christoffel symbols.*)
(*This process can be tedious and error-prone by hand but is straightforward for Mathematica.*)

(*Step 2.1: Inverse Metric*)

(* Calculate the inverse of the metric tensor, g^ij *)
ginv = Simplify[Inverse[g]];
MatrixForm[ginv]

(*Step 2.2: Christoffel Symbols*)

(*The Christoffel symbols are calculated from the derivatives of the metric tensor
 using the formula:*)
(*Γ^k_ij = (1/2) * g^kl * (g_jl,i + g_il,j - g_ij,l)*)

(* First, calculate all the first partial derivatives of the metric *)
dg = D[g, {coords}]
(* This creates a nested list or table with three indices or iterators.
The convention for the nesting order is (like a Python Numpy array) that the 
first index is the outermost loop (slowest varying) and the last index
is the innermost loop (fastest varying).
(* Mathematica's D function adds the new index from the differentiation as the
 last or innermost index. *)

(* Now, calculate the Christoffel symbols Γ^k_ij *)
(* Table[expression, {k, 1, 3}, {i, 1, 3}, {j, 1, 3}]
creates a nested list (an array or tensor) by iterating through the indices
 k, i, and j. 
Because k appears first, it is the outer-most (slowest) index as we explained above. *)
Sum:  the summation must be stated explicitly, and notice we sum
explicitly over l *)
(* The syntax for Sum is:
Sum[expression, {iterator, min, max}] *)
christoffel2 = Table[
   (1/2) * Sum[
     ginv[[k, l]] * (dg[[j, l, i]] + dg[[i, l, j]] - dg[[i, j, l]]),
     {l, 1, 3}
   ],
   {k, 1, 3}, {i, 1, 3}, {j, 1, 3}
   ];

(* Display all 7 non-zero Christoffel symbols to check our work *)
(* Christoffel symbol Γ^r_rr *)
Simplify[christoffel2[[1, 1, 1]]]
(* Christoffel symbol Γ^r_θθ *)
Simplify[christoffel2[[1, 2, 2]]]
(* Christoffel symbol Γ^r_ϕϕ *)
Simplify[christoffel2[[1, 3, 3]]]
(* Christoffel symbol Γ^θ_rθ *)
Simplify[christoffel2[[2, 1, 2]]]
(* Christoffel symbol Γ^θ_ϕϕ *)
Simplify[christoffel2[[2, 3, 3]]]
(* Christoffel symbol Γ^ϕ_rϕ *)
Simplify[christoffel2[[3, 1, 3]]]
(* Christoffel symbol Γ^ϕ_θϕ *)
Simplify[christoffel2[[3, 2, 3]]]

(*Step 2.3: The Riemann Tensor*)

(*Finally, we calculate the Riemann tensor, R^ρσμν, using its definition in 
terms of the Christoffel symbols and their derivatives:*)
(*R^k_lmn = Γ^k_ln,m - Γ^k_lm,n + Γ^k_sm * Γ^s_ln - Γ^k_sn * Γ^s_lm*)

(* First, we need the derivatives of the Christoffel symbols *)
dchristoffel2 = D[christoffel2, {coords}];

(* Now we can assemble the Riemann tensor, R^k_lmn *)
riemann1up3down = Table[
   dchristoffel2[[k, l, n, m]] - dchristoffel2[[k, l, m, n]] +
    Sum[
     christoffel2[[k, s, m]]*christoffel2[[s, l, n]] - 
      christoffel2[[k, s, n]]*christoffel2[[s, l, m]],
     {s, 1, 3}
     ],
   {k, 1, 3}, {l, 1, 3}, {m, 1, 3}, {n, 1, 3}
   ];
(* Find Ricci tensor: Contract the first and third indices of the riemman tensor *)
ricci2down = Table[Sum[riemann1up3down[[s,i,s,j]],
     {s, 1, 3}
     ],
   {i, 1, 3}, {j, 1, 3}
   ];
(*3. Displaying the Result*)

(*The Riemann tensor has 3^4 = 81 components. Most of them are zero. *)

simplifiedAllComponents = Simplify[ricci2down]

(* Find Ricci scalar: Raise the first index, then contract the first and second indices of the ricci tensor *)
ricci1up1down = Simplify[Table[Sum[ginv[[i,s]]*ricci2down[[s,j]],
     {s, 1, 3}
     ],
   {i, 1, 3}, {j, 1, 3}
   ]];
ricciscalar = Sum[ricci1up1down[[s,s]],
     {s, 1, 3}];
(*3. Displaying the Result*)

(*The Ricci scalar. *)

ricciscalar = Simplify[ricciscalar]

(* Einstein tensor: the trace-reverse of the Ricci scalar. *)
einstein2down = Simplify[Table[ricci2down[[i,j]]-(1/2)*g[[i,j]]*ricciscalar,
{i,1,3},{j,1,3}]]
(* *)

The result above corresponds to the Einstein tensor stated in the problem.

Exercise 13.7

Show that the metric with line element 
({{ ((dl)^2) ,  ==(e^(2Λ(r))) ((dr)^2)+(r^2)((dΩ)^2) ,  Schutz Eq. (13.8) }})
is only flat at  if  in

Recall this was used in the derivation of the FRW metric in Schutz . 
For the geometry to be flat we require the Riemann tensor to be zero, see ,

Lowering and then raising the first index in Schutz Eq.(6.71) we easily confirm that

which will simplify our calculations. 
We calculated the Riemann tensor for the line element Schutz Eq. (13.8) in Exercise 13.6 above.

(*1. execute the cell of the previous Exercise 13.6 *)
(* which will calculate the Riemann Tensor from First Principles *)

(*Step 2.3: The Riemann Tensor*)

(*Finally, we calculate the Riemann tensor, R^ρσμν, using its definition in 
terms of the Christoffel symbols and their derivatives:*)
(*R^k_lmn = Γ^k_ln,m - Γ^k_lm,n + Γ^k_sm * Γ^s_ln - Γ^k_sn * Γ^s_lm*)

(* First, we need the derivatives of the Christoffel symbols *)
dchristoffel2 = D[christoffel2, {coords}];

(* Now we can assemble the Riemann tensor, R^k_lmn *)
riemann1up3down = Table[
   dchristoffel2[[k, l, n, m]] - dchristoffel2[[k, l, m, n]] +
    Sum[
     christoffel2[[k, s, m]]*christoffel2[[s, l, n]] - 
      christoffel2[[k, s, n]]*christoffel2[[s, l, m]],
     {s, 1, 3}
     ],
   {k, 1, 3}, {l, 1, 3}, {m, 1, 3}, {n, 1, 3}
   ];
(* Find Ricci tensor: Contract the first and third indices of the riemman tensor *)
riemann4down = Table[Sum[g[[s,k]]*riemann1up3down[[s,l,m,n]],
     {s, 1, 3}
     ],
   {k, 1, 3}, {l, 1, 3},{m, 1, 3},{n, 1, 3}
   ];
(*3. Displaying the Result*)

(*The Riemann tensor has 3^4 = 81 components. Most of them are zero. *)

simplifiedAllComponents = Simplify[riemann4down]

To compare the Riemann tensor above with that found in Scott eqn(13.17) we could type pull out the desired components various ways:

simplifiedAllComponents [[3,2,3,All]]

simplifiedAllComponents [[1,3,1,3]]

We can use Mathematica to find and display only the unique, non-zero components. But first we must convert the nested list of lists into one big list of elements:

simplifiedUniqueComponents = DeleteDuplicates[Flatten[simplifiedAllComponents]]

nonZeroComponents = Select[simplifiedUniqueComponents,UnsameQ[#,0]&]

Mathematica Tip: We cannot simply compare the elements of the list directly to zero, as in the following line of code. This fails because the elements are symbolic expressions and Mathematica doesn’t know if they are 0 or not (we haven’t told it what lambda(r) is yet). The result of such an uncomparable expression is nothing.

nonZeroComponents = Select[simplifiedAllComponents, # != 0 &]

Mathematica apparently knows how to simplify the following:

Simplify[Exp[-F[r]]*Exp[F[r]]]

Mathematica Limit: And yet it didn’t know that the last pair of elements were equivalent to the previous pair in “nonZeroComponents” above.

(*As the output shows, the components of the Riemann tensor are generally **not 
zero**. They depend on the function λ[r] and its derivative.*)
(*To complete the exercise and check for flatness, one must substitute the 
specific form for λ[r] that corresponds to a homogeneous and isotropic space, 
which is given in Schutz's Equation (13.12):*)
(*Exp[2λ[r]] = 1 / (1 + (1/3) k r^2  - A/r )*)

Exercise 13.15

The integral over the black-body spectrum times the energy per photon , where  is Planck’s constant and  is the photon frequency) gives the total radiation energy density

where  is the radiation constant. Here is the integral over all frequencies:

(* Assumptions *)
$Assumptions = T > 0 && kB > 0 && h > 0;

(* Planck spectrum integral over frequency *)
Integrate[
  (8 Pi h ν^3)/(Exp[h ν/(kB T)] - 1),
  {ν, 0, Infinity},
  Assumptions -> $Assumptions
]

Exercise 13.19

From our earlier result in Scott eqn.(13.56), we have

Finding a numerical value for the age of the Universe, we must be careful of the units. An approximate value of the Hubble constant is . Recall, 1 parsec = 30 856 775 814 671.9 km.
This gives :

(* First execute the cell in Chapter "Fundamental constants" *)
t0 = 2/(3*hubbleSI)

UnitConvert[t0,"Years"]

Exercise 13.23

Calculate the redshift of decoupling by assuming that the cosmic microwave radiation has temperature  K today and had the temperature  at decoupling, where eV is the energy needed to ionize hydrogen, see Exer. 13.22(c), [and eV/K is the Boltzmann constant.]

The cosmic microwave radiation has been out of equilibrium with matter since the time of decoupling, and has cooled in proportion to the expansion of the universe,

T = ((T0 R0)/R)

where  and  are the radiation temperature and expansion factor today, see Exer.~13.15.

(* by entering the 2.7 as a Quantity with units "Kelvins", Mathematica can keep 
track of the units *)
T0 = Quantity[2.7,"Kelvins"]

Using the suggestion  for the temperature at decoupling, we find a temperature at the time of decoupling

(* by entering the 13.6 as a Quantity with units "eV", Mathematica telling us 
below that the temperature is in Kelvin *)
T = Quantity[13.6,"eV"]/(20*kBeVperK)

Why the factor of 20? Schutz discusses this in his Exer.~13.22(c), where he gives a more precise value of , which implies a decoupling temperature of .
Using the redshift formula, (Scott eqn. (13.26), Schutz 13.67))

Substituting the temperatures for the scale factors using the first equation above we find:

(* Redshift z at the time of decoupling: *) 
z = T/T0 - 1

Exercise 13.25

Estimate the times earlier than which our uncertainty about the laws of physics prevents us drawing firm conclusions about cosmology as follows.

(a) Deduce that, in the radiation-dominated early universe, where the curvature term depending on the curvature constant  is negligible, the temperature  behaves as

where  is Boltzmann’s constant.

See Scott eqn. (13.79)

(* We have defined precise values of the standard physical constants in Chapter
Fundamentals above.  *)
beta =   (45 * hbarSI^3 * cSI^5 / (32 * Pi^3 * GSI * kBSI^4))^(1/4);
(* Note that when we evaluate beta we get both a numerical value and its units *)
N[beta,4]
(* However, the units seems complicated. We anticipated K/sqrt(s), let's check! *)
UnitConvert[beta,"Kelvins"*"Seconds"^(1/2)]
QuantityMagnitude[beta]

(* Liddle includes the effects of neutrinos with a numerical factor *)
betaLiddle = 
  (45 * hbarSI^3 * cSI^5 / (32 * Pi^3 * GSI * 1.68 * kBSI^4))^(1/4)
UnitConvert[betaLiddle,"Kelvins"*"Seconds"^(1/2)]

(b) Assuming that our knowledge of particle physics is uncertain for GeV, find the earliest time  at which we can have confidence in the physics.

(* We set kB T = 1000 GeV in the relation above for T as a function of t, 
solve for t: *)
t = (beta * kBeVperK/Quantity[10^12,"Electronvolts"])^2
(* Confirm the units are seconds *) 
UnitConvert[t,"Seconds"]

(c) Quantum gravity is probably important when a photon has enough energy  to form a black hole within one wavelength . Show that this gives . This is the Planck temperature. At what time  is this an important worry?

(* We set kB T = sqrt(h c^5/G) see Scott eqn.(13.83)
*)
energyBlackHolePhoton = Sqrt[hSI*cSI^5/GSI]
UnitConvert[energyBlackHolePhoton,"Joules"]
UnitConvert[energyBlackHolePhoton,"eV"]
(*in the relation above for T as a function of t, 
solve for t: *)
t = (beta * kBSI/Sqrt[hSI*cSI^5/GSI])^2
(* Confirm the units are seconds *) 
UnitConvert[t,"Seconds"]

Supplementary Problem 13.7(d)

We need to perform the following integral

Integrate[Sin[x/R]^2, {x,0,r}]