Plasmapy: Add functionality to find fusion reaction rates, cross sections, and energy generation rates

[namurphy]

In laboratory fusion devices and in stellar astrophysics, it's important to be able to calculate quantities associated with different fusion reactions. We will need a way to calculate the following quantities as functions of temperatures and densities:

• Reaction rates (the number of reactions per cubic meter per second)
• Reaction cross sections (in square meters)
• Energy generation (in joules released per cubic meter per second)

We have some of the capability already in the atomic package, including nuclear_binding_energy, energy_from_nuclear_reaction (which might be worth renaming nuclear_reaction_energy), and half_life. In particular, energy_from_nuclear_reaction has code that parses reactions. Some of the reactions will be radioactive decay, so we will need the capability to calculate rates for that too.
It will probably be reasonable for the first pull request to be able to do only 'D + T --> He-4 + n' since getting that will be tricky enough!

There is a huge number of reactions that happen, so it will probably be worth working with applications in this order:

1. Laboratory fusion devices
2. Stellar interiors (especially the CNO cycle and the proton-proton chain)
3. Supernovae nucleosynthesis

Reactions that are important in the laboratory include:

• 'D + T --> He-4 + n'
• 'D + D --> T + p'
• 'D + D --> T + p'
• 'D + He-3 --> He-4 + p'
• 'T + T --> He-4 + 2n'
• 'He-3 + He-3 --> He-4 + 2p'
• 'He-3 + T --> He-4 + p + n'
• 'He-3 + T --> He-3 + D'
• 'D + Li-6 --> 2He-4'
• 'D + Li-6 --> He-3 + He-4 + n'
• 'D + Li-6 --> Li-7 + p'
• 'D + Li-6 --> Be-7 + n'
• 'He-3 + Li-6 --> He-4 + He-3'
• 'p + B-11 --> 3He-4'

I got these from the Wikipedia article on nuclear fusion, it turns out. The articles on the CNO cycle and proton-proton chain should also be helpful. There's even more in the article on stellar nucleosynthesis.

We'll probably want to first assume cases with a Maxwellian particle distribution and have the corresponding input be a temperature. Things will get more complicated if we have a neutron/particle beam so that we can't assume a particular temperature...but we can probably ignore that for now.

There's a chance that we could learn from or use other packages (like PyNE, though I think that may be more on the nuclear engineering side of things which is kind of different).

This will need documentation and tests as well. Reactions that aren't implemented should raise a NotImplementedError.

My sincere thanks!
-Nick

[namurphy]

A review article: Solar fusion cross sections

[jessicalostinspace]

Is there any reason why the input is a string that is being parsed through here? Further development would be much more convenient if the function took in two lists: reactants and products (e.g. nuclear_reaction_energy("D + T --> alpha + n") would become nuclear_reaction_energy(reactants=[D, T], products=[alpha, n]). That way, if you were checking against a database of reactants to products to automate this process, you could easy check by lengths of the arrays => possible reactions.

[jessicalostinspace]

Also, I'm not an expert on this BY ANY MEANS but as far as the initial temperature you mentioned in the description - is that the ignition temperature? If so, shouldn't it default to that temperature depending on the reaction? I haven't taken a physics class in years but was reading more into this to tackle this issue: https://www.ems.psu.edu/~radovic/Chapter14.pdf

[namurphy]

Is there any reason why the input is a string that is being parsed through here? Further development would be much more convenient if the function took in two lists: reactants and products (e.g. nuclear_reaction_energy("D + T --> alpha + n") would become nuclear_reaction_energy(reactants=[D, T], products=[alpha, n]).

@jessicalostinspace - This is a great idea! You're right that we should have this be able to take a list of reactants and a list of products. I put this as a string to maximize readability, but I was operating under the (incorrect!) assumption that people would only care about two or three reactions at a time. This does ignore nucleosynthesis in massive stars and supernovae which have to deal with a lot of reactions. Thank you for the code review!

I wonder if the best way to do this would be to allow both the string reaction input (which would be most readable to those of us who care only about a couple of reactions) as well as your idea for the reactants and products keywords (which will be most useful for anyone who has to iterate). Please feel free to make this change in nuclear.py if you would like!

[namurphy]

Oof, it took me a few minutes to figure out how to add an image in one of these comments!

The temperature is not the ignition temperature but rather the temperature of the plasma at the reaction is occuring in, so the goal is to figure out how much each reaction is happening at different temperatures as in this plot.

[Image description: Fusion reactivity as functions of temperature for deuterium-deuterium, deuterium-tritium, and deuterium-helium-3 reactions. The reactivity for each reaction increases with temperature until it reaches a maximum and then gradually declines. Licensed under CC BY 2.5 under the associated link.]

[jessicalostinspace]

@namurphy would you suggest opening up a ticket for the refactor? or just branching off on this ticket? It should probably be refactored before implementing the logic for reaction rates, etc. They should also be in two different PR's. There won't really be a need to parse a string so a lot of that code can be cleaned up and more succinct. Thank you for your feedback :)

[namurphy]

Branching off this ticket sounds fine to me, and I agree that two separate PRs are needed.

I'm wondering if we should require that reactants and products be keyword-only arguments. That would require that resulting code be more clear and explicit, but perhaps it's not necessary since maybe it would be common sense that reactants would go first in a call like nuclear_reaction_energy(['Be-8'], ['He-4', 'He-4']). Either way sounds good to me.

Again, thank you!

[namurphy]

Another resource is SkyNet which has nucleosynthesis data for astrophysics (see Lippuner & Roberts 2017). This uses a BSD 3-clause license as well.