Energy loss function discrepancies

[astronomike]

The energy loss functions in electron.c and cn_electron.py are different - the green's method assumes the input array comes in as gammas, whereas I had the crank version coming in as energies. The crank version also doesn't choose between the ISRF being switched on/off.

This is the difference between the two methods as it stands (lossOnly set to True):

After a quick change of the above things (converting crank energy loss function to gamma input array, and correctly accounting for ISRF) in cn_electron.py:

    def set_b(self,B,ne,E,ISRF=False):
        #set and return energy loss function [GeV s^-1]
        b = self.constants
        me = 0.511e-3       #[me] = GeV/c^2 
        if ISRF:
            eloss = b['ICISRF']*(E*me)**2 + b['sync']*(E*me)**2*B**2 + b['coul']*ne*(1+np.log(E/ne)/75.0) + b['brem']*ne*E*me
        else:
            eloss = b['ICCMB']*(E*me)**2 + b['sync']*(E*me)**2*B**2 + b['coul']*ne*(1+np.log(E/ne)/75.0) + b['brem']*ne*E*me

        self.b = eloss/me

        return eloss/me

This is the output:

This looks a bit better, but I'm unhappy with a few things. Firstly, just printing out the energy loss vectors calculated with each method shows there's still a fair amount of discrepancy - it's a little difficult to pinpoint this since they have different shapes, and are based on average quantities for greens but not for crank. Do you think this is expected?

Secondly, and the part I am most confused about, is the values for the SB at different frequencies in the crank method just don't behave at all like the greens ones. Here is an animation of this:

The spread of values over the WIMP masses is also a little different. It seems like there is some unexpected dependence on WIMP mass with this implementation that I'm not sure about.

Could you please take a look at this and see if you can help me spot the issue?

Another option that I was thinking of would be to just import the energy loss functions defined in electron.py for use with the crank method as well. Do you think this is preferable? If we can deal with the reshaping of the arguments, we could bypass the differing implementations of these functions.

[astronomike]

I've now noticed that the E_set is getting a factor of mxEff before being passed to the cn solver, so that probably explains the weird WIMP mass behaviour, it was being applied twice with the new change.

[Hyperthetical]

This is the correct loss function, when our E's are unitful

     def set_b(self,B,ne,E,ISRF=False):
         #set and return energy loss function [GeV s^-1]
         b = self.constants
         if ISRF:
             eloss = b['ICISRF']*(E)**2 + b['sync']*(E)**2*B**2 + b['coul']*ne*(1+np.log(E/ne)/75.0) + b['brem']*ne*E
         else:
             eloss = b['ICCMB']*(E)**2 + b['sync']*(E)**2*B**2 + b['coul']*ne*(1+np.log(E/ne)/75.0) + b['brem']*ne*E

         self.b = eloss

         return eloss

I've now noticed that the E_set is getting a factor of mxEff before being passed to the cn solver, so that probably explains the weird WIMP mass behaviour, it was being applied twice with the new change.

The mxEff factor makes it unitful, as the linspace is unitless and relative to mxEff

[astronomike]

Wouldn't we still need a factor of me in the coulomb scattering term, since I think it uses a unitless gamma term?

With your above change I get the following:

It does look like this gets rid of the big discrepancy we were seeing earlier. I'm just not sure about the above point.

[Hyperthetical]

Yes, I shouldn't have removed the me from the Coulomb term. Did you include the eloss/me or just eloss?

[astronomike]

I at first included eloss/me to make that animation. So we can disregard that.

Sorry for the confusion. It seems like the energy loss function was correct (other than the ISRF factor). I'm busy doing some debugging now, but the issue persists and I'm not sure where it's coming from anymore.

[Hyperthetical]

I also get that using eloss/me seems to fix the discrepancy. Is there any reason why this should be the case?

[astronomike]

Not that I can see. The only times that function is used are the following

• to define the energy loss timescale (only used for convergence testing - not relevant to this issue)
• in the alpha coefficients.

I assumed the coefficients would take a unitful energy loss function when I wrote this. The log transformed grid makes this a bit tricky, but I think they should still be unitful. I will double check this.