Add the colour correction

[jhmatthews]

Scattering/radiative transfer in the disc atmosphere can change the radiation temperature of the emitted radiation from an accretion disc - in particular it often makes the shape of the spectrum hotter even if the effective temperature is the same and set by local energy dissipation. This effect is often parameterised via a colour correction $f$ such that

$$B\nu (\nu, T) \to \frac{1}{f^4} B\nu (\nu, fT)$$

[ooh, github has math mode now]

where $f$ is in general a function of at least temperature if not also surface density and other things.

We should implement this colour correction factor as an option in python, which can be done fairly straightforwardly. There are various options. Done et al. 2012 give a simple form of the colour correction in which $f=1$ for $T<3\times10^4~{\rm K}$, and for $T>3\times10^4~{\rm K}$

$f=\left(\frac{T}{3\times10^4~{\rm K}} \right)^{0.82}$

I suggest we implement this as a starting point.

Proposed implementation:

• either a new disc temperature mode called "standard_fcol" or a separate question called "colour_correction". Alternatively, this could be part of "disk radiation type" which could also have bb_fcol. These questions would turn on a new mode stored in geo, presumably.
• new function that returns the colour correction when queried with a temperature. This function could check the mode variable.
• modify the relevant fluxes in the disc calculation, and also the temperature. Probably safest to only apply this for BB at the moment.

[lazygun37]

I actually have an implementation (in python -- the language) of the colour correction factor described by https://arxiv.org/abs/1809.05134

It's a little more complicated, for various reasons, but probably would still be kind of easy to port to Python (the code).

[jhmatthews]

That might be handy. I think I'll start by adding the done colour correction then we can always add more (I'll write it in such a way that we can add different colour correction prescriptions easily.

I realised two things.

• it makes more sense to define this as an additional spectral type rather than in the temperature profile question. You should still have the freedom to apply a colour correction to any temperature profile
• I think I'm write that each annulus will still emit exactly the same luminosity, so I think this actually only matters when choosing the photon frequencies. This is neat because it means that all the disk initialisation and luminosity calculations are identical, it's only the photon generation that needs a factor of "f" in the temperature when it calls the routine that generates a BB distribution of photons. Almost like we planned it all along.

[lazygun37]

OK, here are the key functions. Can obviously provide them as a file also.

def get_rel_corr(a, r): """ Function to calculate the relativistic correction factor from that we need for the surface density (the thing called m_0,rel in Davis & El-Abd, Eq 5).

These correction factors are defined in Novikov & Thorne
Their Eqs 5.4.1 defined the basic quantities we need A-E

And then the surface density expression in their Eq 5.9.10 tells us
how to combine these.

a = dimensionless spin parameter (dimensionless angular momentum)
r = radius in units of Rg
"""

A = 1.0 + a**2./r**2. + 2.0*a**2./r**3.
B = 1.0 + a/r**(3./2.)
C = 1.0 - 3.0/r + 2.0*a/r**(3./2.)
D = 1.0 - 2.0/r + a**2.0/r**2
E = 1.0 + 4.0*a**2.0/r**2.0 - 4.0*a**2.0/r**3.0 + 3.0*a**4.0/r**4.0

rel_corr = A**(-2.0) * B**(3.0) * C**(1./2.) * E * D**(-1.0)

return rel_corr

def get_surf_dens(M_BH, Mdot_acc, a_spin, R, R_g, alpha): """ calculates the surface density in the disk as a function of radius """ r_scaled = R / R_g relcorr = True if relcorr: rel_corr = get_rel_corr(a_spin, r_scaled) else:
rel_corr = 1.0

#print()
#print('rel_corr =', rel_corr)
#I'm going to do this via Eq 5.9.10 from NT, rather than
#from Davis & El-Abd, because I don't see how/why the radiative
#efficiency in their case comes into it.

surf_dens =   20.0 * alpha**(-1) * (M_BH/(3.0*M_sun)) * (Mdot_acc/1.e17)**(-1.0) \
            *  (R/R_g)**(3./2.) * rel_corr

return surf_dens

def get_f(R_g, R, T_eff, a_spin, surf_dens): """ source: https://arxiv.org/pdf/1809.05134.pdf

Function which calculates the spectral hardening factor 

Inputs:
    R_g: gravitational radius 
    R: radius array (of the disk)
    T_eff: effective temperature 

Outputs:
    f: spectral hardening factor

"""

Q = (c**2.0)/((R_g**2.0) * ((R/R_g)**3.0)) 
lQ = np.log10(Q)

#print(Q)
#print(lQ.shape)
lT = np.log10(T_eff)
lsurf_dens = np.log10(surf_dens)

#print('lsurf_dens = ', lsurf_dens)

f =   1.74 + \
    + 0.08*(a_spin/0.9) \
    + 1.06*(lT - 7.0) \
    + 0.14*(lQ - 7.0) \
    - 0.07*(lsurf_dens - 5.0)

#print('a_spin shape ', a_spin.shape)
#print('lQ shape ', lQ.shape)
#print('surf_dens shape ', surf_dens.shape)
#print('lT shape ', lT.shape)

f = np.where(f > 1.0, f, 1.0)
return f

[jhmatthews]

Thanks! I realised that my thinking above about only needing to modify the frequency distribution is wrong, because of banding. The colour correction doesn't change the total luminosity from each annulus, but it does change the band limited luminosity, so as a result we do need to change the routines that calculate the disc luminosity and split it up into rings (especially since the rings are, I think, defined separately for each band).

[jhmatthews]

I've implemented an initial version of the Done colour correction. The results look sensible (below). Here's a test of an initial implementation of the colour correction compared to disc spectra from a python script I have. It looks like the colour correction is working.

I'm not sure if the disagreement in the Wien tail means my python script is inaccurate, or Python the MCRT code is inaccurate - happens without a colour correction so it's not new, and it's entirely possible Python is working perfectly.

This could either be merged now or I could implement some other forms of the colour correction.

[jhmatthews]

This is now included in the code albeit only with the Done model.