I have some troubles with reproducing the result of 1507.00982, in which the effect of DM velocity dispersion on CMB is taken into account.
I modified the file 'equation.f90' corresponding to the paper but cannot reproduce the same result in the code.
Could you please explain where I made a mistake?
I guess I might have messed up the indexing or arranged the equations in the wrong order, but I'm not sure.
I checked the result recovers the standard LCDM model when \beta^2=0 (in my code 'beta_dm'), where the paper says that CMB should be the same as LCDM in this case.
I'm using the most recent CAMB version.
Here is how I modified 'equation.f90'.
First, I added indexes for new components.
Code: Select all
(......)
integer, parameter :: basic_num_eqns = 7 !originally 4(add 3 eqs)
integer, parameter :: ix_etak=1, ix_clxc=2, ix_clxb=3, ix_vb=4 !Scalar array indices for each quantity
integer, parameter :: ix_v_dm=5, ix_v1_dm=6, ix_Pi_dm=7 !for DM velocity & its dispersion perturbation
(......)
Code: Select all
(......)
integer, parameter :: i_clxg=1,i_clxr=2,i_clxc=3, i_clxb=4, &
i_qg=5,i_qr=6,i_vb=7,i_pir=8, i_eta=9, i_aj3r=10,i_clxde=11,i_vde=12, &
! additional initial conditions
i_v_dm=13, i_v1_dm=14, i_Pi_dm=15
(......)
! external vars
real(dl) sigma_0_dm, hplus6eta
(......)
!relationship between sigma_0^2 and beta^2
sigma_0_dm=CP%beta_dm/a2
! Set adiabatic initial conditions
chi=1 !Get transfer function for chi
initv(1,i_clxg)=-chi*EV%Kf(1)/3*x2*(1-omtau/5)
initv(1,i_clxr)= initv(1,i_clxg)
initv(1,i_clxb)=0.75_dl*initv(1,i_clxg)
! here we have modified conditions (eqs 2.62-2.65)
! this was modified
initv(1,i_clxc)=(-initv(1,i_v1_dm)+initv(1,i_clxb)*(1+5*sigma_0_dm/6))/(1-0.5*sigma_0_dm)
!the original was: initv(1,i_clxc)=initv(1,i_clxb)
initv(1,i_qg)=initv(1,i_clxg)*x/9._dl
initv(1,i_qr)=-chi*EV%Kf(1)*(4*Rv+23)/Rp15*x3/27
initv(1,i_vb)=0.75_dl*initv(1,i_qg)
initv(1,i_pir)=chi*4._dl/3*x2/Rp15*(1+omtau/4*(4*Rv-5)/(2*Rv+15))
initv(1,i_aj3r)=chi*4/21._dl/Rp15*x3
initv(1,i_eta)=-chi*2*EV%Kf(1)*(1 - x2/12*(-10._dl/Rp15 + EV%Kf(1)))
! this is auxillary
hplus6eta=-2*initv(1,i_clxb)+6*(-chi*2*EV%Kf(1)*(1 - x2/12*(-10._dl/Rp15 + EV%Kf(1))))
! this were added
initv(1,i_v_dm)=0
initv(1,i_v1_dm)=5*sigma_0_dm*initv(1,i_clxb)/9
initv(1,i_Pi_dm)=-2*sigma_0_dm*hplus6eta/27
(......)
y(ix_etak)= -InitVec(i_eta)*k/2
!get eta_s*k, where eta_s is synchronous gauge variable
! CDM
y(ix_clxc)=InitVec(i_clxc)
! Baryons
y(ix_clxb)=InitVec(i_clxb)
y(ix_vb)=InitVec(i_vb)
! DM velocity added term
y(ix_v_dm)=InitVec(i_v_dm)
y(ix_v1_dm)=InitVec(i_v1_dm)
y(ix_Pi_dm)=InitVec(i_Pi_dm)
Code: Select all
(......)
real(dl) dgrho_de, dgq_de, cs2_de
!vars for nonzero DM velocity dispersion based on 1507.00982
!nonzero DM velocity
real(dl) v_dm, v_dm_dot
!DM velocity dispersion sigma_0^2 defined by eq 3.2 on the paper
real(dl) sigma_0_dm
!DM velocity dispersion perturbation v_1^2 and Pi (eq 2.16&2.26)
real(dl) v1_dm, v1_dm_dot, Pi_dm, Pi_dm_dot
k=EV%k_buf
k2=EV%k2_buf
(......)
!relationship between sigma_0^2 and beta^2
sigma_0_dm=CP%beta_dm/a2
etak=ay(ix_etak)
! CDM variables
clxc=ay(ix_clxc)
! Baryon variables
clxb=ay(ix_clxb)
vb=ay(ix_vb)
! DM velocity dispersion vars
v1_dm=ay(ix_v1_dm)
Pi_dm=ay(ix_Pi_dm)
! DM velocity var
v_dm=ay(ix_v_dm)[/b]
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=State%grhob/a
grhoc_t=State%grhoc*(1+0.5*sigma_0_dm)/a
grhor_t=State%grhornomass/a2
grhog_t=State%grhog/a2
(......)
! 8*pi*a*a*SUM[(rho_i+p_i)*v_i] ; added eq 2.44
dgq=grhob_t*vb + grhoc_t*v_dm/(1+sigma_0_dm/2)
!DM pressure term added
gpres_noDE = gpres_nu + (grhor_t + grhog_t)/3 + sigma_0_dm*grhoc_t/(3+1.5*sigma_0_dm)
(......)
if (.not. EV%is_cosmological_constant) &
call State%CP%DarkEnergy%PerturbationEvolve(ayprime, w_dark_energy_t, &
EV%w_ix, a, adotoa, k, z, ay)
! CDM equation of motion
!original was: clxcdot=-k*z
clxcdot=-(adotoa*(v1_dm-sigma_0_dm*clxc)+k*v_dm+k*z*(1+5*sigma_0_dm/6))/(1+0.5*sigma_0_dm)
ayprime(ix_clxc)=clxcdot
! Baryon equation of motion.
clxbdot=-k*(z+vb)
(......)
else
vbdot=-adotoa*vb+k*delta_p_b-photbar*opacity*(4._dl/3*vb-qg)
end if
ayprime(ix_vb)=vbdot
!DM velocity eq 2.33
!original was: v_dm_dot=0
v_dm_dot=-adotoa*v_dm+k*v1_dm/3+k*Pi_dm
ayprime(ix_v_dm)=v_dm_dot
! DM velocity dispersion eq 2.34&2.35
! pls note that h_dot=2kz & h_dot+6eta_dot=2k*sigma
v1_dm_dot=-2*adotoa*v1_dm-5*k*z*sigma_0_dm/3
ayprime(ix_v1_dm)=v1_dm_dot
Pi_dm_dot=-2*adotoa*Pi_dm-4*k*sigma*sigma_0_dm/9
ayprime(ix_Pi_dm)=Pi_dm_dot
if (.not. EV%no_phot_multpoles) then
! Photon equations of motion
(......)
if (EV%is_cosmological_constant) then
dgrho_de=0
dgq_de=0
end if
! eq 2.45&46
dgpi = grhor_t*pir + grhog_t*pig + grhoc_t*Pi_dm/(1+0.5*sigma_0_dm)
dgpi_diff = grhoc_t*Pi_dm*(0.5*sigma_0_dm-1)/(0.5*sigma_0_dm+1) !sum (3*p_nu -rho_nu)*pi_nu
pidot_sum = grhog_t*pigdot + grhor_t*pirdot + grhoc_t*Pi_dm_dot/(1+0.5*sigma_0_dm)
clxnu =0
if (State%CP%Num_Nu_Massive /= 0) then
call MassiveNuVarsOut(EV,ay,ayprime,a, adotoa, dgpi=dgpi, clxnu_all=clxnu, &
dgpi_diff=dgpi_diff, pidot_sum=pidot_sum)
end if
(......)