Numerical calculations for _Origin of Cosmological Temperature_.ipynb

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])

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])

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);

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);

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);

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)

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)

tW = hbar/mW;
pretty_print(LE(r"t_{\mathrm{W}}=\frac{\hbar}{m_W}="),tW)

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}"))

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)

ratio12 = N(sqrt(3/2))*g^2/lambdaH;
pretty_print(LE(r"\sqrt{\frac32}\frac{g^2}{\lambda}="),ratio12);

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

5.6 Seesaw ratio $\epsilon_W$

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

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}}"))

ratio134=(2*g^2)^(-1/4);
pretty_print(LE(r"\left(\frac{1}{2g^2}\right)^{1/4}="),ratio134)

ratio13=(lambdaH^2/(3*g^6))^(1/4);
pretty_print(LE(r"\left(\frac{\lambda^2}{3g^6}\right)^{1/4}="),ratio13)

ratio14=((3*g^2)/(4*lambdaH^2))^(1/4);
pretty_print(LE(r"\left(\frac{3g^2}{4\lambda^2}\right)^{1/4}="),ratio14)

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)

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)

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)

8. Cosmological temperature

kT = mH/N((6*pi)^(1/2));
pretty_print(LE(r"k_B T ="),kT)

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

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

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

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}"))

pretty_print(LE(r"\hat a_{\mathrm{EW}} ="),\
             sqrt(ratio6),\
            LE(r"\:(2 E_{\hat a})^{1/4}"))

thatEW = asinh(ratio6)/2;
pretty_print(LE(r"\hat t_{\mathrm{EW}} ="),thatEW)

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

TEW = inverse_jacobi('cn', e^(-thatEW), 0.5);
pretty_print(LE(r"T_{\mathrm{EW}} ="),TEW,LE(r"\;\epsilon_a"))