© ESO, 2014

1. Introduction

A key nucleosynthesis mechanism for producing heavy elements beyond the iron peak is rapid neutron capture or the r-process (Burbidge et al. 1957; Cameron 1957). In spite of a large volume of observational data (Sneden et al. 2008), where the r-process occurs remains an open question. Explosive, neutron-rich environments provide the ideal conditions for r-processes to occur. The predominant astrophysical sites being studied as possible locations for r-process: the high-entropy winds (HEW) from proto-neutron stars (see Qian & Woosley 1996; Farouqi et al. 2010) and ejected matter from neutron star mergers (see Freiburghaus et al. 1999a; Goriely et al. 2011). However, both of these scenarios face significant hurdles that must be overcome. The long time-scale for neutron star merger events (Faber & Rasio 2012) limits this scenario’s ability to explain the r-process element enrichment of metal-poor stars (Sneden et al. 2003). HEW models have been shown to be very sensitive to the chosen physical conditions of the winds, and elaborate hydrodynamic models have yet to prove that the HEW scenario can provide the necessary environment for significant r-process to occur (e.g. Hoffman et al. 2008; Janka et al. 2008; Roberts et al. 2010; Fischer et al. 2010; Wanajo et al. 2011). The theoretical quark nova has also been proposed as a potential r-process site (Jaikumar et al. 2007). The explosive, neutron-rich environment of a quark nova provides an intriguing avenue of study for astrophysical r-process (Jaikumar et al. 2007). It is important to note that the observed nuclear abundance of r-process elements in metal-poor stars (Sneden et al. 2003) and abundance data of certain radionuclides found in meteorites (Qian & Wasserburg 2008) point to the likelihood of multiple r-process sites (Truran et al. 2002). Certainly there is much to be learned about the astrophysical r-process and its necessary conditions. To help drive this study we have developed r-Java (Charignon et al. 2011), which is a cross-platform, flexible r-process code that is transparent and freely available for download by any interested party¹.

Fig. 1 Schematic representation of the functionality of r-Java 2.0. See Table 1 for description of symbols.

The purpose of this article is to introduce the second version of r-Java (r-Java 2.0), discuss the new features that it contains, and display a selection of simulation results. Prior to delving into the details of r-Java 2.0, it would be enlightening to briefly examine the capabilities and limitations of the first version of the code. As discussed in Charignon et al. (2011), our aim for r-Java was to create an easy-to-use, cross-platform r-process code that avoids the “black-box” pitfalls that plague many scientific codes. To achieve this goal, r-Java was developed with an intuitive graphical user interface (GUI) and we provide extensive documentation and user tutorials on our website. The first version of r-Java was predicated upon the waiting point approximation (WPA), which assumes an equilibrium between neutron captures and photo-dissociations. The practical effect of using the WPA is that the relative abundance along isotopic chains depends only on neutron density (n_n), temperature (T), and neutron separation energy (S_n) (For more details see Charignon et al. 2011, Eq. (12)). Through the WPA, the number of coupled differential equations that must be solved is reduced from thousands to about a hundred (the number of isotopic chains in the reaction network), greatly reducing computational costs.

Using the WPA comes at the expense of generality, in that the assumption is only valid in high-temperature and neutron density environments (typically only considered for T₉> 2 ,where T₉ is in units of 10⁹ K and n_n> 10²⁰ cm^-3). Within the context of the WPA r-Java 1.0 provided users the ability to run r-process simulations for a wide variety of scenarios, including neutron irradiation of static targets and the dynamic expansion of r-process sites. A key feature of r-Java that has been maintained through the release of the second version is the ability for users to easily make changes to nuclear inputs between simulation runs. Other flexibilities of r-Java include the ability to specify the amount of heating from neutrinos and turn on/off the various processes and to choose the density profile of the expanding material.

This paper is organized as follows. Section 2 provides an overview of the advancements made in r-Java 2.0 and discusses the default rates included with the code. The nuclear statistical equilibrium (NSE) module is discussed in Sect. 3. The new fission methodology is detailed in Sect. 4. The effect of β-delayed neutron emission is discusssed in Sect. 5. Section 6 compares the full reaction network calculation with the WPA approach. Comparative analysis between r-Java 2.0 and other full reaction network codes is carried out in Sect. 7. Finally in Sect. 8 a summary is provided along with a look to future work on r-Java.

2. Overview of r-Java 2.0

There are major developments made to r-Java since the original release in 2011. Most noteworthy is that r-Java 2.0 is now capable of solving a full reaction network and is no longer solely reliant on the WPA. The latest release also contains a more accurate handling of fission, as well as the implementation of β-delayed neutron emission of up to three neutrons. Another expansion to r-Java is the ability for the user to specify the astrophysical environment of the r-process, which determines the methodology used to evolve the density and temperature. Furthermore, the nuclear statistical equilibrium (NSE) module has been expanded to include the effect of Coulomb screening. Finally any nuclear reaction can be turned on or off, which allows the user to investigate individual processes. An organizational chart displaying the functionality of r-Java 2.0 can be seen in Fig. 1.

As seen in Fig. 1, r-Java 2.0 contains several distinct modules. NSE, WPA network, and the full network constitute the research modules that are used by scientists to study the r-process. Complementary components to r-Java 2.0 are the teaching modules, aimed for use in the classroom at the graduate and undergraduate levels, these modules allow for investigating individual nucleosynthesis processes. The fission module calculates the mass fragmentation of neutron-induced fission. The user is able to choose the target nucleus (or nuclei), vary the incident neutron energy, and adjust four parameters related to the potential energy of the fragmentation channels (to be discussed in Sect. 4). The remaining teaching modules; β-decay, α-decay, photo-disocciation, and neutron capture all act in similar manners. The user can specify an initial abundance of nuclei and then investigate how varying the rates affects the final abundances for different physical conditions.

2.1. Rates

The default rates and cross-sections are based on the Hartree-Fock-Bogolubov 21 (HFB21) mass model (Samyn et al. 2002) as calculated by the reaction code talys (Goriely et al. 2008). The publicly available Maxwellian-averaged neutron capture cross-sections and corresponding photo-dissociation rates are provided on a temperature grid that extends from 10⁶ K to 10¹⁰ K (Goriely et al. 2008). Because neutron capture cross-sections and photo-dissociation rates can change by many orders of magnitude between temperature grid points, a simple cubic spline interpolation was insufficient. A unique interpolation method was developed that does not fall victim to the overshooting and correction of a normal cubic spline. To avoid adding uncertainty by extrapolating the photo-dissociation rates and neutron capture cross-sections, the extremel values of the temperature grid provide the temperature bounds for the r-process calculations in r-Java 2.0.

The β⁻ decay half-lives and probability of β-delayed neutron emission of up to three neutrons are considered in r-Java 2.0 by making use of the calculations by Möller et al. (2003). For complete consistency, these rates should be calculated using the HFB21 mass model, however to our knowledge such a calculation has yet to be carried out. Alpha decay half-lives are calculated based on an empirical formula dependent on the ejected alpha particle kinetic energy (Lang 1980). Figure 2 shows a schematic representation of all the processes that are incorporated in the full reaction network calculation.

Fig. 2 Reactions incorporated in the full reaction network calculation in r-Java 2.0 are represented schematically for a given isotope (Z,A). These reactions: neutron-capture (n,γ), photo-dissociation (γ,n), β-decay, beta-delayed neutron emission (βdn), α-decay, and fission.

Fig. 3 Fastest rate plotted given a temperature of 1 × 109 K and a neutron density of 1 × 1030 cm-3. The contour lines indicate when the probability of β-delayed emission of n neutrons reaches 50%. The neutron drip line and the locations of the proton and neutron magic numbers are denoted with black solid lines. The location of the stable nuclei are denoted by the black squares. A color version of this figure is available in the online article.

Fig. 4 Same as Fig. 3 but with a temperature of 3 × 109 K and a neutron density of 1 × 1020 cm-3. A color version of this figure is available in the online article.

Fig. 5 Same as Fig. 3 but with a temperature of 1 × 109 K and a neutron density of 1 × 1020 cm-3. A color version of this figure is available in the online article.

Since r-Java 2.0 makes use of temperature-dependent neutron capture cross-sections, photo-dissociation rates and neutron-induced fission cross-sections at any given temperature, and neutron density, the dominant transmutational process for each nuclei could be different. Figures 3–5 consider all available processes in r-Java 2.0 and display the dominant process for each nuclei in our network at three different neutron density and temperature combinations; Fig. 3 – ($\log {}^{\mathrm{(}}\mathit{n}_{\mathrm{n}}{}^{\mathrm{)}}\mathrm{=}\mathrm{30}\mathit{,}{\mathit{T}}_{\mathrm{9}}\mathrm{=}\mathrm{1}$), Fig. 5 – ($\log {}^{\mathrm{(}}\mathit{n}_{\mathrm{n}}{}^{\mathrm{)}}\mathrm{=}\mathrm{20}\mathit{,}{\mathit{T}}_{\mathrm{9}}\mathrm{=}\mathrm{1}$) and Fig. 4 – ($\log {}^{\mathrm{(}}\mathit{n}_{\mathrm{n}}{}^{\mathrm{)}}\mathrm{=}\mathrm{20}\mathit{,}{\mathit{T}}_{\mathrm{9}}\mathrm{=}\mathrm{3}$). For the high neutron density, low temperature scenario shown in Fig. 3 neutron capture or neutron-induced fission are the dominant channels for most nuclei in the network. Photo-dissociation is the strongest rate only for the most neutron-rich isotopes of each element. The waiting points at the N = 82, 126 and 184 closed shells can be seen as steps along the interface between neutron capture and photo-dissociation the latter being the dominant process. In the high temperature, low neutron density case displayed in Fig. 4 photo-dissociation becomes the most probable channel for the majority of nuclei in the network. α-decay dominates neutron-induced fission for some neutron-poor isotopes of heavy elements, in particular in a region around the N = 126 closed shell. For the low neutron density, low temperature example shown in Fig. 5β-decay is the dominant channel for a band of nuclei that stretches nearly the entire length of the network. Odd-even effects can be seen as along many isotopic chains β-decay and photo-dissociation alternate as the dominant process. The region in which α-decay dominates is similar to but more robust than the high temperature low neutron density case seen in Fig. 4. The fertile and fissile regime is dominated by neutron-induced fission except in a region in which spontaneous fission is the dominant decay channel. This region of spontaneous fission instability can also be seen in Fig. 4.

A fundamental tenet followed during the development of r-Java 2.0 was to maximize the flexibility afforded to the user. To this end, built into the r-Java 2.0 interface is a module dedicated to displaying and editing the nuclear parameters. The user of r-Java 2.0 can modify any parameter (mass or β-decay rate for instance) in between simulation runs, without having to restart the program. This allows users of r-Java 2.0 to quickly and easily test the effect of changing nuclear properties on r-process abundances. The choice of Java as the language for developing our nucleosynthesis code was made to ensure that r-Java 2.0 could be used across all platforms. Special attention was paid to designing a graphical user interface that is intuitive and easy-to-use.

2.2. Getting started with r-Java 2.0

As discussed above, maximizing flexibility was paramount when developing r-Java 2.0. This extends beyond the nuclear inputs and to the astrophysical parameters that govern how the temperature and density of the system evolve. With r-Java 2.0 we endeavoured to create r-process software that could be applied to any potential astrophysical r-process site. For this purpose the user of r-Java 2.0 can choose from a set of astrophysical sites which provide unique density evolutions and related input parameters. The choices for astrophysical sites are high-entropy winds around a proto-neutron star, and ejecta from a neutron star merger or the ejecta from a quark nova. The details of the specific physics implemented for each of the astrophysical sites will be discussed in a forthcoming paper. If one chooses to go beyond the aforementioned astrophysical sites, a custom density evolution can be selected. The user of r-Java 2.0 is free to define any dynamical evolution for the density or choose a static r-process site. For the remainder of this work we will consider custom density evolution equations.

With the nuclear and physical parameters chosen and before an r-process simulation can be run the user must determine the initial abundances of the r-process site. This is handled differently for the WPA and full network modules. When entering a WPA simulation, the user may specify the initial electron fraction (Y_e) and element, then based on this information and the initial temperature r-Java 2.0 computes the starting neutron density (n_n,0) and isotopic abundances using Maxwell-Boltzmann statistics. If the user chooses a full network simulation the initial mass fractions must be specified after which r-Java 2.0 calculates the starting Y_e,0 and n_n,0 ensuring that baryon number and charge are conserved.

Once the initial condition are determined, the r-process code follows the algorithm detailed in Fig. 6 and runs until the user-specified duration is met or one of the stopping criterion is satisfied. The minimum temperature stopping criterion is determined by the neutron capture and photo-dissociation temperature grid.

Fig. 6 Schematic representation of a single timestep in our full network code. { Y0(Z1,A1),Y0(Z2,A2),... } denotes the set of initial nuclei abundances, ρ0 is the initial mass density, ρ(t) defines the density evolution, T0 is the initial temperature, Ye,0 denotes the initial electron fraction, and nn,0 is the initial neutron density. First the neutron decay is computed before the reaction network is solved using the Crank-Nicholson algorithm. Next the fission contribution is calculated along with the new physical parameters. If the changes in abundance or nn are too large, the timestep is reattempted with dt = 0.1dt. Adaptive timesteps are used to maximize dt.

3. Nuclear statistical equilibrium

This release of r-Java includes a refinement to the NSE module whereby the user can now choose to include the effect of Coulomb screening. Under NSE the nuclei abundances are uniquely determined by three parameters; Y_e, mass density (ρ), and T. In the most conventional sense, when a system that follows Maxwell-Boltzmann statistics is said to be in NSE the particle number density of nuclei i, which contains Z protons and N neutrons (where mass number A = Z + N) is given by (e.g. Pathria 1977) ${\mathit{n}}_{\mathit{i}}\mathrm{=}{\mathit{g}}_{\mathit{i}}{\left(\frac{\mathrm{2}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathit{T}}{{\mathit{h}}^{\mathrm{2}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\exp \left(\frac{{\mathit{\mu }}_{\mathrm{i}}\mathrm{+}{\mathrm{B}}_{\mathrm{i}}}{\mathrm{k}\hspace{0.17em}\mathrm{T}}\right)$(1)where T represents the temperature of the system, k is Boltzmann’s constant, h Planck’s constant, and B_i, μ_i, g_i denote the binding energy, chemical potential, and statistical weight, respectively. When the Coulomb correction is applied, μ_C,tot is added in the exponential where ${\mathit{\mu }}_{\mathrm{C}\mathit{,}\mathrm{tot}}\mathrm{=}\mathit{Z}\hspace{0.17em}{\mathit{\mu }}_{\mathrm{C}\mathit{,}\mathrm{p}}\mathrm{-}{\mathit{\mu }}_{\mathrm{C}}\mathrm{\left(}\mathit{Z,A}\mathrm{\right)}\mathit{.}$(2)This correction to the chemical potential arises from the Coulomb contribution to the free energy, which becomes significant for heavier nuclei (μ_C,p is the Coulomb potential of a bare proton). Our methodology for calculating μ_C(Z,A) is similar to that of Goriely et al. (2011) and is given by ${\mathit{\mu }}_{\mathrm{C}}\mathrm{\left(}\mathit{Z,A}\mathrm{\right)}\mathrm{=}\mathit{k}\hspace{0.17em}\mathit{T}\hspace{0.17em}{\mathit{f}}_{\mathrm{C}}\mathrm{\left(}{\mathrm{\Gamma }}_{\mathit{i}}\mathrm{\right)}$(3)where f_C(Γ_i) is the Coulomb free energy per ion in units of kT. For a Coulomb liquid, f_C(Γ_i) can be expressed as (Haensel et al. 2007) $\begin{array}{ccc}{\mathit{f}}_{\mathrm{C}}\mathrm{(}{\mathrm{\Gamma }}_{\mathit{i}}\mathrm{)}& \mathrm{=}& {\mathit{A}}_{\mathrm{1}}\hspace{0.17em}\sqrt{{\mathrm{\Gamma }}_{\mathit{i}}\mathrm{(}{\mathit{A}}_{\mathrm{2}}\mathrm{+}{\mathrm{\Gamma }}_{\mathit{i}}\mathrm{)}}\mathrm{-}{\mathit{A}}_{\mathrm{1}}\mathrm{\times }{\mathit{A}}_{\mathrm{2}}\ln \left(\sqrt{{\mathrm{\Gamma }}_{\mathit{i}}\mathit{/}{\mathit{A}}_{\mathrm{2}}}\mathrm{+}\sqrt{\mathrm{1}\mathrm{+}{\mathrm{\Gamma }}_{\mathit{i}}\mathit{/}{\mathit{A}}_{\mathrm{2}}}\right)\\ & & \mathrm{+}\mathrm{2}{\mathit{A}}_{\mathrm{3}}\left(\sqrt{{\mathrm{\Gamma }}_{\mathit{i}}}\mathrm{-}\arctan \left(\sqrt{{\mathrm{\Gamma }}_{\mathit{i}}}\right)\right)\mathrm{+}{\mathit{B}}_{\mathrm{1}}\left({\mathrm{\Gamma }}_{\mathit{i}}\mathrm{-}{\mathit{B}}_{\mathrm{2}}\ln \left(\mathrm{1}\mathrm{+}\frac{{\mathrm{\Gamma }}_{\mathit{i}}}{{\mathit{B}}_{\mathrm{2}}}\right)\right)\\ & & \mathrm{+}\frac{{\mathit{B}}_{\mathrm{3}}}{\mathrm{2}}\ln \left(\mathrm{1}\mathrm{+}\frac{{\mathrm{\Gamma }}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{B}}_{\mathrm{4}}}\right)\end{array}$(4)with A₁ = −0.9070, A₂ = 0.62954, A₃ = 0.27710, B₁ = 0.00456, B₂ = 211.6, B₃ = −0.0001, B₄ = 0.00462. When a user chooses to include the Coulomb correction, r-Java 2.0 will only do so if the Coulomb liquid approximation is valid, which is to say that the Coulomb coupling parameter, ${\mathrm{\Gamma }}_{\mathit{i}}\mathrm{=}\frac{{\mathit{Z}}^{\mathrm{2}}\hspace{0.17em}{\mathit{e}}^{\mathrm{2}}}{{\mathit{a}}_{\mathit{i}}\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathit{T}}$(5)where a_i is the ion-sphere radius, is smaller than the melting value Γ_m = 175.0 ± 0.4 (Potekhin & Chabrier 2000).

Fig. 7 NSE abundance distribution subject to the following physical conditions: temperature of 1 × 1010 K, mass density of 2 × 1011 g cm-3, and electron fraction of 0.3. The red dashed line denotes a calculation that includes the effects of Coulomb interactions, while for the black solid line the Coulomb interactions were ignored.

The effect of including Coulomb screening can be seen in Fig. 7, which displays an overlay of two NSE abundances both considering the same temperature (T = 1 × 10¹⁰ K), mass density (ρ = 2 × 10¹¹ g cm^-3), and electron fraction (Y_e = 0.3), the only difference being whether Coulomb screening is included. Coulomb screening can allow for the formation of a significant number of heavier elements. As the example in Fig. 7 shows, with the Coulomb correction to the chemical potential included, a peak appears at approximately A = 124, which is absent in the case where Coulomb screening is ignored.

Since the r-process requires an explosive astrophysical site, there is a likelihood that the material that will undergo r-process will have begun in NSE (Goriely et al. 2011). To accommodate such r-process scenarios, r-Java 2.0 gives the user the option of running the NSE module and setting the resulting nuclei abundance as the initial abundance for an r-process simulation. Currently under development is a charged particle reaction network module that will be incorporated in a future release of r-Java.

4. Fission

The previous release of r-Java instituted a simple maximum Z and A approach to fission. After the reaction network was solved, species with higher values of Z or A than the imposed limit were split into two smaller species (see Charignon et al. 2011, for more details). For r-Java 2.0, the users are given the option of turning off fission, choosing the same cut-off approach as in the previous r-Java release, or choosing a more realistic treatment that includes spontaneous, neutron-induced, and β-delayed fission. Spontaneous fission rates are computed using the logic presented by Kodama & Takahashi (1975), and the β-delayed fission probabilities were taken from Panov et al. (2005). Fission barrier heights and neutron-induced fission rates provided as defaults in r-Java 2.0 are calculated by Goriely et al. (2009) based on the HFB14 mass model.

In r-Java 2.0 for the full fission treatment, three mass fragmentation channels are considered and neutron evaporation is explicitly handled for each fission event. The probability that the fission will follow a symmetric scission or one of the two standard channels is determined by integrals over the level density up to the available energy at the saddle point (Benlliure et al. 1998; Schmidt & Jurado 2010). The first standard channel (SI) results in the heavier fission fragment containing 82 neutrons, and for the second standard channel (SII) the heavier fission fragment contains approximately 88 neutrons. The likelihood of the fission event following a particular standard channel is parameterized by the relative strength (CI and CII) and depth (δVI and δVII) of the corresponding valleys in the potential energy landscape at scission. For r-Java 2.0 the strength and depth of the standard channel parameters are found through fitting observed fission fragmentation distributions for a range of nuclei between ²³²Th and ²⁴⁸Cm (Chadwick et al. 2006). The remaining fissile and fertile nuclei use the standard channel parameter values of ²³⁵U as the default values, however these parameters can be adjusted using the fission module of r-Java 2.0. The mass fragmentation distributions for ²³²Th, ²³⁵U, and ²⁴⁰Pu are displayed in Fig. 8. For reference in Fig. 8, the results calculated using the fission module in r-Java 2.0 are compared to the results of the GEF model (Schmidt & Jurado 2010), as well as observations (Chadwick et al. 2006).

Fig. 8 Top: fission fragment mass distribution resulting for neutron-induced fission of 232Th by 1.5 MeV neutrons. Middle: fission fragment mass distribution resulting for neutron-induced fission of 235U by 1.5 MeV neutrons. Bottom: fission fragment mass distribution resulting for neutron-induced fission of 240Pu by 1.5 MeV neutrons.

Fig. 9 Comparison of the full fission methodology to the mass cut-off approach. Two different initial neutron-to-seed ratios (top: Yn/Yseed ~ 137, bottom: Yn/Yseed ~ 186) are considered while all other parameters remain the same (see section 4 of text for details). In both panels the red dashed line denotes the final abundance of a simulation that used the mass cut-off approach, while the black solid line represents the full fission treatment. The relevant magic numbers are highlighted with a fine vertical black line.

Fig. 10 Overlay of the abundances after having allowed the system to decay back to stability (black solid line) and at the end of the r-process (red dashed line) for the same two initial neutron-to-seed ratios simulations shown in Fig. 9. Top: Yn/Yseed = 137. Bottom: Yn/Yseed = 186. The relevant magic numbers are highlighted with a fine vertical black line.

The result of the mass fragmentation calculation for each fissionable parent is that a probability distribution of potential daughter pairs is found. In r-Java 2.0 the probability for each daughter species is multiplied by the parent fission rate and is incorporated as the daughter production rate in the network calculation.

A comparison of r-process final abundance distributions for both the mass cut-off and full fission treatment can be seen in Fig. 9. For the cut-off fission methodology, a maximum mass of A = 272 was used, and each fissioning nuclei splits into two daughter species. The full fission treatment uses the three fission processes discussed above, as well as the fission fragmentation calculation. For each simulation displayed in Fig. 9, the same initial abundance of iron-group nuclei were used, starting from the same initial temperature of 1.0 × 10⁹ K. The initial mass density was 10¹¹ g cm^-3 for each simulation run, which followed the same density profile, ρ(t) = ρ₀/ (1 + 1.5 t/τ)² with an expansion timescale (τ) of 0.003 s. The only variation between simulation runs shown in the top and bottom panels of Fig. 9 is the neutron-to-seed ratio (Y_n/Y_seed). For the top panel Y_n/Y_seed = 137 was used and the bottom panel displays the r-process yield of a more neutron-rich simulation run, which began with Y_n/Y_seed = 186. These two parameter sets were chosen to highlight the differences between the two fission methodologies.

In the smaller Y_n/Y_seed scenario, shown in the top panel of Fig. 9, the r-process is just capable of breaking through the N = 184 magic number. In this case the full fission run has been able to cross over a region of instability at about A ~ 280 and has produced a small peak of super-heavies at about A ~ 290. Aside from the small super-heavy peak, the results of both the full fission and cut-off methodology are largely the same.

For the larger Y_n/Y_seed scenario seen in the bottom panel of Fig. 9, the fission cut-off approach overproduces nuclei at A ~ 130 by over ten times compared to the full fission treatment. This overproduction is due to the increased fission recycling caused by forcing all nuclei heavier than A ~ 272 to undergo fission. For this Y_n/Y_seed, the full fission simulation run produces a super-heavy peak of the order of 10^-7. The results of these simulation runs do not speak to the long-term stability of the super-heavy nuclei produced but rather shows the large variation between the two fission methodologies at the point of neutron freeze-out, which is to say that when the neutron to r-process product ratio drops below one, (Y_n/Y_r< 1).

The final abundances once the systems are allowed to decay to stability can be seen in Fig. 10. Fission recycling gives rise to nearly all the nuclei abundances below A ~ 150 seen in both cases. The distribution of fission recycled nuclei is similar in both cases because fission is occurring from the same region. The shape of the fission contribution found using r-Java 2.0 coincides with the findings of Petermann et al. (2008), who used the statistical code ALBA to calculate the fission yield for each fission event. For comparison, the abundances at the moment r-process stops is included in Fig. 10 for both neutron-to-seed simulation runs. The robust fission calculations included in r-Java 2.0 provides an accurate assessment of the role of fission recycling in the r-process.

Fig. 11 Effect of β-delayed neutron emission on nuclei abundance. The black line denoting an r-process simulation that included β-delayed neutron emission, and for the red dashed line that process was omitted. The results plotted in this figure, as well as in Figs. 12 and 13, are from simulation runs that were identical with the exception of whether or not β-delayed neutron emission was included. Top: the nuclei abundances at the moment the neutron-to-seed ratio drops below one. Bottom: the nuclei abundances after decay to stability. The relevant magic numbers are highlighted with a fine vertical black line.

5. Beta-delayed neutron emission

To study the effects of β-delayed neutrons on the r-process, we compare two simulation runs that are identical except for whether β-delayed neutron emission is included. The top panel of Fig. 11 displays a comparison of the abundance distributions at the end of the r-process, which for this study was defined to be once the neutron-to-r-process products ratio (Y_n/Y_r) drops below one. The emission of β-delayed neutrons acts to smooth out the variability in nuclei distribution compared to that of the case without β-delayed neutrons. The peak at A ~ 188 is shifted slightly heavier with the inclusion of β-delayed neutrons and also the abundance of nuclei with mass greater than A = 200 is increased. The lower panel of Fig. 11 shows the final nuclei abundance distribution once the systems are allowed to decay to stability. For the simulation that did not include β-delayed neutron emission, the nuclei abundance distribution below A ≃ 209 remains virtually unchanged from the time r-process stops to that of stability. However, when β-delayed neutrons are included, the decay to stability causes further reduction in the variability of the abundance distribution, and a shifting of the peaks towards lower mass.

Fig. 12 Evolution of neutron density until the r-process is terminated. The black line denotes an r-process simulation that included β-delayed neutron emission, and for the red dashed line that process was omitted.

In Fig. 12 the evolution of neutron density during the r-process is compared between the simulations with and without β-delayed neutron emission. The β-delayed neutrons act to keep the neutron density higher for longer then in the case in which β-delayed neutron emission was ignored. By bolstering the neutron density, β-delayed neutron emission allows the r-process to proceed more readily to heavier elements, an effect that can can be seen in the top panel of Fig. 11.

The abundances of nuclei at the end of the r-process plotted on the (N,Z) plane can be seen in Fig. 13 (top panel displays the case where β-delayed neutron emission was ignored and the bottom panel the case with its inclusion). For the simulation run that included β-delayed neutrons, the r-process accesses a broader range (along lines of constant Z) of nuclei, reaching closer to the valley of stability. This broadening effect caused by β-delayed neutrons is most noticeable around the N = 82 and 126 closed shells. The ability of β-delayed neutron emission to allow for matter flow past the N = 126 closed shell can be seen in Fig. 13 as the breadth of populated nuclei and abundance in the region past N = 126 is increased in the case where β-delayed neutron emission is included.

6. Full reaction network

To expand beyond the WPA, reactions that stay within an isotopic chain, namely neutron-capture and photo-dissociation, must be included in the network calculation. This means that rather than solving a system of equations the size of which is determined by the number of isotopic chains (110) as in the WPA case, for the full network case an equation for every nuclei must be included (a total of 8055). The computational cost of this addition is significant since finding a solution to a reaction network scales as N³ where N is the number of coupled differential equations. However there are methods that can be invoked to mitigate this cost; we take advantage of the fact that each nuclei in the network is only coupled to another nuclei if there is an adjoining reaction (i.e. nuclei (Z,A) is coupled to both (Z+2, A+4) and (Z-2, A-4) via α decay). This is effectively utilizing the sparseness of the reaction rate matrix, which alleviates memory load issues and speeds up runtime. We solve the fully implicit network using the Crank-Nicholson method. The rate of thermonuclear energy released (or absorbed) is calculated using the methodology laid out by Hix & Meyer (2006).

Fig. 13 Nuclei abundances at the moment the neutron-to-seed ratio drops below one plotted on the (N,Z) plane. Stable nuclei, the location of the proton and neutron closed shells, and the neutron drip line are included for reference. Top: simulation that did not include β-delayed neutron emission. Bottom: simulation including β-delayed neutron emission.

Figure 14 highlights the importance of the imposed stopping criteria on network calculations through a comparison of the results from the WPA to that of the full network. For the results plotted in both panels of Fig. 14, the same initial conditions and expansion profiles² were chosen. The simulations considered begin from an iron seed with Y_e,0 = 0.16, T₀ = 4 × 10⁹ K, and ρ₀ = 10¹⁰ g cm^-3. In the top panel of Fig. 14, the calculations are stopped when the temperature falls below 2 × 10⁹ K, an imposed cut-off based on the work of Cowan et al. (1983). The nuclei distribution in the WPA simulation is peaked at A = 80 with lower abundances of nuclei up to A ~ 120 and then a precipitous drop in abundance for heavier nuclei. In the case where the full reaction network calculation was stopped once the temperature fell to 2 × 10⁹ K (displayed in the top panel of Fig. 14), there is good agreement to the WPA calculation. The shape of the nuclei abundance distribution is the same for both network calculations, with the full reaction network producing a slightly greater abundance of heavy nuclei. However, the results displayed in the top panel of Fig. 14 are not indicative of the full potential of the r-process for this chosen environment, since as the temperature drops below the imposed minimum cut-off, the neutron density still remains high (n_n ~ 10³⁰ cm^-3). For the bottom panel of Fig. 14, the minimum temperature stopping criterion for the WPA was lowered to 10⁹ K, and for this case both the WPA and full reaction network calculations halt at neutron freeze-out (Y_n/Y_r = 1). Once again both network calculations display similar nuclei abundance distributions. The results of the WPA reflect the (n,γ) ⇋ (γ,n) equilibrium, which is not as accurate as the full treatment. This is manifested as deeper troughs in nuclei abundance, especially around the A = 190 peak, and greater variability for the lower mass nuclei. The smoother distribution in the full reaction network results is also due to the inclusion of β-delayed neutron emission.

The users of r-Java 2.0 are afforded the option of choosing the stopping criteria for r-process calculations: minimum temperature and neutron density for the WPA network and Y_n/Y_r for both networks.

Fig. 14 Comparison of the simulation results from the WPA (red dashed line) with that of the full network (black solid line). Top: nuclei abundances when the temperature drops to 2 × 109 K. Bottom: the nuclei abundances at neutron freeze-out, see text for details of stopping criteria. The relevant magic numbers are highlighted with a fine vertical black line.

7. Test cases

As part of the testing phase of the development of r-Java 2.0 we attempted to reproduce the results from three other full network r-process codes; the Clemson University nucleosynthesis code (Jordan & Meyer 2004) which will furthermore be referred to as the Clemson code, and the Basel University nucleosynthesis code (Freiburghaus et al. 1999a), to be referred to as the Basel code and the nucleosynthesis code developed at Université Libre de Bruxelles (Goriely et al. 2011) which will be called the Bruxelles code for the remainder of this article. While a complete apples-to-apples comparison was not tenable, the results of our tests showed good agreement with each of the three codes studied.

7.1. Clemson nucleosynthesis code

For the Clemson code comparison seen in Fig. 15, we endeavoured to reproduce the results shown in Figs. 7 and 8 of Jaikumar et al. (2007). We found that the initial abundance was not very important for either case because the neutron-to-seed ratio was high enough that any influence from the initial abundance was washed away by the r-process. The top panel of Fig. 15 shows the results of a fast expansion r-process site and the bottom panel a slow expansion (corresponding to Figs. 7 and 8 from Jaikumar et al. 2007, respectively). In the fast expansion case both r-Java 2.0 and the Clemson code show that the r-process is not capable of proceeding past the A = 130 magic number. In the slow expansion case, the environment remains favourable for the r-process much longer, and the final abundance for both r-Java 2.0 and the Clemson code contains peaks shifted to the heavy side of the A = 130 and A = 190 observed solar peaks. The differences between the final abundances from r-Java 2.0 and the Clemson code seen in both cases can be credited to the fact that the two codes use different mass models (the Clemson code used the finite range droplet model and r-Java 2.0 HFB21), which has been shown to affect the r-process abundance yield (e.g. Farouqi et al. 2010).

Fig. 15 Top: comparison of the final abundances from r-Java 2.0 (red dashed line) and the Clemson nucleosynthesis code (black solid line) for a fast expansion r-process site. Bottom: a comparison of the final abundances from r-Java 2.0 with two different initial neutron-to-seed ratios (Yn/Yseed ~ 1100 denoted by the green dotted line, Yn/Yseed ~ 1300 by the red dashed line) and the Clemson nucleosynthesis code (black solid line) for a slow expansion r-process site. The relevant magic numbers are highlighted with a fine vertical black line.

7.2. Basel nucleosynthesis code

To compare to an updated version of the Basel code, we pushed to reproduce the abundances shown in Fig. 10 of Farouqi et al. (2010), which considers the HFB17 mass model. As described there the r-process network begins at the termination of the charged particle network, thus we used the abundance per mass number at the end of the charged-particle network displayed in their Fig. 5 to determine our initial seed nuclei distribution for comparison. Having only the abundance per mass number information, we had to choose which nuclei to set each abundance to in order to build our initial seed nuclei. We made the assumption that the system is in (n,γ) ⇋ (γ,n) equilibrium at the beginning of the r-process based on the initial conditions used in Farouqi et al. (2010) of T = 3 × 10⁹ K and n_n = 10²⁷ cm^-3. Then for abundance at each mass number plotted in Fig. 5 of their work, we set it to the isotope that most closely matched the predictions of the nuclear Saha equation. For each different entropy simulation that we ran, these abundances were uniformly scaled such that they produced the correct seed abundance as shown in Fig. 3 of Farouqi et al. (2010). The initial neutron abundance was then determined from the neutron-to-seed ratio stated in their Table 5. The use of the same initial abundance distribution for each simulation run by r-Java 2.0, which may not have been the case for their simulations, is the largest potential source of discrepancy in this comparative analysis.

Consistent with Farouqi et al. (2010), we used Y_e = 0.45 and started our simulations with an initial temperature of 3 × 10⁹ K. We followed the same constant entropy methodology described there to evolve temperature and density. In this scenario the temperature evolves adiabatically, and the entropy is assumed to be radiation-dominated, which allows for the inference of the evolution of matter density. The time dependence of the temperature and matter density (ρ₅ is in units of 10⁵ g cm^-3) are thus governed by the following equations, $\begin{array}{ccc}{\mathit{T}}_{\mathrm{9}}\mathrm{(}{\mathit{t}}^{\mathrm{)}}& \mathrm{=}& {\mathit{T}}_{\mathrm{9}}\mathrm{(}\mathit{t}\mathrm{=}{\mathrm{0}}^{\mathrm{)}}\frac{{\mathit{R}}_{\mathrm{0}}}{{\mathit{R}}_{\mathrm{0}}\mathrm{+}{\mathit{v}}_{\exp }\hspace{0.17em}\mathit{t}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{5}}\mathrm{(}{\mathit{t}}^{\mathrm{)}}& \mathrm{=}& \mathrm{1.21}\frac{{\mathit{T}}_{\mathrm{9}}^{\mathrm{3}}}{\mathit{S}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{7}}{\mathrm{4}}\frac{{\mathit{T}}_{\mathrm{9}}^{\mathrm{2}}}{\mathrm{(}{\mathit{T}}_{\mathrm{9}}^{\mathrm{2}}\mathrm{+}{\mathrm{5.3}}^{\mathrm{)}}}\right)\mathit{,}\end{array}$where R₀ = 130 km and v_exp = 7500 km s^-1.

To maintain consistency with the Basel code, we terminated the r-process once the neutron-to-seed ratio dropped below one, and the abundances shown in Fig. 16 are after decay back to stability.

The top left-hand panel of Fig. 16, which displays the results of the S = 175 simulation runs, shows the best agreement between r-Java 2.0 and the Basel code of all the cases tested. Both r-Java 2.0 and the Basel code show a final nuclei abundance that predominantly ranges from 70 <A< 135. The results from each code displays a peak below the A = 130 magic number, however the Basel code peak is shifted heavier with respect to that of r-Java 2.0.

Fig. 16 Comparison of r-process abundance yields as calculated by r-Java 2.0 (black solid line) and the Basel nucleosynthesis code (red dashed line). The relevant magic numbers are highlighted with a fine vertical black line. For each panel a different entropy was assumed, which changes the initial neutron-to-seed ratio as well as the evolution of the density, see text for details. Top-left: simulation run assuming the entropy of the wind is S = 175. Top-right: simulation run assuming the entropy of the wind is S = 195. Bottom-left: simulation run assuming the entropy of the wind is S = 236. Bottom-right: simulation run assuming the entropy of the wind is S = 280. See text for details of initial conditions.

The S = 195 simulation results (displayed in the top right-hand panel of Fig. 16) from the Basel code and r-Java 2.0 are both dominated by a peak at the A = 130 magic number. The differences between the final abundances from the two codes for this entropy are consistent with differing initial abundances. That the results from r-Java 2.0 display a more distinct peak at the A = 80 magic number is consistent with the simulation run of r-Java 2.0 starting with more nuclei below the A = 80 magic number. This would lead to nuclei piling up at A = 80 for r-Java 2.0, which would not be the case for the Basel code simulation run. The difference in initial abundance also has an effect on the heavy side of the final abundance distribution. With more nuclei initially between the A = 80 and A = 130 observed solar peaks, the r-process simulation run of the Basel code is more capable of pushing through the A = 130 magic number to higher masses. As for the r-Java 2.0 simulation, once the r-process pushes through the A = 80 magic number, nuclei will pile up on the light-side of the A = 130 magic number. By the time the r-process reaches the A = 130 peak in the r-Java 2.0 run, the neutron density will have dropped too low to significantly push past the A = 130 magic number. The result of this is the increased production of nuclei on the lower mass side of the A = 130 peak for the r-Java 2.0 simulation run with respect to that of the Basel code and a longer high-mass tail in the Basel code simulation.

Similar to the S = 195 case, the presence of nuclei below the A = 80 magic in the S = 236 r-Java 2.0 simulation (seen in the lower-left panel of Fig. 16) leads to the final abundance containing nuclei around A = 80, which is not the case for the Basel code results. Once again this difference is consistent with different initial abundances for the two runs. Neglecting the relatively small abundance for 80 ≲ A ≲ 125 in the r-Java 2.0 results, the final distributions of the S = 236 simulation runs for both codes are consistent with peaks around A = 80, 165, and 190. The A = 190 peak in the Basel code simulation is stronger and shifted towards heavier masses with respect to that of r-Java 2.0, which can be attributed to fact that the r-Java 2.0 simulation had more nuclei stuck below the A = 80 magic number.

The S = 280 simulation runs seen in the lower right-hand panel of Fig. 16 shows the same basic features for both codes. The final nuclei abundance for both codes contains strong peaks at A = 130 and A = 195, with the r-Java 2.0 results displaying stronger peaks. The increased abundance of Th and U at stability in the Basel code simulation run could be due to the initial abundance differences discussed for the S = 195 and S = 236 cases or due to different definitions of stability. For the r-Java 2.0 simulations, the systems decayed for 13 Gyr or until the percent change in any nuclei abundance was less than 1 × 10^-15.

7.3. Université Libre de Bruxelles nucleosynthesis code

For our comparison to the Bruxelles code we attempted to reproduce the abundances after the decompression displayed in Fig. 10 of Goriely et al. (2011). As discussed in Goriely et al. (2011), the initial abundances used for the r-process simulation are important because in this scenario the initial neutron-to-seed ratio is roughly 5, and the r-process is only capable of shifting the abundances toward heavier nuclei without dramatically altering the relative shape of the abundance distribution. Goriely et al. (2011) provide the initial abundances used for the r-process simulation, which are calculated under NSE with Coulomb interactions included. A comparison of the initial abundances of the Bruxelles code and r-Java 2.0 can be seen in the top panel of Fig. 17. The peaks roughly centred at A = 80 and 125 as calculated by r-Java 2.0 are higher than for the Bruxelles code, while the intermediate-mass region is more abundant in the Bruxelles calculation. The NSE calculation performed by r-Java 2.0 assumes Maxwell-Boltzmann statistics, while the Bruxelles code used Fermi-Dirac, which accounts for the differences in abundances. While the nuclear physics used in the Bruxelles code is the most similar to r-Java 2.0 of all the codes studied, we had to implement an analytic approximation to the density evolution used by Goriely et al. (2011). To compare to the Bruxelles code, we chose the density profile shown in Eq. (8) $\mathit{\rho }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{(}\mathit{a}\hspace{0.17em}\mathit{t}\mathit{/}{\mathit{\tau }}^{\mathrm{)}}}^{\mathit{b}}}\right)}^{\mathit{c}}\mathit{,}$(8)where a, b, and c are free parameters. A value of 3 × 10^-4 s was used for the expansion timescale (τ), which is consistent with the one used by Goriely et al. (2011).

Fig. 17 Top: initial abundances used for the comparison of r-process simulations from r-Java 2.0 and the Bruxelles nucleosynthesis code. The red dashed line denotes r-Java 2.0 and the black solid line the Bruxelles code. Bottom: final abundances from r-Java 2.0 considering two different density evolution profiles (red dashed line and green dotted line) compared to that of the Bruxelles code (black solid line). See text for details of simulations. The relevant magic numbers are highlighted with a fine vertical black line.

Fig. 18 Top: final (black solid line) and initial (red dashed line) abundances as calculated by r-Java 2.0 for comparison to the Bruxelles code. Bottom: final (black solid line) and initial (red dashed line) abundances as calculated by the Bruxelles code. See text for details of simulations. The relevant magic numbers are highlighted with a fine vertical black line.

A comparison of two different sets of free parameters used in the density profile of r-Java 2.0 to the final abundances of the Bruxelles code can be seen in the bottom panel of Fig. 17. As expected, the differences in initial abundances are carried through to the final nuclei abundances with r-Java 2.0 displaying higher peaks at approximately A = 85 and 130 with the intermediate mass region more strongly produced in the Bruxelles code simulation. To show that in both codes the r-process has the same effect on abundances in Fig. 18, the final and initial abundances are overplotted for each code. For both r-Java 2.0 and the Bruxelles nucleosynthesis code the r-process acts to shift the peaks towards heavier nuclei.

8. Summary and conclusions

This paper has discussed the nuclear physics incorporated in r-Java 2.0; providing cutting-edge fission calculations, β-delayed neutron emission of up to three neutrons and neutron capture and photo-dissociation rates from one of the most sophisticated mass models (HFB21). Nevertheless, it is the ability to change any parameter quickly and easily that makes r-Java 2.0 a powerful tool for studying nuclear astrophysics. r-Java 2.0 is capable of solving a full r-process reaction network containing over 8000 nuclei and can do so both accurately and efficiently, with a typical full reaction network simulation completed in minutes. The scientific aim of this release of r-Java is to study r-process nucleosynthesis in the expansion phase (T ≲ 3 × 10⁹ K, e.g. Howard et al. 1993) and NSE at high temperature (T ≳ 4 × 10⁹ K, e.g. Truran et al. 1966). We are currently developing a charged-particle reaction network module that will be incorporated into a future version of r-Java.

With a more realistic treatment of fission we have added to r-Java 2.0 the ability to investigate the role of fission recycling in the r-process. In the past by simply using the mass cut-off approach the mistake of going to too high of a neutron density was masked by the fact that fission recycling would not allow the r-process to proceed beyond the cut-off. This presents in the r-process abundance at neutron freeze-out in two ways; an underproduction of super-heavy nuclei and the overproduction of nuclei around the A = 130 magic number. With the fission methodology implemented here the super-heavy regime (A> 270) can be studied using r-Java 2.0. The preliminary study undertaken here supports the findings of Petermann et al. (2012), where super-heavy nuclei (A ~ 290) can be formed by the r-process. The super-heavies subsequently decay in seconds.

The emission of β-delayed neutrons can act to maintain a sufficiently high neutron density to allow for the r-process to reach heavier elements. The effect of β-delayed neutron emission is also significant during the decay to stability once the r-process has stopped. They act to smooth out the nuclei distribution on the path to stability and shifts the abundances to lower masses. Their role may in some cases not be as direct as just stated. The β-delayed neutrons can alter the r-process path, thereby accessing nuclei that would more readily capture neutrons and cause the neutron density to drop more rapidly than if they were ignored. This must be studied in more detail, and with r-Java 2.0 the user can quickly and easily investigate the effect of β-delayed neutrons on r-process abundances.

By performing a comparative study between r-Java 2.0 and three other full network r-process codes, we have found good agreement between the codes, however undertaking this analysis has highlighted the potential pitfalls of comparing the results from different codes. Factors such as choice of mass model, evolution methodology of physical parameters, code stopping criteria, and precision can contribute to variations in r-process abundances that are artifacts of the nucleosynthesis code structure rather than of the physical scenarios being studied. This comparative analysis highlights the universality of r-Java 2.0, which by allowing the user to customize both the nuclear and astrophysical parameters, is capable of reproducing the results of other nucleosynthesis codes.

The development of r-Java 2.0 was done in a way that maximizes the flexibility of the software, allowing for the adjustment of any nuclear or physical property both quickly and easily. The choice of Java as the programming language allowed for including an easy to use GUI that is cross-platform compatible. Beyond its applicability to scientific study, the goal of r-Java 2.0 was to make it accessible in a teaching capacity by ensuring it is easy to use and allowing for the investigation of individual processes.

In the follow-up paper to this work we will turn our attention to the astrophysical side of the r-process that is well covered by r-Java 2.0. Built into the interface of r-Java 2.0 is the option of defining a custom density evolution or of selecting one of three proposed astrophysical r-process sites: high-entropy winds around protoneutron stars (example studies: Woosley & Hoffman 1992; Qian & Woosley 1996; Thompson et al. 2001; Farouqi et al. 2010), ejecta from neutron star mergers (Freiburghaus et al. 1999b; Goriely et al. 2011, and others), or ejecta from quark novae (Jaikumar et al. 2007). For each of the proposed astrophysical sites, r-Java 2.0 consistently calculates the temperature and density evolution, the details of which will be discussed in this upcoming paper. By including the physics of different astrophysical sites in one piece of r-process software, we have provided a common platform for comparing the r-process abundances of different astrophysical sites.

Acknowledgments

This work is supported by the Natural Sciences and Engineering Research Council of Canada. NK acknowledges support from the Killam Trusts.

References

All Tables

All Figures

	Fig. 1 Schematic representation of the functionality of r-Java 2.0. See Table 1 for description of symbols.
In the text

	Fig. 2 Reactions incorporated in the full reaction network calculation in r-Java 2.0 are represented schematically for a given isotope (Z,A). These reactions: neutron-capture (n,γ), photo-dissociation (γ,n), β-decay, beta-delayed neutron emission (βdn), α-decay, and fission.
In the text

Fig. 3 Fastest rate plotted given a temperature of 1 × 109 K and a neutron density of 1 × 1030 cm-3. The contour lines indicate when the probability of β-delayed emission of n neutrons reaches 50%. The neutron drip line and the locations of the proton and neutron magic numbers are denoted with black solid lines. The location of the stable nuclei are denoted by the black squares. A color version of this figure is available in the online article.

In the text

	Fig. 4 Same as Fig. 3 but with a temperature of 3 × 109 K and a neutron density of 1 × 1020 cm-3. A color version of this figure is available in the online article.
In the text

	Fig. 5 Same as Fig. 3 but with a temperature of 1 × 109 K and a neutron density of 1 × 1020 cm-3. A color version of this figure is available in the online article.
In the text

Fig. 6 Schematic representation of a single timestep in our full network code. { Y0(Z1,A1),Y0(Z2,A2),... } denotes the set of initial nuclei abundances, ρ0 is the initial mass density, ρ(t) defines the density evolution, T0 is the initial temperature, Ye,0 denotes the initial electron fraction, and nn,0 is the initial neutron density. First the neutron decay is computed before the reaction network is solved using the Crank-Nicholson algorithm. Next the fission contribution is calculated along with the new physical parameters. If the changes in abundance or nn are too large, the timestep is reattempted with dt = 0.1dt. Adaptive timesteps are used to maximize dt.

In the text

	Fig. 7 NSE abundance distribution subject to the following physical conditions: temperature of 1 × 1010 K, mass density of 2 × 1011 g cm-3, and electron fraction of 0.3. The red dashed line denotes a calculation that includes the effects of Coulomb interactions, while for the black solid line the Coulomb interactions were ignored.
In the text

	Fig. 8 Top: fission fragment mass distribution resulting for neutron-induced fission of 232Th by 1.5 MeV neutrons. Middle: fission fragment mass distribution resulting for neutron-induced fission of 235U by 1.5 MeV neutrons. Bottom: fission fragment mass distribution resulting for neutron-induced fission of 240Pu by 1.5 MeV neutrons.
In the text

Fig. 9 Comparison of the full fission methodology to the mass cut-off approach. Two different initial neutron-to-seed ratios (top: Yn/Yseed ~ 137, bottom: Yn/Yseed ~ 186) are considered while all other parameters remain the same (see section 4 of text for details). In both panels the red dashed line denotes the final abundance of a simulation that used the mass cut-off approach, while the black solid line represents the full fission treatment. The relevant magic numbers are highlighted with a fine vertical black line.

In the text

	Fig. 10 Overlay of the abundances after having allowed the system to decay back to stability (black solid line) and at the end of the r-process (red dashed line) for the same two initial neutron-to-seed ratios simulations shown in Fig. 9. Top: Yn/Yseed = 137. Bottom: Yn/Yseed = 186. The relevant magic numbers are highlighted with a fine vertical black line.
In the text

Fig. 11 Effect of β-delayed neutron emission on nuclei abundance. The black line denoting an r-process simulation that included β-delayed neutron emission, and for the red dashed line that process was omitted. The results plotted in this figure, as well as in Figs. 12 and 13, are from simulation runs that were identical with the exception of whether or not β-delayed neutron emission was included. Top: the nuclei abundances at the moment the neutron-to-seed ratio drops below one. Bottom: the nuclei abundances after decay to stability. The relevant magic numbers are highlighted with a fine vertical black line.

In the text

	Fig. 12 Evolution of neutron density until the r-process is terminated. The black line denotes an r-process simulation that included β-delayed neutron emission, and for the red dashed line that process was omitted.
In the text

	Fig. 13 Nuclei abundances at the moment the neutron-to-seed ratio drops below one plotted on the (N,Z) plane. Stable nuclei, the location of the proton and neutron closed shells, and the neutron drip line are included for reference. Top: simulation that did not include β-delayed neutron emission. Bottom: simulation including β-delayed neutron emission.
In the text

	Fig. 14 Comparison of the simulation results from the WPA (red dashed line) with that of the full network (black solid line). Top: nuclei abundances when the temperature drops to 2 × 109 K. Bottom: the nuclei abundances at neutron freeze-out, see text for details of stopping criteria. The relevant magic numbers are highlighted with a fine vertical black line.
In the text

Fig. 15 Top: comparison of the final abundances from r-Java 2.0 (red dashed line) and the Clemson nucleosynthesis code (black solid line) for a fast expansion r-process site. Bottom: a comparison of the final abundances from r-Java 2.0 with two different initial neutron-to-seed ratios (Yn/Yseed ~ 1100 denoted by the green dotted line, Yn/Yseed ~ 1300 by the red dashed line) and the Clemson nucleosynthesis code (black solid line) for a slow expansion r-process site. The relevant magic numbers are highlighted with a fine vertical black line.

In the text

Fig. 16 Comparison of r-process abundance yields as calculated by r-Java 2.0 (black solid line) and the Basel nucleosynthesis code (red dashed line). The relevant magic numbers are highlighted with a fine vertical black line. For each panel a different entropy was assumed, which changes the initial neutron-to-seed ratio as well as the evolution of the density, see text for details. Top-left: simulation run assuming the entropy of the wind is S = 175. Top-right: simulation run assuming the entropy of the wind is S = 195. Bottom-left: simulation run assuming the entropy of the wind is S = 236. Bottom-right: simulation run assuming the entropy of the wind is S = 280. See text for details of initial conditions.

In the text

Fig. 17 Top: initial abundances used for the comparison of r-process simulations from r-Java 2.0 and the Bruxelles nucleosynthesis code. The red dashed line denotes r-Java 2.0 and the black solid line the Bruxelles code. Bottom: final abundances from r-Java 2.0 considering two different density evolution profiles (red dashed line and green dotted line) compared to that of the Bruxelles code (black solid line). See text for details of simulations. The relevant magic numbers are highlighted with a fine vertical black line.

In the text

	Fig. 18 Top: final (black solid line) and initial (red dashed line) abundances as calculated by r-Java 2.0 for comparison to the Bruxelles code. Bottom: final (black solid line) and initial (red dashed line) abundances as calculated by the Bruxelles code. See text for details of simulations. The relevant magic numbers are highlighted with a fine vertical black line.
In the text