Minor Changes to Neutrino Cooling Rate

[gchauhan42]

I have been working on a project (along with Mathieu Renzo and Cecilia Lunardini) involving the effect of Beyond-the-Standard-Model neutrino physics on stellar evolution. This requires changing the existing neutrino cooling implementation in MESA. To achieve this, I used the flag (s% other_neu) in run_star_extras.f90. My modified neutrino cooling rate is based directly on the standard MESA neutrino cooling rate implemented in mod_neu.f90 at the main GitHub page for MESA (https://github.com/MESAHub/mesa/tree/main/neu/private).

I want to report two issues that I found during this time :

1. Hardcoded Constant in the Plasma-Neutrino process : In mod_neu.f90, all cooling rates are organized in the form of a thermal integral times some function of pre-defined factors labeled as t1,t2,t3,t4,t5,t6 (see lines 38-43). These pre-defined factors are determined by a combination of cv, cvp and ca, cap (vector and axial-vector couplings, see lines 34-37). Vector coupling is a function of weinberg_theta angle (line 34). Most of the cooling processes (photo-neutrino, pair-neutrino, bremsstrahlung-neutrino, recombination-neutrino) use these pre-defined factors defined globally. Therefore, changing the weinberg_theta or any of the (cv, cvp, ca, cap), will affect all these neutrino cooling rates.

However, I found that the dependence on (cv, cvp) for the plasma-neutrino process is explicitly hard-coded : value of 0.93153d0 (lines 1582,1583). The expression for this pre-factor is in equation(4.1) of Itoh1996. If evaluated for weinberg_theta used in Itoh1996 (see eqns 2.2-2.4) and set n=2, we get exactly 0.93153. Therefore, if any global change is done for weinberg_theta and/or num_neu_fam, the change is not reflected in the plasma-neutrino process. Moreover, in const_def.f90, definition of weinberg_theta = 0.22290d0 (line 111). While plasma-neutrino process (the hard-coded constant) uses Itoh et. al.'s value of weinberg_theta = 0.2319.

Fix is to define a new factor (around line 37-38)--> real(dp), parameter :: tfac0 = cvcv + (num_neu_fam-1.0d0) * (cvpcvp)

Original line 1582 : a2 = 0.93153d0 * 3.0d21 * input% xl9
Original line 1583 : a3 = 0.93153d0 * 3.0d21 * 9.0d0input% xl8input% xldt

Fix for line 1582 : a2 = tfac0 * 3.0d21 * input% xl9
Fix for line 1582 : a3 = tfac0 * 3.0d21 * 9.0d0input% xl8input% xldt

2. Test conditions are for T/Rho instead for temp/den : This major issue pertains to variable assignment and test conditions defined in the code block (lines 230-280). The code correctly allows either T or logT, and either Rho or logRho, but later it tests if (T <= 0) and if (Rho <= 0). If MESA calls the hook with T = arg_not_provided and only logT supplied, or similarly for density, those checks will fire on the sentinel value and wrongly abort the routine. Those checks should use temp and den instead (line 240 and line 266).

Original line 240 : if (T <= 0d0) then
Fix for line 240 : if (temp <= 0d0) then

Original line 246 : logtemp = log10(T)
Fix for line 246 : logtemp = log10(temp)

The original lines are wrong because when MESA passes T = arg_not_provided and only supplies logT instead. In that case, T is a sentinel, not the actual temperature. Similarly for density,

Original line 266 : if (Rho <= 0) then
Fix for line 266 : if (den <= 0d0) then

Original line 272 : logden = log10(Rho)
Fix for line 272 : logden = log10(den)

The code after the above-mentioned changes does show some differences compared to the older results, but they seem quite minimal.

[fxt44]

I found that the dependence on (cv, cvp) for the plasma-neutrino process is explicitly hard-coded : value of 0.93153d0

yes. i own this bug as its even present in my originating source code. thanks for the clear fix!