CAMB: trouble with reproducing the result of 1507.00982

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
Jiwon Park
Posts: 9
Joined: July 09 2019
Affiliation: Soongsil University

CAMB: trouble with reproducing the result of arxiv:1507.00982

Post by Jiwon Park » May 20 2021

Dear all,

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
    
    (......)
    
In the subroutine for scalar initial conditions, I added the corresponding conditions for new components.

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)
    
And to solve new equations, the following were added in the subroutine 'derivs.'

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
        
    (......)

Post Reply