Status of Positron capture bug?

[chrisCastillo-physics]

In line 101 in src/weakrates/interface.f90, I see a stop condition if the emissivity from positron capture is attempted:

stop "Positron capture effective q is bugged and not currently working,&
please turn it off in requested_interactions.inc."

I was wondering the current status of this bug and how I could attempt to fix it myself.

[evanoconnor]

The issue is just a coding one at this point. It can be done.

The following lines:

rbeta = return_weakrate(weakratelib,A,Z,t9,lrhoYe,idxtable,4)
rcap = return_weakrate(weakratelib,A,Z,t9,lrhoYe,idxtable,5)
rnu = return_weakrate(weakratelib,A,Z,t9,lrhoYe,idxtable,6)
qec_eff = weakratelib%tables(idxtable)%nuclear_species(weakratelib%tables(idxtable)%nucleus_index(A,Z),1)
avgenergy(1) = rnu/(rcap + rbeta)
avgenergy(2) = qec_eff

needs to be correctly implemented for positron capture. I believe there is an issue here where, for positron capture, there A and Z are incorrect and the wrong rate is returned.

Another issue is the effective spectrum used is the approximate rates is based on electron capture, not positron capture, this would have to be modified.

I would be happy to work through this with you. I have done it before, but only in a hacked way.

[chrisCastillo-physics]

Great, although I must profess my fortran knowledge is a little lacking.
For example, when looking at the definition of return_weakrate it uses a ton of "this" references as well as "%" that makes it much harder to parse through compared to the rest of Nulib.

As far as debugging/rewriting what would be the first step? Isolate this piece of code and check which rates returned match the tables?

[evanoconnor]

Yes, the syntax of the weak rates library is a bit different from the rest of NuLib, this is because it comes from a different author.

As a starting point, pick a density, temperature, and electron fraction that match a point of the electron capture tables. First try and reproduce the electron capture rate for a particular nucleus then check the positron capture rate.

This is where you should see the first issue. The following is the first few lines of the LMP electron capture table:

neg. daughter Sc46 z=21 n=25 a=46 Q= -0.8653
pos. daughter Ca46 z=20 n=26 a=46 Q= 0.8650
+++ Sc46 --> Ca46 +++ --- Ca46 --> Sc46 ---
t9 lrho uf lbeta+ leps- lrnu lbeta- leps+ lrnubar

The entries for electron capture on Sc46 (leps-) and positron capture on Ca46 (leps+) are in the same table entry, but are for different parent nuclei. I think the place to correct this in subroutine microphysical_electron_capture in interface.F90 by appropriately choosing the right A and Z for the corresponding species based on the value of neutrino_species.

[chrisCastillo-physics]

Weak rates are only implemented with the Hempel EOS. This gives A,Z and respective abundances as a function of purely thermodynamic variables. So to directly compare an entry in the table with Nulib output requires commenting out the prepopulation of abundances (the call to nuclei_distribution_Hempel in microphysical_electron_capture ) and just making microphysical_electron_capture a function of ints A and Z.

Reactions are referenced by the A,Z for the parent in electron capture. So I am pretty sure to fix this is as easy as changing all Z to Z+1 in the function emissivity_from_weak_interaction_rates (when (neutrino_species.eq.2)):

rbeta = return_weakrate(weakratelib,A,Z+1,t9,lrhoYe,idxtable,4)
rcap = return_weakrate(weakratelib,A,Z+1,t9,lrhoYe,idxtable,5)
rnu = return_weakrate(weakratelib,A,Z+1,t9,lrhoYe,idxtable,6)
qec_eff = -weakratelib%tables(idxtable)%nuclear_species(weakratelib%tables(idxtable)%nucleus_index(A,Z+1),1)

So I implemented the above, and it runs and gives an output, however I am having issues comparing the output with the tables.

[evanoconnor]

One quick way to test a single rate is to reset all the mass fractions to 0 except 1.

I think the A,Z+1 has to occur higher up the chain of subroutines. If the A,Z rate exists for electron capture, it does not mean the A,Z+1 rate for electron capture rate exists (which would contain the positron capture rate for A,Z).

[chrisCastillo-physics]

I believe it is sufficient to add to interface.f90 at line 457:
if (neutrino_species .eq. 2) then Z = Z+1 endif
As tables are referenced by A,Z for electron capture, this ensures we are in the right table, and then we grab the correct positron capture columns in emissivity_from_weak_interaction_rates.