🗎 Code for the
static Big Bang spacetime diagrams
(* | Solver for the critical value of ΩΛ in a closed universe || yukterez.net | *) Ωk = 1-Ωm-Ωr-ΩΛ; pr = Ωr/3; pm = 0; H² = Ωr/a^4+Ωm/a^3+Ωk/a^2+ΩΛ; ȧ² = H² a^2; ä = a (ΩΛ-(Ωm+pm)/a^3/2-(Ωr+pr)/a^4/2); Ωr = 3/10; Ωm = 11/10; f = N[Reduce[ȧ²==ä==0 && Ωr+Ωm+ΩΛ>1 && Ωr>=0 && Ωm>=0 && ΩΛ>0 && a>1, ΩΛ]] A = f[[1,2]] ΩΛ = f[[2,2]]/.a->A plot[x_]:=Plot[x, {a, A-1, A+1}, Frame -> True, GridLines -> {{A}, {}}] "H²"->plot[H²] "ȧ²"->plot[ȧ²] "ä"->plot[ä]
(* | Evolution of a closed FLRW Universe that reaches an unstable equilibrium | *) set = {"GlobalAdaptive", "MaxErrorIncreases"->100, Method->"GaussKronrodRule"}; (* Integration Rule *) n = 100; (* Recursion Depth *) int[f_, {x_, xmin_, xmax_}] := (* Integral *) NIntegrate[f, {x, xmin, xmax}, Method->set, MaxRecursion->n, WorkingPrecision->wp]; wp = MachinePrecision; (* Working Precision *) im = 200; (* Image Size *) ηmax = 800; pmax = 800; (* Plot Range *) amax = Root[4+11#-10#^2-33#^3+8#^4&, 2, 0]; (* Maximal Scale Factor *) tmax = 444 Gyr; (* Integration Limit *) c = 299792458 m/sek; (* Lightspeed *) G = 667384*^-16 m^3 kg^-1 sek^-2; (* Newton Constant *) Gyr = 10^7*36525*24*3600 sek; (* Billion Years *) Glyr = Gyr*c; (* Billion Lightyears *) Mpc = 30856775777948584200000 m; (* Megaparsec *) ρc[H_] := 3H^2/8/π/G; (* Critical Density *) ρΛ = ρc[H0] ΩΛ; (* Dark Energy Density *) kg = m = sek = 1; (* SI Units *) ΩR = 3/10; (* Radiation Proportion including Neutrinos *) ΩM = 11/10; (* Matter Proportion including Dark Matter *) ΩΛ = 3/32+(77 Root[4+11#-10#^2-33#^3+8#^4&, 2, 0])/16+ (203 Root[4+11#-10#^2-33#^3+8#^4&, 2, 0]^2)/160- (11 Root[4+11#-10#^2-33#^3+8#^4&, 2, 0]^3)/20; (* Dark Energy Proportion *) ΩT = ΩR+ΩM+ΩΛ; (* Total Density over Critical Density *) ΩK = 1-ΩT; (* Curvature Density *) rK = c/H0/Sqrt[-ΩK]; (* Curvature Radius *) H0 = 67150 m/Mpc/sek; (* Hubble Constant *) H[a_] := H0 Sqrt[ΩR/a^4+ΩM/a^3+ΩK/a^2+ΩΛ] (* Hubble Parameter *) sol = Quiet[NDSolve[{A'[t]/A[t] == H[A[t]], A[0] == 1*^-15, WhenEvent[Abs[A[t]] == amax, tMax=t; "StopIntegration"]}, A, {t, 0, tmax}, MaxSteps->∞, WorkingPrecision->wp]]; â[t_] := Evaluate[(A[t]/.sol)[[1]]]; (* Scale Factor a by Time t *) a[t_] := If[t<tmax, â[t], amax]; If[t>tmax, amax, â[t]]; (* Optimized a of t *) т[a_] := int[1/A/H[A], {A, 0, a}]; (* Time t by Scale Factor a *) rP[t_] := a[t] int[c/a[т], {т, 0, t}]; (* Proper Particle Horizon by t *) rp[a_] := a int[c/A^2/H[A], {A, 0, a}]; (* Proper Particle Horizon by a *) RP[t_] := int[c/a[т], {т, 0, t}]; (* Comoving Particle Horizon by t *) Rp[a_] := int[c/A^2/H[A], {A, 0, a}]; (* Comoving Particle Horizon by a *) rE[t_] := Nothing; (* Proper Event Horizon by t *) re[a_] := Nothing; (* Proper Event Horizon by a *) RE[t_] := Nothing; (* Comoving Event Horizon by t *) Rε[a_] := Nothing; (* Comoving Event Horizon by a *) rL[t0_, t_] := a[t] int[c/a[т], {т, t, t0}]; (* Proper Light Cone by t *) rl[a0_, a_] := a int[c/A^2/H[A], {A, a, a0}]; (* Proper Light Cone by a *) RL[t0_, t_] := int[c/a[т], {т, t, t0}]; (* Comoving Light Cone by t *) Rl[a0_, a_] := int[c/A^2/H[A], {A, a, a0}]; (* Comoving Light Cone by a *) rH[t_] := c/H[a[t]]; (* Proper Hubble Radius by t *) rh[a_] := c/H[a]; (* Proper Hubble Radius by a *) RH[t_] := c/H[a[t]]/a[t]; (* Comoving Hubble Radius by t *) Rh[a_] := c/H[a]/a; (* Comoving Hubble Radius by a *) t0 = Quiet[Re[t/.FindRoot[a[t]-1, {t, 10 Gyr}]]]; ti = t Gyr; τi = τ Gyr; "t0"->t0/Gyr "Gyr" (* Current Time *) tmax = Quiet[Re[t/.FindRoot[a[t]-amax, {t, 10 Gyr}]]]; "tmax"->tmax/Gyr "Gyr" ηH = Quiet[Interpolation[Join[{{0, 0}}, (* Hubble Radius by Conformal Time *) ParallelTable[ {Rp[amax (Sin[π a/amax/2])^3]/Glyr, Rh[amax (Sin[π a/amax/2])^3]/Glyr}, {a, amax/2/im, amax-amax/2/im, amax/2/im}], {{amax, Infinity}}, {{2 amax, Infinity}}]]]; rpN = rp[1]/Glyr; "PROPER DISTANCES, f(t)" pt = Quiet[Plot[ {rH[τi]/Glyr, rP[τi]/Glyr, π rK a[τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->pmax/pmax, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot1[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[ {rL[ti, τi]/Glyr, -rL[ti, τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], pt]], 90 Degree]}}]]; Do[Print[plot1[t]], {t, {tmax/Gyr}}] plot2 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[1/amax a[τ Gyr] n^4/250, {n, 4, 55, 1}]], {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(t)" ct = Quiet[Plot[ {RH[τi]/Glyr, RP[τi]/Glyr, π rK/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot3[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[ {RL[ti, τi]/Glyr, -RL[ti, τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], ct]], 90 Degree]}}]]; Do[Print[plot3[t]], {t, {tmax/Gyr}}] plot4 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 100, pmax, 100}]], {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "CONFORMAL DIAGRAM, f(η)" cη = Quiet[Plot[ {ηH[Ct], Ct, π rK/Glyr}, {Ct, 0, ηmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot9[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[ {η-Ct, Ct-η}, {Ct, 0, ηmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], cη]], 90 Degree]}}]]; Do[Print[plot9[η]], {η, {RP[tmax]/Glyr}}] plot10 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 100, pmax, 100}]], {Ct, 0, ηmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] s[text_] := Style[text, FontFamily->"Lucida Console", FontSize->36]