[A] Changes in module mgvariables in equations_ppf.f90

real(dl) :: fR0, fRn ! the |f_{R0}| and n parameters in the Hu-Sawicki model

[B1] Changes in subroutine derivs in equations_ppf.f90

  Add the following just after ‘use mgvariables’, which is at the very beginning of this subroutine.

use constants, only : c

  In the parameter definition part (also at the beginning of this subroutine), add the following,

real(dl) mubar, Rbar, R0bar, Rdot, mudot

  Just after ‘if (tempmodel /= 0) then’, change the next line to,

if (model==1 .or.model==4 .or.model==5.or.model==6 .or. model==8) then ! model 8 is the Hu-Sawicki model (model 7 is reserved for DGP)

  Just after model 6 block [if (model==6) then ...... end if], add the following piece of code for Hu-Sawicki,

if (model==8) then ! Hu-Sawicki

omm=CP%omegab+CP%omegac
Rbar=3.d0*((CP%H0/c*1.d3)**2)*(omm*(a**(-3))+4.d0*(1.d0-omm))
R0bar=3.d0*((CP%H0/c*1.d3)**2)*(omm+4.d0*(1.d0-omm))
mubar=1.d0/3.d0/(1.d0+fRn)
mubar=mubar*Rbar/fR0
mubar=mubar*((Rbar/R0bar)**(1.d0+fRn))
mubar=sqrt(mubar)

MG_mu = 4.d0/3.d0-1.d0/3.d0*(mubar*a)**2/(k2+(mubar*a)**2)
mudot=1.d0/3.d0/(1.d0+fRn)/fR0/((R0bar)**(1.d0+fRn))
mudot=sqrt(mudot)
mudot=mudot*(1.d0+fRn/2.d0)*(Rbar**(fRn/2.d0))*Rdot

MG_mudot = -2.d0/3.d0*k2/((k2+(mubar*a)**2)**2)

MG_gamma = 1.d0-2.d0*k2/(3.d0*(mubar*a)**2+4.d0*k2)

end if

[B2] Do the same changes in subroutine output in equations_ppf.f90

[C]  Changes in inidriver.F90

Add the Hu-Sawicki model just after the model 6 block [if (model==6) then ...... end if]

else if (model ==8) then ! Hu-Sawicki

write(*,*) '|fR0|=',fR0, 'fRn=',fRn

[D]  Add the following stuff to params.ini

#model= 8 : Hu-Sawicki

and the related parameters

|fR0|=1E-4

fRn=1.0

[E]  Recompile the code and run the code using the following settings (the Planck parameters) in params.ini

model = 8
GRtrans= 0.1
use_physical   = T

ombh2          = 0.022161

omch2          = 0.11889

omnuh2         = 0

omk            = 0

hubble         = 67.77

initial_power_num         = 1

pivot_scalar              = 0.05

pivot_tensor              = 0.05

scalar_amp(1)             = 2.21381e-9

scalar_spectral_index(1)  = 0.9611

scalar_nrun(1)            = 0

transfer_high_precision = T

transfer_kmax           = 20

transfer_num_redshifts  = 1

transfer_interp_matterpower = T

transfer_redshift(1)    = 0

If everything works, you should be able to reproduce my figure at the top of this page. Note that it is slightly different from Fig 4 in the Hu-Sawicki paper
because different cosmological parameters used.

Add the Hu-Sawicki model to MGCAMB 