## Running a network calculation with simple hydrodynamics

In my recent posts, I have been running calculations with the network/run_entropy.cpp examples code. These calculations have had the parameter τ set to infinity, so the density has remained fixed. In this post, I will show how to use the code to run network calculations with simple one-zone hydrodynamics.

The starting point is to consider a ball of radius R. We consider the ball to have uniform properties, such as the density, which can change as the ball's radius expands or contracts with time. The run_entropy.cpp code solves a set of three similtaneous first order differential equations for the three quantities x[0], x[1], and x[2]. Together these quantities form a state vector x. The three components of x are, first, the scaled radius

where R0 is the initial radius of the ball. The second component is the scaled velocity

The third component is the entropy per nucleon s:

To solve for the time evolution, we must provide the time derivatives for these quantities, which are encapsulated in a state vector dxdt. The time derivative of x[0] is the velocity:

The time derivative of x[1] is the scaled acceleration:

The time derivative of x[2] is the entropy generation (or loss) rate:

Since the ball is of uniform density, we can find the density in terms of x[0]:

These are the general equations. To solve a particular problem, we must supply an acceleration equation and the initial values of the components of x. The default problem in NucNet Tools is for an exponentially expanding density. In this case, we find the time dependence of x[0] as

From this we can supply the acceleration equation:

In previous calculations, such as these, I set x[0] = 1, x[1] = 0, x[2] to the initial entropy. I also set τ = ∞. With this value of τ, the acceleration is zero, and the ball remains at fixed radius (and density). In this calculation, I will set τ = 0.1 seconds. This will cause the radius of the ball to grow exponentially on a time scale of 0.3 seconds.

I will choose to run calculations that include screening and NSE corrections to the entropy. I set up the code appropriately in the my_examples/network directory, and I edit the ../../data_pub/zone.xml file to read

```<zone_data>
<zone>
<optional_properties>
<property name="tend">2.</property>
<property name="tau">0.1</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">yes</property>
<property name="use nse correction">yes</property>
</optional_properties>
<mass_fractions>
<nuclide name="he4">
<z>2</z>
<a>4</a>
<x>1</x>
</nuclide>
</mass_fractions>
</zone>
</zone_data>
```

Note that τ is set to 0.1. Also note how x[0] is set to 1 and x[1] is set to 0 through the tags. The code computes the initial x[2] from the input conditions. I run a calculation with and without entropy generation by typing

./run_entropy ../../data_pub/my_net.xml ../../data_pub/zone.xml my_output_sdot.xml --nuc_xpath "[z <= 50]" --reac_xpath "[not(product[contains(.,'neutrino')])]"

and

./run_entropy ../../data_pub/my_net.xml ../../data_pub/zone.xml my_output_no_sdot.xml --nuc_xpath "[z <= 50]" --reac_xpath "[not(product[contains(.,'neutrino')])]" --sdot_nuc_xpath none

Note that I am still not including weak interactions in these calculations.

Once the calculations are done, I use the analysis/print_properties.cpp code to obtain the time, temperature, and entropy per nucleon, as before. This plot shows the graph of the total entropy per nucleon as a function of time in the calculation:

Of course the calculation without entropy generation has constant entropy throughout the calculation. In the calculation with entropy generation, almost all the entropy generation occurs in the first nanosecond. The temperature as a function of time looks similar to our previous calculations except that after ~0.01 seconds, the expansion of the ball begins and the density and temperature fall:

If I plot the final elemental abundances in the calculations, I get this figure:

The abundances are dominated by Ni and He. The abundances in the calculation with entropy generation show more light species, which also gives a higher abundance of some of the heavier species, than the calculation without entropy generation. This is in accordance with the former's higher final entropy.

If I compare the elemental abundances to equilibrium, I can create entropy generation and no entropy generation movies. These movies show that the elemental abundances achieve nuclear statistical equilibrium (NSE) but then fall out of NSE as the material expands and cools. The abundances pass through quasi-equilibrium (QSE) before freezing out.

The above considerations show a network calculation with simple hydrodynamics. A user can make a different calculation by supplying an appropriate acceleration function (the default is in user/hydro_helper.cpp). To do this, the user would create his or her own hydro_helper.cpp/h files (preferably in my_user) with the proper acceleration. The user would then set that acceleration function in his or her own version of run_entropy.cpp (preferably in my_examples/network).

I cut off the calculations above at a final time of 2 seconds. The reason is that, for longer times, the density and temperature got into a region in which the default screening function went beyond the strong screening regime. To treat this, I would simply provide a screening function that appropriately handles this low temperature regime. I leave that as an exercise for the interested user.

Posted by 2016-08-14

Anonymous