In [89]:
import evolve_background
from numpy.polynomial import Polynomial, Legendre
from scipy.interpolate import interp1d, splrep
In [90]:
Mp = 2.435e27 ##eV
k_B = 8.617332478e-5 ##eV/K
H0 = 1.494e-24 ##70 km/s/Mpc in eV
rho_crit = 3*(H0*Mp)**2
m_chi = 70 * H0
m_phi = 1.0/20.0 * H0
z0 =30
params = {
'verbose': False,
'step_density': 1e3,
'backwards_evolve': False, ##Flag to say whether these are initial conditions today, or far into the past.
'z0':z0, ##redshift of your initial conditions (reheating)?? Shouldn't matter much as long as it's far into radiation era, and t(z0) << period of DM oscillator
'z_end':0, ##Redshift where you end the integration.
'rho_b0':.0002 * sqrt(3) * (Mp*H0)**2 * (1+z0)**3,##Baryons. Some thought goes into what these should be. What we know is that \Omega_b h^2 = 0.022, but we don't know what h = H_0/100 is in this model.
##Should be enough information, though, to work it out in terms of the other numbers here...not quite seeing how right now though.
'rho_c0':.0028 * sqrt(3) * (Mp * H0)**2 * (1+z0)**3, ##Dark Matter (non oscillating). Get this from \Omega_c h^2 = 0.120
'rho_nu0':7.0/8.0 * 3.046 * (4.0/11.0)**(4/3) * pi**2 * (k_B * 2.73)**4/15.0 * (1+z0)**4, ##Neutrinos. Get this from \Omega_{\nu} h^2 = \frac{7}{8} N_{\nu} \left(\frac{4}{11}\right)^{4/3} \Omega_{\gamma} h^2, where N_{\nu} = 3.046 including 1 loop.
'rho_gamma0':pi**2 * (k_B * 2.73)**4/15.0 * (1+z0)**4,##This is just the energy density of a black body spectrum with temperature T_{CMB}. Easy to figure out, but will depend on units,
##and where you're setting your initial conditions (when is reheating?).
'V_phi':Polynomial([0, 0, m_phi**2, 0]), ##Can add more terms in the effective field theory here. Polynomial expands in powers of phi.
'V_chi':Polynomial([0, 0, m_chi**2, 0]), ##Also can add here. For true cosine potentials, switch to Legendre expansion
'phi0': sqrt(0.010) * sqrt(sqrt(3.0)) * (Mp*H0)/(m_phi), ##value of phi field at initial conditions (reheating)?
'chi0': sqrt(0.1) * sqrt(sqrt(3.0)) * (Mp*H0)/(m_chi) ##value of chi field. Code assumes initial velocities are given by slow roll; chidot0 = -V'/3H, etc.
}
print 'chi0 is ', params['chi0']
print 'phi0 is ', params['phi0']
rhoc = params['rho_c0']
rhophi = params['V_phi'](params['phi0'])
rhochi = params['V_chi'](params['chi0'])
print 'densities are:\n\t matter: %(rhoc)g \n\t phi: %(rhophi)g \n\t chi: %(rhochi)g \n\t at z = %(z0)g'%{'rhoc':rhoc,'rhophi':rhophi,'rhochi':rhochi,'z0':z0}
In [91]:
chi = evolve_background.run(**params)
In [92]:
params['chi0'] = 0.0005
params['phi0'] = sqrt(0.007) * sqrt(sqrt(3.0)) * (Mp*H0)/(m_phi)
lcdm = evolve_background.run(**params)
In [93]:
loglog(1/(1.0+chi['z']), chi['Vchi'] + 0.5 * chi['chidot']**2,label = r'$\rho_{\chi}$')
loglog(1/(1.0+chi['z']), chi['Vphi'] + 0.5 * chi['phidot']**2,label = r'$\rho_{\phi}$')
loglog(1/(1.0 + lcdm['z']), lcdm['Vphi'] + 0.05 * lcdm['phidot']**2, linestyle = 'dashed', color = 'green')
z = chi['z']
loglog(1/(1.0+z), params['rho_c0']*(1+z)**3/(1+z0)**3 + chi['Vchi'] + 0.5 * chi['chidot']**2, label = r'$\rho_m$')
loglog(1/(1.0+z), params['rho_c0']*(1+z)**3/(1+z0)**3, linestyle = 'dashed', color = 'red')
legend()
xlim(.1,1)
xlabel('scale factor', fontsize = 20)
ylabel('energy density', fontsize = 20)
Out[93]:
In [107]:
Hlcdm = interp1d(lcdm['z'], lcdm['H'])
Hchi = interp1d(chi['z'], chi['H'])
D_Alcdm = interp1d(lcdm['z'], lcdm['D_A'], s=0)
D_Achi = interp1d(chi['z'], chi['D_A'], s=0)
plot(lcdm['z'],D_Alcdm(lcdm['z']))
In [82]:
z = linspace(0,z0,100000)
semilogx(z, Hchi(z)/Hlcdm(z) - 1)
ylabel(r'$H(z)/H(z)^{\rho_{\chi} = 0}$', fontsize = 20)
#ylim(5e-6,8e-6)
gca().yaxis.set_major_formatter(FormatStrFormatter('%g'))
xlabel('redshift', fontsize = 20)
Out[82]:
In [88]:
z = linspace(0, z0, 10000)
loglog(z, Hchi(z)-Hlcdm(z), color = 'green', linestyle = 'dashed')
loglog(z, Hlcdm(z))
loglog(z, Hchi(z))
Out[88]:
In [74]:
z = linspace(0,z0,100000)
semilogx(z, D_Achi(z)/D_Alcdm(z) - 1)
ylabel(r'$D_A/D_A^{\rho_{\chi} = 0}$', fontsize = 20)
#ylim(5e-6,8e-6)
gca().yaxis.set_major_formatter(FormatStrFormatter('%g'))
xlabel('redshift', fontsize = 20)
Out[74]:
In [75]:
semilogx(lcdm['z'], lcdm['wtot'],label = r'$\chi = 0$' )
semilogx(chi['z'], chi['wtot'],label = r'$\chi \ne 0$')
semilogx(chi['z'], chi['wchi'], linestyle = 'dashed', label = r'$w_{\chi}$')
legend(loc = 2)
ylim(-1.1, 1.1)
xlabel('redshift', fontsize = 20)
ylabel(r'$w_{\rho + \chi}$', fontsize = 20)
Out[75]:
In [76]:
output = lcdm
global_vars = params
output['wdark'] = ((-output['Vchi'] + output['chidot']**2/2 -output['Vphi'] + output['phidot']**2/2) /
(output['chidot']**2/2 + output['Vchi'] + output['phidot']**2/2 + output['Vphi'] + output['chidot']**2/2 + global_vars['rho_c0'] * (1+output['z'])**3/(1+z0)**3))
semilogx(output['z'], output['wdark']+1, label = r'$\rho_{\chi} = 0$')
output = chi
output['wdark'] = ((-output['Vchi'] + output['chidot']**2/2 -output['Vphi'] + output['phidot']**2/2) /
(output['chidot']**2/2 + output['Vchi'] + output['phidot']**2/2 + output['Vphi'] + output['chidot']**2/2 + global_vars['rho_c0'] * (1+output['z'])**3/(1+z0)**3))
semilogx(output['z'], output['wdark']+1, label = r'$\rho_{\chi} \ne 0$')
legend()
xlabel('redshift', fontsize = 20)
ylabel(r'$w_{\chi + \phi + DM}+1$', fontsize = 20)
Out[76]:
In [ ]: