In [229]:
import evolve_background
from numpy.polynomial import Polynomial, Legendre
from scipy.interpolate import interp1d, splrep, splev, InterpolatedUnivariateSpline
In [305]:
Mp = 2.435e27 ##eV
k_B = 8.617332478e-5 ##eV/K
H0 = 1.494e-33 ##70 km/s/Mpc in eV
rho_crit = 3*(H0*Mp)**2
m_chi = H0
m_phi = 1.0/200.0 * H0
z0 =100
params = {
'verbose': False,
'step_density': 5e2,
'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.009) * sqrt(sqrt(3.0)) * (Mp*H0)/(m_phi), ##value of phi field at initial conditions (reheating)?
'chi0': sqrt(0.2) * 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 [306]:
chi = evolve_background.run(**params)
print chi['H'][::-1]
In [296]:
params['chi0'] = 0.0005
params['rho_c0'] = .0035 * sqrt(3) * (Mp * H0)**2 * (1+z0)**3
lcdm = evolve_background.run(**params)
print lcdm['H'][::-1]
In [307]:
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), .0028 * sqrt(3) * (Mp * H0)**2 *(1+z)**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/z0,1)
xlabel('scale factor', fontsize = 20)
ylabel('energy density', fontsize = 20)
Out[307]:
In [308]:
Hlcdm = interp1d(lcdm['z'][::-1], lcdm['H'][::-1])
Hchi = interp1d(chi['z'][::-1], chi['H'][::-1])
DAlcdm = interp1d(lcdm['z'][::-1], lcdm['D_A'][::-1])
DAchi = interp1d(chi['z'][::-1], chi['D_A'][::-1])
In [309]:
figure()
z = linspace(0,z0,100000)
semilogx(z, Hchi(z)/Hlcdm(z) - 1)
#semilogx(z, DAchi(z)/DAlcdm(z) - 1)
ylabel(r'$\delta 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[309]:
In [316]:
z = linspace(0,10,1000)
figure()
semilogx(z, Hlcdm(z), label = r'$\rho_{\chi} = 0$')
semilogx(z, Hchi(z), label = r'$\rho_{\chi} \ne 0$')
legend()
xlabel('redshift', fontsize = 20)
ylabel(r'$H(z)$', fontsize = 20)
Out[316]:
In [311]:
z = linspace(0,z0,100000)
semilogx(z, DAchi(z)/DAlcdm(z) - 1)
ylabel(r'$\delta 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[311]:
In [312]:
z = linspace(0,10,100000)
figure()
loglog(z, DAlcdm(z), label = r'$\rho_{\chi} = 0$')
loglog(z, DAchi(z), label = r'$\rho_{\chi} \ne 0$')
legend()
xlabel('redshift', fontsize = 20)
ylabel(r'$D_A(z)$', fontsize = 20)
Out[312]:
In [313]:
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[313]:
In [314]:
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[314]:
In [50]:
In [ ]: