Brief description

I found the evolution of the normalized axion density rho{ax}/rho{ax,initial} derived from loga_table and grhoax_table has a significant deviation from the analytical solution in radiation-dominated case (m_ax\~10^-22eV). The normalized density is \~25% lower than the analytical solution at a\~a_osc (Fig 1). The deviation is expected to be tiny because the universe is well radiation-dominated at the epoch involved. This issue is independent of how the initial axion field is determined due to the normalization. A full picture of the deviation can be obtained if one evolves the background equations toward much later phase. The normalized density as a function of scale factor in log-scale plot is shifted by a horizontal displacement (Fig 2).

I have found out that the inconsistencies in the axion_background module cause the issue. The result of axionCAMB and the analytical solution match perfectly if proper modifications of the code are made (Fig 3 and Fig 4).

Details

My input parameters

I used these parameters

    use_physical = T
    ombh2 = 0.0
    omch2 = 0.07595
    omnuh2 = 0
    omk = 0
    hubble = 70

    use_axfrac=F
    omaxh2 = 0.07595
    m_ax = 1.493183773790898e-22

    temp_cmb = 2.7255
    massless_neutrinos = 3.046
    massive_neutrinos = 0

The parameters above give Omega_m=0.31, f_ax=f_cdm=0.5, m_ax/H_0=1.0e11, and Omega_r=8.538e-5.

Comparison of the result of axionCAMB and the analytical solution

I added the following statements after the calling of w_evolve in the driver inidriver_axion.F90 to print the contents of loga_table and grhoax_table.

      do i = 1, ntable
          print *, 10.0_dl ** (P%loga_table(i)), 10.0_dl ** (P%grhoax_table(i))
      end do

      stop 'axion background evolution test done'

The output gives the computed axion density as a function of scale factor.

The analytical solution of the evolution of the normalized axion density in radiation-dominated case is

rho{ax}/rho{ax,initial}=(Gamma(5/4))^2*(mt/2)^(-1/2)[(J{1/4}(mt))^2+(J{5/4}(mt))^2]

where Gamma is the gamma function and J is Bessel function of the first kind. m is the axion mass and t is cosmic time (not conformal).

Fig 1 shows the comparison of the result of axionCAMB and the analytical solution. The relative error of axionCAMB in the normalized density is about 25% at a\~a_osc provided that the analytical solution is exact.

Fig1 The deviation of the normalized axion density (in linear scale).

A notable feature appears if the background equations are evolved toward much later phase. The normalized density as a function of scale factor in log-scale plot is shifted by a horizontal displacement (Fig 2).

Fig2 The "shift feature" of the result of axionCAMB (in log scale). m=30H_osc is used just to evolve the background equations toward much later phase.

Inconsistencies in the axion background module

The shift feature shown in Fig 2 implies the time scale is "zoomed". The module axion_background uses dimensionless variables for time and other quantities. I therefore inferred that these variables may be introduced in an inconsistent way.

Here are some code lines related to axion density in the w_evolve subroutine of the axion_background module

      forall(i=1:ntable)

*     !...some code...

      grhoax_table_internal(i)=((v_vec(2,i)/a_arr(i))**2.0d0+(maxion_twiddle*v_vec(1,i))**2.0d0)
      end forall
*     !This is the axion energy density *h^2/(3H_0^2/8\pi G), where h is the dimensionless Hubble Parameter

The code above introduces a dimensionless axion density defined as rho_ax*h^2/rho_crit0 (rho_crit0 is the critical density at present). This definition is consistent with the other parts of w_evolve where axion density is involved.

The definition, however, is not compatible with the background equations in the derivs subroutine of axion_background. The definition implies the dimensionless axion field v_1 and its derivative v_2 are given by

phi=sqrt(3/(4pi G))v_1/h
phi_dot=sqrt(3/(4pi G)H_0v_2/h

In terms of v_1 and v_2, the background equations thus should be

dv_1/da=v_2/(a*H/H_0)
dv_2/da=-2v_2/a-(maxion_twiddle)^2*a*v_1/(H/H_0)

where maxion_twiddle is the dimensionless axion mass m_ax/H_0. The background equations in derivs are

      dvt_da(1)=v(2)/(a*(lhr))
      dvt_da(2)=-2.0d0*v(2)/(a)-(maxion_twiddle**2.0d0)*a*v(1)/(lhr)

The variable lhr should be H/H_0, but it is actually H/(100km/s/Mpc) in the code.

Solution to the issue

The issue can be solved by modifying the background equations in the derivs subroutine of the axion_background module, which is shown as follows.

*     !hsq is the square of H_0/(100km/s/Mpc) in the "derivs" subroutine
*     !Multiply the relevant terms of the original background equations by dsqrt(hsq)
      dvt_da(1)=v(2)/(a*(lhr))*dsqrt(hsq)
      dvt_da(2)=-2.0d0*v(2)/(a)-(maxion_twiddle**2.0d0)*a*v(1)/(lhr)*dsqrt(hsq)

Fig 3 and Fig 4 show that the result of axionCAMB and the analytical solution match perfectly after modifying the background equations.

Fig3 The evolution of the normalized axion density after modifying the background equations (in linear scale).

Fig4 The evolution of the normalized axion density after modifying the background equations (in log scale). m=30H_osc is used just to evolve the background equations toward much later phase.

Thanks astralsight for bringing this to my attention! I need to check my detailed notes to make sure I agree with all the definitions and relevant factor of little h being carried around, but it looks like you have been very careful so I will get back to you soon. Once nice thing is this lets us trust background evolution for much longer!

Thanks astralsight for bringing this to my attention! I need to check my detailed notes to make sure I agree with all the definitions and relevant factor of little h being carried around, but it looks like you have been very careful so I will get back to you soon. Once nice thing is this lets us trust background evolution for much longer!

Hi. Sorry to bother you again after several months. Recently I am trying to get more accurate matter power spectra by the axionCAMB code which is slightly modified based on my solution to the previously descriebed issues of the axion_background module. I expect the results of axionCAMB after my modifications won't change much for axions with m_{ax}>>H_0, because they begin to coherently oscillate at very high redshifts and thus the axion_background module is only needed in a short time. But I found that notable deviations between the results of matter power spectrum after the modifications and those computed by the original program can occur at some scales. The typical value of the relative deviation is \~20% at k\~(0.02-0.2) hMpc^-1 for m_{ax}/H0=10^6 and \~10%-30% at k\~(10-20) hMpc^-1 for m{ax}/H_0=10^11 (Omegam=0.31, h=0.7, Omega{ax}/Omega{DM}=0.9). I also noticed the matter power spectrum computed by the modified code has very deep trough which is several orders of magnitude lower than the original result at a specific k value in the two models respectively. So far what I can confirm is that my solution given this issue report can eliminate the inconsitencies between the theoretical and computed results of the evolution of the normalized axion density rho{ax}/rho_{ax,ini}. In view of the deviations seen in matter power spectra, is there any incompleteness of my solution? Do I need to modify the other modules of axionCAMB?

Dear astralsight - I also need to look into this, but since the potential error you found can be thought of as an improperly normalized mass, perhaps the sound speed (in the perturbation equations) is worth re-examining as well - in particular its parametric dependence on the mass. I too will look into this. Let me know what you find!

Dear astralsight - I also need to look into this, but since the potential error you found can be thought of as an improperly normalized mass, perhaps the sound speed (in the perturbation equations) is worth re-examining as well - in particular its parametric dependence on the mass. I too will look into this. Let me know what you find!

Thank you for reminding me of a more complete solution to the issues of the axion_background module. My previous solution is confined to only modifying the background ODEs and ignores the other aspects. I examined the source code in axion_background.f90 in more details, and found all equations (including the adiabatic sound speed cs2_table) are essentially consistent if we treat the dimensionless axion mass maxion_twiddle as m_ax/(100km/s/Mpc in eV). But the value of maxion_twiddle is m_ax/(H_0 in eV) according to the original intention of code design, which leads to the unexpected results. Therefore, we simply need to multiply the maxion_twiddle variable by h=H_0/100km/s/Mpc in its initialization statement to solve the issues. The new solution has been tested and it gives consistent background evolution. I also examined the source code related to the perturbation equations in equations_ppf.f90. The code defines its own varible EV%m_ax which is initilized to m_ax/(H_0 in eV). This definition is consistent with the calculation of the another sound speed csquared_ax. So it seems that we don't need further modifications in equations_ppf.f90. I found the matter power spectrum has much smaller changes as expected after applying the new solution. Maybe it is the final solution?

astralsight, would you be willing to share contact info? - I think you may be onto the solution here, but I want to check that this can't be just addressed at the input level (interpreting mass correctly when comparing with analytic solutions). We (Doddy, Renee, and I), are also interested in implementing the axionCAMB changes on a much more recent version of camb and would be interested in collaborating. Our emails are dgrin@haverford.edu, hlozek@dunlap.utoronto.ca, david.j.marsh@kcl.ac.uk

astralsight, would you be willing to share contact info? - I think you may be onto the solution here, but I want to check that this can't be just addressed at the input level (interpreting mass correctly when comparing with analytic solutions). We (Doddy, Renee, and I), are also interested in implementing the axionCAMB changes on a much more recent version of camb and would be interested in collaborating. Our emails are dgrin@haverford.edu, hlozek@dunlap.utoronto.ca, david.j.marsh@kcl.ac.uk

My email is qdzh_yclhq1@163.com, associated with my github account and also used in local git repositories. The topic of my dissertation involves matter power spectrum in ultralight axion model. I had considered developing my own code to compute the power spectrum in axion models. For simplicity I planned to firstly build an equivalent version of old CAMB (only for the Boltzmann code) using the C language (that I am more familar with ) and then inplement axions in it. Therefore it may be easy to introduce the new features of the latest CAMB when needed. Focusing only on the Boltzmann code, there seems to be few modifications compared with the old one, but the code structure changes a lot. The object oriented programming (OOP) paradigm is widely used in the new CAMB, which leads to many "global variables" (only public variable is inheritable in a Fortran class). This may cause great difficulties in code development. Due to the time consumption, I paused the relevant works but I have learned a lot about the inplementaion details of CAMB and its physics. Now I am using axionCAMB for my research project. I would be glad to contact with you in the future days.

dgrin1 / axionCAMB

Discrepancy between the computed axion background tables and the analytical solution at radiation epoch #6

Brief description

Details

My input parameters

Comparison of the result of axionCAMB and the analytical solution

Inconsistencies in the axion background module

Solution to the issue