In a previous post, I ran a network calculation at constant entropy. In a subsequent post, I computed the entropy generation rate for that network calculation. The modeled system generated entropy during its expansion; thus, holding the entropy constant is not a consistent treatment. In this calculation, I will include entropy generation in my network calculation.
To do so, I use the examples/network code run_entropy.cpp in NucNet Tools. This code evolves, in addition to the network abundances, a vector x with three components. The first component x is a scale factor. The second component x is the time derivative of the first component: x = dx/dt. The third component x is the total entropy per nucleon. The user supplies expressions for dx/dt and for the density and temperature in terms of x and dxdt. dx/dt = x, as noted before. dx/dt = d2x/dt2 , while dx/dt is the entropy generation rate.
The default functions are defined in user/hydro_helper.c/h. Here, the density is given by ρ(t) = ρ0 x-3 = ρ0 e-t/τ, where τ is the e-folding time scale for the density. This means that x = et/3τ. The temperature is found from inverting the entropy. dx/dt = x / 3τ. To use other expansions, one simply defines x, dxdt, ρ(t), and the temperature appropriately.
In this post, I will consider a system that is not expanding; thus, τ = ∞. I let the system begin at a temperature T9 = T / 109 K = 2 with the initial abundances all 4He. I will not include weak interaction rates, and I will use a network that includes species up to Z = 60 (neodymium). The system should evolve into nuclear statistical equilibrium (NSE).
I begin by updating NucNet Tools and compiling. I type
I then edit ../../data_pub/zone.xml until it reads
<zone_data> <zone> <optional_properties> <property name="tend">1.</property> <property name="tau">inf</property> <property name="munuekT">-inf</property> <property name="t9_0">2.</property> <property name="steps">5</property> <property name="rho_0">3.e9</property> <property name="x" tag1="0">1</property> <property name="x" tag1="1">0</property> <property name="use screening">no</property> <property name="use nse correction">no</property> </optional_properties> <mass_fractions> <nuclide name="he4"> <z>2</z> <a>4</a> <x>1</x> </nuclide> </mass_fractions> </zone> </zone_data>
Notice for this calculation I do not use screening or the reverse ratio correction. Notice also that I set the initial x = 1 and the initial x = 0. Because τ = ∞, this will keep the density at the fixed value of 3 x 109 g/cc throughout the calculation. I execute the code by typing
./run_entropy ../../data_pub/my_net.xml ../../data_pub/zone.xml my_output.xml --nuc_xpath "[z <= 60]" --reac_xpath "[not(product[contains(.,'neutrino')])]"
The nuclear XPath expression limits the network to Z up to 60 while the reaction XPath expression excludes weak reactions. The calculation takes a few hours to run.
Once the calculation finishes, I can plot the temperature versus time. I make sure the analysis codes are complied by typing
I then type
../analysis/print_properties my_output.xml time t9 > props.txt
If I plot column 3 of props.txt versus column 2, I get this graph
Here I see that the temperature rises dramatically at about 1 nanosecond from T9 = 2 to the final value of T9 = 10.318, at which point the system has attained NSE and no longer evolves. I can compute the entropy per nucleon during the calculation by compiling and running the relevant codes. I type
../thermo/compute_thermo_quantity my_output.xml "total entropy per nucleon" > s.txt
I plot column 2 of s.txt versus column 2 of props.txt to get this figure:
Here I see that the entropy rises from its initial value of 2.194 kB, where kB is Boltzmann's constant, to 2.830 kB. Nuclear reactions generated 0.636 kB units of entropy during the evolution. I compute the entropy generation rate by typing
../thermo/compute_sdot my_output.xml > sdot.txt
I plot column 4 of sdot.txt versus column 3 of props.txt to get this figure:
The entropy generation rate is high throughout the evolution until equilibrium is attained at which point it plummets.
Finally, I graph the elemental abundances. I type
../analysis/compare_equil my_output.xml > equil.txt
From these data, I make the elemental abundances movie. Here I see that the abundances indeed evolve into equilibrium. Once they attain NSE, they would only change by weak reactions, which have been disabled in this calculation.