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_lim=1.2
ah_plot=plot((1/2)*(ah^2-ah^4), (ah,-ah_lim,ah_lim),aspect_ratio=12)
#ah_plot+= circle((1.07,E_ah),0.02,aspect_ratio=0.1,fill=true)
ah_plot+=text(r'$E_{\hat{a}}$', (-1.35,E_ah),fontsize=8)
#ah_plot+=text(r'$E_{\hat{a}}=\frac{1}{2}\epsilon_{\mathrm{w}}^4 E_b$',
#              (-1.6,E_ah),fontsize=8)
#ah_plot+=text(r'$E_{\hat{a}}=\frac{1}{2}\epsilon_{\mathrm{w}}^4 E_b$',
#              (1.4,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'$\frac{1}{\sqrt{2}}$',(1/sqrt(2),-0.015),fontsize=6,aspect_ratio=1)
#ah_plot+=text(r'$-\frac{1}{\sqrt{2}}$',(-1/sqrt(2),-0.015),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'$0$',(-4.5,0),fontsize=2,)
#ah_plot+=text(r'$V_{\hat{a}}$',(-1.5,-0.1),fontsize=14)
ah_plot+=text(r'$V_{\hat{a}}=\frac{1}{2}(\hat{a}^2-\hat{a}^4)$',(0,-0.04),fontsize=12)
#ah_plot+=text(r'$\hat{a}$',(0,-0.05),fontsize=14)
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('ah_plot_2.pdf',dpi=200,axes=False)
show(ah_plot,axes=False,dpi=200,figsize=[3.2,2.4])

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=plot((1/2)*(b^2-1)^2, (b,-3,3),aspect_ratio=.5)
#b_plot+= circle((2.7,E_b),0.25,aspect_ratio=1,fill=true)
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+=arrow((-1.97,E_b-2),(-1.9,E_b-4),width=.5,arrowsize=2)
b_plot.save('b_plot.pdf',dpi=200,axes=False)
show(b_plot,axes=False,dpi=200,figsize=[3.2,2.4])

5.1 Fundamental constants

%display latex
LE = lambda latex_string: LatexExpr(latex_string);
# pretty_print(LE("c ="),c);

declare units as variables

# 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

# fundamental constants (NIST 2018)
# c = 299792458 * m * s^(-1)
# hbar = 1.054571817  * 10^(-34)* J * s
# e = 1.602176634  * 10^(-19) * C 
# G = 6.67430*10^(-11)*m^3*kg^(-1)* s^(-2)
# kB = 1.380649 * 10^(-23) * J * K-1
#
#
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);

$$\newcommand{\Bold}[1]{\mathbf{#1}}c = \frac{299792458 \, \mathit{meters}}{s}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}e = \left(1.60217663400000 \times 10^{-19}\right) \, C$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\hbar = \left(1.05457181700000 \times 10^{-34}\right) \, J s$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}k_{B} = \frac{\left(1.38064900000000 \times 10^{-23}\right) \, J}{K}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}G = \frac{\left(6.67430000000000 \times 10^{-11}\right) \, m^{3}}{\mathit{kg} s^{2}}$$

$$\newcommand{\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

#  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)
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"\hbar ="),hbar);
pretty_print(LE(r"k_{B} ="),kB);
pretty_print(LE(r"G ="),G);
pretty_print(LE(r"\kappa ="),kappa);

$$\newcommand{\Bold}[1]{\mathbf{#1}}\hbar = \left(6.58211956547607 \times 10^{-25}\right) \, \mathit{GeV} s$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}k_{B} = \frac{\left(8.61733326214518 \times 10^{-14}\right) \, \mathit{GeV}}{K}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}G = \frac{\left(4.41583261432942 \times 10^{-63}\right) \, s}{\mathit{GeV}}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\kappa = \frac{\left(1.10981978405276 \times 10^{-61}\right) \, s}{\mathit{GeV}}$$

5.2 Standard Model coupling constants from PDG

# Particle Data Group (2018,2019)
#
# GFermi =  1.1663787*10^(-5) *GeV^(-2);
# mW = 80.379*GeV;
# mH = 125.35*GeV;
# 
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);

$$\newcommand{\Bold}[1]{\mathbf{#1}}G_F = \frac{0.0000116637870000000}{\mathit{GeV}^{2}}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}m_W = 80.3790000000000 \, \mathit{GeV}$$

$$\newcommand{\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)
#pretty_print(LE(r"g^2 = \left(\frac{2 m_W}{v }\right)^2="),g^2);
lambdaH = mH/hbar_v;
pretty_print(LE(r"\lambda = \frac{m_H}{\hbar v }="),lambdaH)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\hbar v = 2^{-1/4} G_F^{-1/2}= 246.219650794137 \, \mathit{GeV}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}v ^{-1} = \left(2.67327142421274 \times 10^{-27}\right) \, s$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}g = \frac{2 m_W}{\hbar v }= 0.652904833068782$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\lambda = \frac{m_H}{\hbar v }= 0.508082923505546$$

5.3 Gravitational and weak time scales $t_{\mathrm{grav}}$, $t_{\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)

$$\newcommand{\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)

$$\newcommand{\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 $\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}"))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{\hbar}\mathcal{E}_{0}= 2.84118595562545 \: t_{\mathrm{W}}^{-4}$$

5.5 Seesaw time scale $t_{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)

$$\newcommand{\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);

$$\newcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac32}\frac{g^2}{\lambda}= 1.02756853595816$$

tI*c;
pretty_print(LE(r"t_{I} c="),tI*c);

$$\newcommand{\Bold}[1]{\mathbf{#1}}t_{I} c= 0.0764308558042792 \, \mathit{meters}$$

5.6 Seesaw ratio $\epsilon_W$

epsilonW = (kappa*hbar/(2*g^2*tI^2))^(1/4);
pretty_print(LE(r"\epsilon_{W}="),epsilonW)

$$\newcommand{\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}}"))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\epsilon_{W}= 1.04068084127939 \:\left(\frac{t_{\mathrm{grav}}}{t_{I}}\right)^{1/2}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\epsilon_{W}= 1.02662576744880 \:\frac{t_{\mathrm{grav}}}{t_{\mathrm{W}}}$$

$$\newcommand{\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)

$$\newcommand{\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)

$$\newcommand{\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)

$$\newcommand{\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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{6\pi^{2}}{g^2}= 138.915667118982$$

7.5 $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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}K(1/\sqrt{2}) = \frac{\Gamma(1/4)^{2}}{4 \pi^1/2}= 1.85407467730137$$

7.6 $\langle \mathrm{cn}^2 \rangle$ for $k=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)

$$\newcommand{\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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}k_B T = 28.8142120659094 \, \mathit{GeV}$$

pretty_print(LE(r"T ="),kT/kB) 

$$\newcommand{\Bold}[1]{\mathbf{#1}}T = 3.34375046077032 \times 10^{14} \, K$$

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

pretty_print(LE(r"\frac{\epsilon_W}{\sqrt{2}} ="),epsilonW/N(sqrt(2))) 

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\epsilon_W}{\sqrt{2}} = 2.39597956199416 \times 10^{-17}$$

9.2 $\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}"))

$$\newcommand{\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}"))

$$\newcommand{\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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\hat t_{\mathrm{EW}} = 0.125799532989201$$

tEW = thatEW * tI;
pretty_print(LE(r"t_{\mathrm{EW}} ="),tEW)

$$\newcommand{\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"))

$$\newcommand{\Bold}[1]{\mathbf{#1}}T_{\mathrm{EW}} = 0.501068706214232 \;\epsilon_a$$