## Calculating nuclear statistical equilibrium

The two previous posts explored nuclear binding energies and nuclear partition functions. We can now explore nuclear statistical equilibrium (NSE), which broadly defined, is a free energy minimum with some number of constraints. A system held at constant temperature and volume will evolve until it reaches its free energy minimum or equilibrium. When it does, it evolves no further unless some constraint on the system (such as its volume) changes.

For a system of nuclei at fixed temperature and density (that is, volume), nuclear reactions among the various nuclear species cause changes in the nuclear abundances. If the system has enough time, it reaches its free energy minimum, which is an NSE. In fact, the system typically proceeds through a hierarchy of statistical equilibria by adding or removing constraints (see, for example, the conclusion of our paper on quasi-equilibrium nucleosynthesis.
The important computational point is that, once the system attains an equilibrium, it is possible to compute the abundances without knowledge of reaction rates. We simply need to supply binding energies and partition functions and the values for other constraints (such as the ratio of the number of total neutrons to protons) in the system.

Our library of codes libnuceq permits a user to compute nuclear statistical equilibria. The easiest way to use the library is within NucNet Tools. Here are the steps I carried out. I updated NucNet Tools and compiled (if you haven't yet downloaded NucNet Tools, see the previous
post)

cd nucnet-tools-code

svn up

cd build

make -f Makefile.libnuceq all_libnuceq

When I did this, I got the following message:

../vendor/libstatmech/0.7/src/Libstatmech.c: No such file or directory

This is because there was an update to libstatmech since I last downloaded NucNet Tools codes. (You may not get such an error depending on when you last updated or downloaded the codes, but if you do, follow the next two steps.) To correct this, I ran

make -f Makefile.libnuceq cleanall_libnuceq

make -f Makefile.libnuceq all_libnuceq

This correctly compiled all the libnuceq examples. I then copied the example data over to the data_pub directory:

make -f Makefile.libnuceq libnuceq_data

I could now run an NSE calculation. To calculate the nuclear statistical equilibrium at a temperature of T9 = T / 109 K = 4, a density of 107 g/cc, and an electron fraction Ye (think of this as the total number of protons per nucleon in the system) of 0.5 (thus, equal numbers of neutrons and protons), I did the following:

cd ../libnuceq

./compute_nse ../data_pub/example_nuc.xml 4. 1.e7 0.5

The following results printed to the screen:

```Name: h1   Mass fraction: 2.211535e-03
Name: he4   Mass fraction: 7.972111e-04
Name: c12   Mass fraction: 1.807406e-10
Name: o16   Mass fraction: 1.234431e-09
Name: mg24   Mass fraction: 1.843703e-08
Name: si27   Mass fraction: 6.542801e-10
Name: si28   Mass fraction: 3.048839e-04
Name: si29   Mass fraction: 1.973449e-08
Name: p29   Mass fraction: 2.731223e-07
Name: si30   Mass fraction: 1.440652e-10
Name: p30   Mass fraction: 5.696343e-08
Name: s30   Mass fraction: 6.899328e-09
Name: p31   Mass fraction: 7.087226e-08
Name: s31   Mass fraction: 2.663897e-07
Name: s32   Mass fraction: 7.459238e-04
Name: s33   Mass fraction: 1.574008e-07
Name: cl33   Mass fraction: 3.450929e-07
Name: s34   Mass fraction: 5.996569e-09
Name: cl34   Mass fraction: 1.180223e-07
Name: ar34   Mass fraction: 9.254037e-09
Name: cl35   Mass fraction: 3.851327e-07
Name: ar35   Mass fraction: 2.779366e-07
Name: ar36   Mass fraction: 7.411405e-04
Name: ar37   Mass fraction: 2.305543e-07
Name: k37   Mass fraction: 9.956854e-08
Name: ar38   Mass fraction: 3.006487e-08
Name: k38   Mass fraction: 2.118344e-07
Name: ca38   Mass fraction: 1.896060e-09
Name: k39   Mass fraction: 1.932706e-06
Name: ca39   Mass fraction: 2.682514e-07
Name: ca40   Mass fraction: 2.216010e-03
Name: ca41   Mass fraction: 3.954164e-07
Name: sc41   Mass fraction: 6.120064e-08
Name: ca42   Mass fraction: 9.927769e-09
Name: sc42   Mass fraction: 7.686254e-09
Name: sc43   Mass fraction: 2.998909e-08
Name: ti43   Mass fraction: 1.031361e-09
Name: ti44   Mass fraction: 3.153511e-05
Name: ti45   Mass fraction: 3.313821e-07
Name: v45   Mass fraction: 7.648868e-09
Name: ti46   Mass fraction: 6.207394e-07
Name: v46   Mass fraction: 4.559378e-08
Name: ti47   Mass fraction: 5.094593e-10
Name: v47   Mass fraction: 3.685273e-06
Name: cr47   Mass fraction: 1.527548e-08
Name: v48   Mass fraction: 6.088644e-08
Name: cr48   Mass fraction: 7.988084e-04
Name: v49   Mass fraction: 1.022015e-08
Name: cr49   Mass fraction: 7.038331e-05
Name: mn49   Mass fraction: 3.186681e-07
Name: cr50   Mass fraction: 1.629572e-04
Name: mn50   Mass fraction: 5.631778e-06
Name: fe50   Mass fraction: 8.450341e-10
Name: cr51   Mass fraction: 2.941954e-07
Name: mn51   Mass fraction: 7.082886e-04
Name: fe51   Mass fraction: 1.262526e-06
Name: cr52   Mass fraction: 3.187597e-08
Name: mn52   Mass fraction: 1.717499e-05
Name: fe52   Mass fraction: 2.753519e-02
Name: mn53   Mass fraction: 8.102160e-06
Name: fe53   Mass fraction: 3.219663e-03
Name: co53   Mass fraction: 2.325983e-06
Name: mn54   Mass fraction: 2.357572e-09
Name: fe54   Mass fraction: 1.768221e-02
Name: co54   Mass fraction: 1.555930e-04
Name: ni54   Mass fraction: 3.066136e-09
Name: fe55   Mass fraction: 2.871695e-05
Name: co55   Mass fraction: 4.502576e-02
Name: ni55   Mass fraction: 1.215131e-05
Name: fe56   Mass fraction: 6.301001e-07
Name: co56   Mass fraction: 3.221973e-04
Name: ni56   Mass fraction: 8.680341e-01
Name: co57   Mass fraction: 2.118667e-05
Name: ni57   Mass fraction: 2.133487e-02
Name: cu57   Mass fraction: 3.882116e-06
Name: co58   Mass fraction: 3.850608e-09
Name: ni58   Mass fraction: 7.646524e-03
Name: cu58   Mass fraction: 1.615764e-05
Name: zn58   Mass fraction: 1.030284e-10
Name: ni59   Mass fraction: 7.415517e-06
Name: cu59   Mass fraction: 1.064034e-04
Name: zn59   Mass fraction: 9.458987e-09
Name: ni60   Mass fraction: 1.544011e-07
Name: cu60   Mass fraction: 1.168337e-06
Name: zn60   Mass fraction: 1.082488e-05
Name: cu61   Mass fraction: 1.192561e-07
Name: zn61   Mass fraction: 5.294452e-07
Name: zn62   Mass fraction: 6.331839e-07
Name: ga62   Mass fraction: 4.970190e-10
Name: zn63   Mass fraction: 1.071218e-09
Name: ga63   Mass fraction: 1.401122e-09

Global results for a T9 = 4.000000, rho = 1.000000e+07 g/cc, Ye = 0.500000 NSE:
1 - Xsum = -9.814372e-14
Number of heavy nuclei per nucleon = 1.790243e-02

Neutron chemical potential / kT = -35.039307
Proton chemical potential / kT = -15.792865
```

It is evident that most of the mass (87%) in such an NSE resides in 56Ni This makes sense since 56Ni is the species with equal numbers of neutrons and protons that has the largest binding energy per nucleon. Other diagnostics are 1 - Xsum, Number of heavy nuclei per nucleon, and the neutron and proton chemical potentials divided by Boltzmann's constant times the temperature T. The sum of all mass fractions should be unity, so 1 - Xsum should always be small (of the order 10-13 or less). The number of heavy nuclei per nucleon is the sum of the abundances of species with Z > 5 divided by the total abundance of nucleons. The neutron and proton chemical potentials are useful for analyzing an equilibrium because the fundamental NSE condition

μ(Z,A) = Z μp + (A - Z) μn

relates the chemical potential of species with atomic number Z and mass number A in an NSE to the chemical potentials of the proton μp and neutron μn.

If you have followed these steps, try NSE at other conditions (temperature, density, or Ye). You can also compute related equilibria with other libnuceq example codes. See the libnuceq tutorials.

Posted by 2012-09-02
• Anonymous - 2013-09-19

We are running your code and finding the results to be very instructive. For certain temperature, density and Ye values the output will not contain some elements ( for example we see carbon but boron and nitrogen will be missing) All three are present for slightly different conditions. It appears to me that there is a lower limit on the calculated mass fractions which governs this. Is that correct?

• Bradley S. Meyer - 2013-09-19

There is no limit on the computed mass fractions. The limit is simply on the print out.

To change that, look for the lines in nucnet-tools-code directory vendor/libnuceq/x.x/examples/compute_nse.c, where x.x is your installed version of libnuceq, that read

```/*#########################################################################
// defines
//#######################################################################*/

#define D_MIN_X_PRINT  1.e-10   /*  Minimum mass fraction to print out. */
```

This tells the code to print out only those species with mass fraction greater than 1.e-10, though all mass fractions are computed.

Change the 1.e-10 to a smaller value to print out species with mass fractions smaller than the default 1.e-10. For example, to print out all values, change to

```/*#########################################################################
// defines
//#######################################################################*/

#define D_MIN_X_PRINT  0.      /*  Minimum mass fraction to print out. */
```

Recompile the code by typing (in the build directory)

make -f Makefile.libnuceq compute_nse

and rerun (in the libnuceq directory). You can do the same for the other example codes.

Remember the example codes are to illustrate the libnuceq API. You should always feel free to modify them according to your purposes.

Last edit: Bradley S. Meyer 2013-09-20
• Comment has been marked as spam.
Undo

You can see all pending comments posted by this user  here

Anonymous - 2016-09-23

I am trying to reproduce the standard solar model values using this post
i.e we should run " ./compute_nse ../data_pub/ZZZ.xml 0.015 1.5e2 0.5"
This would mean that we have to calculate the partition function using a previous post Studying the nuclear partition function and create a new input file similar to ../data_pub/example_nuc.xml.
Firstly, does this make sens and if yes how can I do that (create my ZZZ.xml file)?

• Bradley S. Meyer - 2016-09-24

Thanks for your post. Of course the abundances in the center of the Sun are not in nuclear statistical equilibrium (NSE)--the nuclear reaction timescales at T9 = 0.015, ρ = 150 g/cc, and Ye = 0.5 are so long that achieving NSE would require a time longer than the Sun will live. Moreover, if the abundances were in NSE, the center of the Sun would be at an energy minimum and any further reactions would require net energy input to proceed.

Neverthess, you can compute the abundances that would result if the abundance distribution in the center of the Sun were miraculously to convert to NSE. To do that, you can use the JINA nuclear data file. To do this, type

cd nucnet-tools-code/examples/network

make data

This will download the data, and in particular, the file data_pub/my_net.xml, which has the masses, spins, and partition function information for the species of interest. Now you can run the NSE code. Type

cd ../../build

make -f Makefile.libnuceq all_libnuceq

This creates the libnuceq example codes. Now run the calculation by typing

cd ../libnuceq

./compute_nse ../data_pub/my_net.xml 0.015 150 0.5

You will find that the abundances are completely dominated by 56Ni, which is the species with the highest binding energy per nucleon for Ye = 0.5. Of course 56Ni is radioactive, so that species would immediately start decaying (which will change Ye) and the system would evolve until it reached 56Fe--for more information, see the libnuceq technical report low-temperature nuclear statistical equilibrium.

If you need to create your own nuclear data, you can use the libnucnet example code create_nuc_xml_from_text to convert appropriately formatted text data to XML. This post might also be of interest.

I hope all this helps. Best wishes.

Anonymous