Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

Supplemental material for "Origins of cosmological temperature"

Views: 349
Kernel: SageMath (system-wide)

Numerical calculations for "Origins of Cosmological Temperature"

This SageMath notebook performs numerical calculations for the paper Origins of Cosmological Temperature and the supplemental note Calculations for "Origins of Cosmological Temperature".

It makes graphs of the two anharmonic potentials and does some arithmetic.

Section headings are as in the supplemental note.

3.5 Action in dimensionless variables

ah = var('ah') E_ah = 0.15 ah_lim=1.08 ah_plot=plot((1/2)*(ah^2-ah^4), (ah,-ah_lim,ah_lim),aspect_ratio=12) ah_plot+=text(r'$E_{\hat{a}}$', (-1.35,E_ah),fontsize=8) ah_plot+=text(r'$0$',(-1.3,0),fontsize=6) ah_plot+=text(r'$\frac{1}{8}$',(-1.32,0.125),fontsize=6,aspect_ratio=1) ah_plot+=text(r'$0$',(0,-0.006),fontsize=6,aspect_ratio=1) ah_plot+=text(r'$1$',(0.95,-0.006),fontsize=6,aspect_ratio=1) ah_plot+=text(r'$-1$',(-1.0,-0.006),fontsize=6,aspect_ratio=1) ah_plot+=text(r'$V_{\hat{a}}=\frac{1}{2}(\hat{a}^2-\hat{a}^4)$',(0,-0.04),fontsize=12) ah_plot+= plot(0.125,(ah,-1.25,1.25),linestyle=":") ah_plot+= plot(0,(ah,-1.25,1.25),linestyle=":") ah_plot+= plot(E_ah,(ah,-1.25,1.25),linestyle=":") ah_plot+=arrow((1.11,-0.08),(1.13,-0.08-.02),width=.5,arrowsize=1.5) ah_plot+=arrow((-1.13,-0.08-.02),(-1.11,-0.08),width=.5,arrowsize=1.5) ah_plot+=arrow((0.97,1/16),(1.00,1/16-.02),width=.5,arrowsize=1.5) ah_plot+=arrow((-1.00,1/16-.02),(-0.97,1/16),width=.5,arrowsize=1.5) ah_plot+=arrow((-.34,1/16),(-.27,1/16-.017),width=.5,arrowsize=1.5) ah_plot+=arrow((.27,1/16-.017),(.34,1/16),width=.5,arrowsize=1.5) ah_plot.save('plot_ah.pdf',dpi=200,axes=False) show(ah_plot,axes=False,dpi=200,figsize=[3.2,2.4])
Image in a Jupyter notebook
b = var('b') E_b = 12 b_plot=plot((1/2)*(b^2-1)^2, (b,-2.5,2.5),aspect_ratio=.5) b_plot+=text(r'$E_b$',(-3.1,E_b),fontsize=10) b_plot+=text(r'$0$',(-3.0,0),fontsize=6) b_plot+=text(r'$V_b=\frac{1}{2}(b^2-1)^2$',(0,E_b/2),fontsize=12) b_plot+=text(r'$b$',(0,-1),fontsize=12) b_plot+= plot(0,(b,-2.7,2.7),linestyle=":") b_plot+= plot(E_b,(b,-2.7,2.7),linestyle=":") b_plot+=text(r'$1$',(1.0,-0.4),fontsize=6) b_plot+=text(r'$0$',(0.0,-0.4),fontsize=6) b_plot+=text(r'$-1$',(-1.09,-0.4),fontsize=6) b_plot+=arrow((-1.97,E_b-2),(-1.9,E_b-4),width=.5,arrowsize=2) b_plot.save('plot_b.pdf',dpi=200,axes=False) show(b_plot,axes=False,dpi=200,figsize=[3.2,2.4])
Image in a Jupyter notebook

5.1 Fundamental constants

%display latex LE = lambda latex_string: LatexExpr(latex_string);

declare units as variables

s = var('s', domain='positive'); assume(s,'real'); GeV = var('GeV', domain='positive'); assume(GeV,'real'); J = var('J', domain='positive'); assume(J,'real'); m = var('m', domain='positive'); assume(m,'real'); meters = var('meters', domain='positive'); assume(meters,'real'); kg = var('kg', domain='positive'); assume(kg,'real'); K = var('K', domain='positive'); assume(K,'real'); C = var('C', domain='positive'); assume(C,'real');

fundamental constents from NIST 2018

c = 299792458 * meters * s^(-1); e_charge = 1.602176634 * 10^(-19) * C; hbar = 1.054571817*10^(-34)*J*s; kB = 1.380649*10^(-23)*J*K^(-1); G = 6.67430*10^(-11)*m^3*kg^(-1)* s^(-2); kappa = N(8*pi)*G; # pretty_print(LE(r"c ="),c); pretty_print(LE(r"e ="),e_charge); pretty_print(LE(r"\hbar ="),hbar); pretty_print(LE(r"k_{B} ="),kB); pretty_print(LE(r"G ="),G); pretty_print(LE(r"\kappa ="),kappa);
c=299792458meterss\renewcommand{\Bold}[1]{\mathbf{#1}}c = \frac{299792458 \, \mathit{meters}}{s}
e=(1.60217663400000×1019)C\renewcommand{\Bold}[1]{\mathbf{#1}}e = \left(1.60217663400000 \times 10^{-19}\right) \, C
=(1.05457181700000×1034)Js\renewcommand{\Bold}[1]{\mathbf{#1}}\hbar = \left(1.05457181700000 \times 10^{-34}\right) \, J s
kB=(1.38064900000000×1023)JK\renewcommand{\Bold}[1]{\mathbf{#1}}k_{B} = \frac{\left(1.38064900000000 \times 10^{-23}\right) \, J}{K}
G=(6.67430000000000×1011)m3kgs2\renewcommand{\Bold}[1]{\mathbf{#1}}G = \frac{\left(6.67430000000000 \times 10^{-11}\right) \, m^{3}}{\mathit{kg} s^{2}}
κ=(1.67743454782835×109)m3kgs2\renewcommand{\Bold}[1]{\mathbf{#1}}\kappa = \frac{\left(1.67743454782835 \times 10^{-9}\right) \, m^{3}}{\mathit{kg} s^{2}}

use c=1 units with unit of energy = GeV

m = c^(-1)*meters; J = e_charge^(-1) * C * 10^(-9) * GeV kg = J*s^2*m^(-2) def conv(*args): return [arg.subs(m=m,kg=kg,J=J) for arg in args] [hbar,kB,G,kappa] = conv(hbar,kB,G,kappa) pretty_print(LE(r"\hbar ="),hbar); pretty_print(LE(r"k_{B} ="),kB); pretty_print(LE(r"G ="),G); pretty_print(LE(r"\kappa ="),kappa);
=(6.58211956547607×1025)GeVs\renewcommand{\Bold}[1]{\mathbf{#1}}\hbar = \left(6.58211956547607 \times 10^{-25}\right) \, \mathit{GeV} s
kB=(8.61733326214518×1014)GeVK\renewcommand{\Bold}[1]{\mathbf{#1}}k_{B} = \frac{\left(8.61733326214518 \times 10^{-14}\right) \, \mathit{GeV}}{K}
G=(4.41583261432942×1063)sGeV\renewcommand{\Bold}[1]{\mathbf{#1}}G = \frac{\left(4.41583261432942 \times 10^{-63}\right) \, s}{\mathit{GeV}}
κ=(1.10981978405276×1061)sGeV\renewcommand{\Bold}[1]{\mathbf{#1}}\kappa = \frac{\left(1.10981978405276 \times 10^{-61}\right) \, s}{\mathit{GeV}}

5.2 Standard Model coupling constants from PDG (2018, 2019)

GFermi = 1.1663787*10^(-5)*GeV^(-2); mW = 80.379*GeV; mH = 125.10*GeV; # pretty_print(LE(r"G_F ="),GFermi); pretty_print(LE(r"m_W ="),mW); pretty_print(LE(r"m_H ="),mH);
GF=0.0000116637870000000GeV2\renewcommand{\Bold}[1]{\mathbf{#1}}G_F = \frac{0.0000116637870000000}{\mathit{GeV}^{2}}
mW=80.3790000000000GeV\renewcommand{\Bold}[1]{\mathbf{#1}}m_W = 80.3790000000000 \, \mathit{GeV}
mH=125.100000000000GeV\renewcommand{\Bold}[1]{\mathbf{#1}}m_H = 125.100000000000 \, \mathit{GeV}
hbar_v = N(2^(-1/4))*GFermi^(-1/2) pretty_print(LE(r"\hbar v = 2^{-1/4} G_F^{-1/2}="),hbar_v) v = hbar_v/hbar pretty_print(LE(r"v ^{-1} ="),1/v); g = 2*mW/hbar_v; pretty_print(LE(r"g = \frac{2 m_W}{\hbar v }="),g) lambdaH = mH/hbar_v; pretty_print(LE(r"\lambda = \frac{m_H}{\hbar v }="),lambdaH)
v=21/4GF1/2=246.219650794137GeV\renewcommand{\Bold}[1]{\mathbf{#1}}\hbar v = 2^{-1/4} G_F^{-1/2}= 246.219650794137 \, \mathit{GeV}
v1=(2.67327142421274×1027)s\renewcommand{\Bold}[1]{\mathbf{#1}}v ^{-1} = \left(2.67327142421274 \times 10^{-27}\right) \, s
g=2mWv=0.652904833068782\renewcommand{\Bold}[1]{\mathbf{#1}}g = \frac{2 m_W}{\hbar v }= 0.652904833068782
λ=mHv=0.508082923505546\renewcommand{\Bold}[1]{\mathbf{#1}}\lambda = \frac{m_H}{\hbar v }= 0.508082923505546

5.3 Gravitational and weak time scales tgravt_{\mathrm{grav}}, tWt_{\mathrm{W}}

tgrav = sqrt(kappa*hbar); pretty_print(LE(r"t_{\mathrm{grav}} = (\hbar\kappa)^{1/2}=\ (8\pi)^{1/2} t_{P}="),N((8*pi)^(1/2)),LE(r"t_{P}\\="),tgrav)
tgrav=(κ)1/2=(8π)1/2tP=5.01325654926200tP=(2.70277015574135×1043)s\renewcommand{\Bold}[1]{\mathbf{#1}}t_{\mathrm{grav}} = (\hbar\kappa)^{1/2}=(8\pi)^{1/2} t_{P}= 5.01325654926200 t_{P}\\= \left(2.70277015574135 \times 10^{-43}\right) \, s
tW = hbar/mW; pretty_print(LE(r"t_{\mathrm{W}}=\frac{\hbar}{m_W}="),tW)
tW=mW=(8.18885475743176×1027)s\renewcommand{\Bold}[1]{\mathbf{#1}}t_{\mathrm{W}}=\frac{\hbar}{m_W}= \left(8.18885475743176 \times 10^{-27}\right) \, s

5.4 The scalar field energy density E0\mathcal{E}_{0}

E0 = hbar*lambdaH^2*v^4/8; ratio1 = tW^4*E0/hbar; pretty_print(LE(r"\frac{1}{\hbar}\mathcal{E}_{0}="),\ ratio1,LE(r"\: t_{\mathrm{W}}^{-4}"))
1E0=2.84118595562545tW4\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{\hbar}\mathcal{E}_{0}= 2.84118595562545 \: t_{\mathrm{W}}^{-4}

5.5 Seesaw time scale tIt_{I}

tI = (3/(kappa*E0))^(1/2); ratio2 = tI*tgrav/tW^2; pretty_print(LE(r"t_{I}="), ratio2,LE(r"\:\frac{t_{\mathrm{W}}^{2}}{t_{\mathrm{grav}}}"),\ LE(r"="),tI)
tI=1.02756853595816tW2tgrav=(2.54945892615748×1010)s\renewcommand{\Bold}[1]{\mathbf{#1}}t_{I}= 1.02756853595816 \:\frac{t_{\mathrm{W}}^{2}}{t_{\mathrm{grav}}} = \left(2.54945892615748 \times 10^{-10}\right) \, s
ratio12 = N(sqrt(3/2))*g^2/lambdaH; pretty_print(LE(r"\sqrt{\frac32}\frac{g^2}{\lambda}="),ratio12);
32g2λ=1.02756853595816\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac32}\frac{g^2}{\lambda}= 1.02756853595816
tI*c; pretty_print(LE(r"t_{I} c="),tI*c);
tIc=0.0764308558042792meters\renewcommand{\Bold}[1]{\mathbf{#1}}t_{I} c= 0.0764308558042792 \, \mathit{meters}

5.6 Seesaw ratio ϵW\epsilon_W

epsilonW = (kappa*hbar/(2*g^2*tI^2))^(1/4); pretty_print(LE(r"\epsilon_{W}="),epsilonW)
ϵW=3.38842679174089×1017\renewcommand{\Bold}[1]{\mathbf{#1}}\epsilon_{W}= 3.38842679174089 \times 10^{-17}
ratio3 = epsilonW*tW/tgrav; ratio4 = epsilonW*tI/tW; ratio34 = sqrt(ratio3*ratio4); pretty_print(LE(r"\epsilon_{W}="),ratio34, LE(r"\:\left(\frac{t_{\mathrm{grav}}}{t_{I}}\right)^{1/2}")) pretty_print(LE(r"\epsilon_{W}="),ratio3, LE(r"\:\frac{t_{\mathrm{grav}}}{t_{\mathrm{W}}}")) pretty_print(LE(r"\epsilon_{W}="),ratio4, LE(r"\:\frac{t_{\mathrm{W}}}{t_{I}}"))
ϵW=1.04068084127939(tgravtI)1/2\renewcommand{\Bold}[1]{\mathbf{#1}}\epsilon_{W}= 1.04068084127939 \:\left(\frac{t_{\mathrm{grav}}}{t_{I}}\right)^{1/2}
ϵW=1.02662576744880tgravtW\renewcommand{\Bold}[1]{\mathbf{#1}}\epsilon_{W}= 1.02662576744880 \:\frac{t_{\mathrm{grav}}}{t_{\mathrm{W}}}
ϵW=1.05492833683428tWtI\renewcommand{\Bold}[1]{\mathbf{#1}}\epsilon_{W}= 1.05492833683428 \:\frac{t_{\mathrm{W}}}{t_{I}}
ratio134=(2*g^2)^(-1/4); pretty_print(LE(r"\left(\frac{1}{2g^2}\right)^{1/4}="),ratio134)
(12g2)1/4=1.04068084127939\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{1}{2g^2}\right)^{1/4}= 1.04068084127939
ratio13=(lambdaH^2/(3*g^6))^(1/4); pretty_print(LE(r"\left(\frac{\lambda^2}{3g^6}\right)^{1/4}="),ratio13)
(λ23g6)1/4=1.02662576744880\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{\lambda^2}{3g^6}\right)^{1/4}= 1.02662576744880
ratio14=((3*g^2)/(4*lambdaH^2))^(1/4); pretty_print(LE(r"\left(\frac{3g^2}{4\lambda^2}\right)^{1/4}="),ratio14)
(3g24λ2)1/4=1.05492833683428\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{3g^2}{4\lambda^2}\right)^{1/4}= 1.05492833683428

5.7 Units of action for the two oscillators

ratio5=6*N(pi)^2/g^2; pretty_print(LE(r"\frac{6\pi^{2}}{g^2}="),ratio5)
6π2g2=138.915667118982\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{6\pi^{2}}{g^2}= 138.915667118982

7.5 K(1/2)K(1/\sqrt{2})

K = N(gamma(1/4)^2/(4*pi^(1/2))); pretty_print(LE(r"K(1/\sqrt{2}) = \frac{\Gamma(1/4)^{2}}{4 \pi^{1/2}}="),K)
K(1/2)=Γ(1/4)24π1/2=1.85407467730137\renewcommand{\Bold}[1]{\mathbf{#1}}K(1/\sqrt{2}) = \frac{\Gamma(1/4)^{2}}{4 \pi^{1/2}}= 1.85407467730137

7.6 cn2\langle \mathrm{cn}^2 \rangle for k=1/2k=1/\sqrt2

cn2ave = N(pi/2)/K^2; pretty_print(LE(r"\langle \mathrm{cn}^2 \rangle=\frac{2}{\pi^{2}}\frac{1}{K^{2}}="),cn2ave)
cn2=2π21K2=0.456946581044464\renewcommand{\Bold}[1]{\mathbf{#1}}\langle \mathrm{cn}^2 \rangle=\frac{2}{\pi^{2}}\frac{1}{K^{2}}= 0.456946581044464

8. Cosmological temperature

kT = mH/N((6*pi)^(1/2)); pretty_print(LE(r"k_B T ="),kT)
kBT=28.8142120659094GeV\renewcommand{\Bold}[1]{\mathbf{#1}}k_B T = 28.8142120659094 \, \mathit{GeV}
pretty_print(LE(r"T ="),kT/kB)
T=3.34375046077032×1014K\renewcommand{\Bold}[1]{\mathbf{#1}}T = 3.34375046077032 \times 10^{14} \, K

9.1 Solution for a^\hat a in co-moving time

pretty_print(LE(r"\frac{\epsilon_W}{\sqrt{2}} ="),epsilonW/N(sqrt(2)))
ϵW2=2.39597956199416×1017\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\epsilon_W}{\sqrt{2}} = 2.39597956199416 \times 10^{-17}

9.2 a^EW\hat a_{\mathrm{EW}}

ratio6 =N(3^(1/2)*pi/(8*K^2))*(2*mW/mH); pretty_print(LE(r"\hat a^2_{\mathrm{EW}} ="),\ LE(r"\frac{3^{1/2}\pi}{8 K^2}\frac{2m_W}{m_H}(2 E_{\hat a})^{1/2}="),\ ratio6,LE(r"\:(2 E_{\hat a})^{1/2}"))
a^EW2=31/2π8K22mWmH(2Ea^)1/2=0.254261938075174(2Ea^)1/2\renewcommand{\Bold}[1]{\mathbf{#1}}\hat a^2_{\mathrm{EW}} = \frac{3^{1/2}\pi}{8 K^2}\frac{2m_W}{m_H}(2 E_{\hat a})^{1/2}= 0.254261938075174 \:(2 E_{\hat a})^{1/2}
pretty_print(LE(r"\hat a_{\mathrm{EW}} ="),\ sqrt(ratio6),\ LE(r"\:(2 E_{\hat a})^{1/4}"))
a^EW=0.504243927157457(2Ea^)1/4\renewcommand{\Bold}[1]{\mathbf{#1}}\hat a_{\mathrm{EW}} = 0.504243927157457 \:(2 E_{\hat a})^{1/4}
thatEW = asinh(ratio6)/2; pretty_print(LE(r"\hat t_{\mathrm{EW}} ="),thatEW)
t^EW=0.125799532989201\renewcommand{\Bold}[1]{\mathbf{#1}}\hat t_{\mathrm{EW}} = 0.125799532989201
tEW = thatEW * tI; pretty_print(LE(r"t_{\mathrm{EW}} ="),tEW)
tEW=(3.20720742285761×1011)s\renewcommand{\Bold}[1]{\mathbf{#1}}t_{\mathrm{EW}} = \left(3.20720742285761 \times 10^{-11}\right) \, s
TEW = inverse_jacobi('cn', e^(-thatEW), 0.5); pretty_print(LE(r"T_{\mathrm{EW}} ="),TEW,LE(r"\;\epsilon_a"))
TEW=0.501068706214232  ϵa\renewcommand{\Bold}[1]{\mathbf{#1}}T_{\mathrm{EW}} = 0.501068706214232 \;\epsilon_a