Free Access
Issue
A&A
Volume 531, July 2011
Article Number A79
Number of page(s) 11
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201117008
Published online 16 June 2011

© ESO, 2011

1. Introduction

The rapid neutron capture process (r-process) is believed to be the mechanism for the nucleosynthesis of about half of the stable nuclei heavier than iron (Burbidge et al. 1957; Cameron 1957). Explosive and neutron-rich astrophysical environments present ideal conditions for the r-process to take place. Possible candidate sites discussed in the literature include the neutrino-driven, neutron-rich wind from proto-neutron stars (Qian & Woosley 1996; Qian 2003), prompt explosions of collapsed stellar cores, (Sumiyoshi et al. 2001; Wanajo et al. 2003), neutron star decompression (Meyer 1989; Goriely et al. 2004), tidal disruption in binary merger events (Freiburghaus et al. 1999a), and outflows in gamma-ray bursts (Surman et al. 2006) etc. Most importantly, abundance data on r-process elements in metal-poor stars (Sneden et al. 2003) and certain radionuclides in meteorites (Qian & Wasserburg 2007) point toward the distinct possibility of multiple r-process sites, so that more than one of these sites may contribute to the observed abundance pattern of the r-process elements (Truran et al. 2002).

In 1985, an important work by Bethe & Wilson brought neutrino-driven winds from type II supernovae (SNe) to the forefront of the discussion about astrophysical sites for the r-process. Since then, much progress has been made in the modelling of type II SNe and neutrino winds of nascent neutron stars. Some examples of such work are Woosley et al. (1994), Takahashi et al. (1994), Qian & Woosley (1996), Cardall & Fuller (1997), Otsuki et al. (2000), Wanajo et al. (2001), Thomson et al. (2001). However, we have not arrived at a single self-consistent astrophysical model that reproduces the observed abundance of r-process elements and that at the same time is consistent with spectroscopic and meteoritic data. The most natural explanation is that the observed r-process pattern follows from a superposition of neutron capture events with differing neutron-to-seed ratios and exposure timescales. Producing the third peak requires extreme values of entropy and a dynamic timescale that both challenge current models of type II SNe explosions (Qian & Wasserburg 2007). Neutron star mergers provide a much higher neutron-to-seed ratio than type II SNe, but occur so infrequently that they fail to explain the early enrichment of r-process elements relative to iron observed in metal-poor stars (Qian 2000). Clearly, further study is needed to better understand the r-process conditions, so to this end we present a fairly general and flexible r-process code that is transparent and freely available to the nuclear physics and astrophysics communities.

thumbnail Fig. 1

Schematic representation of r-Java’s capabilities and tools (see text for meaning of symbols). More details can be found in the manual available on the quark nova project website: http://quarknova.ucalgary.ca/.

Open with DEXTER

The motivation behind developing our own r-process code came in part from recognizing the difficulty in having access to an openly shared cross-platform program where any interested researcher can examine and modify the dynamical evolution equations and the nuclear input. We perceive a distinct need for this, as rare-isotope experiments such as FRIB later in this decade will make available new and improved nuclear data for reaction rates relevant to the r-process. Besides this advantage, the r-process is of wide interest to the astrophysics community, who might wish to explore different dynamical environments for the r-process and not be as concerned, so to speak, about the input nuclear physics. To avoid a kind of “black box” situation arising from limited exchange of numerical codes, we decided to build a new code that is flexible enough to incorporate different conditions and has detailed documentation. The code is capable of simulating nuclear statistical equilibrium (NSE), neutron bombardment of a static target, as well as the expansion (including the relativistic case) of neutron-rich material. In this work, we discuss all of the above and also apply this code to the unique environment of an ejected and expanding neutron star crust. For the code we designed a graphical user interface (GUI) which allows for real-time interpretation and analysis of simulations through data and graphs. We have made the r-process application freely available in Java-based format on the quark nova project website1, where one can also find input nuclear data and supporting documentation on most aspects of our code. We request users to test it, and welcome questions and feedback. Our aim is to make our nucleosynthesis code, which we call r-Java, and its workings as transparent as possible so that any scientific user can adapt it to their needs.

This paper is organized as follows. In Sect. 2, we mention our sources for the nuclear input within r-Java and display results for NSE that are in agreement with previous benchmarks. In Sects. 3 and 4, we present the simplified reaction network (waiting-point approximation) and simulation results for heavy-element production by static and dynamic neutron capture. In Sect. 5, the conclusions from our simulations are discussed and we look ahead to future work.

2. The R-process code: r-Java

Fundamental to our code is the waiting-point approximation which sets lower limits on the temperature of 2 × 109 K and neutron density of 1020 cm-3 for all our simulations (Cameron et al. 1983b). Each of the modules of r-Java; namely NSE, static and dynamic have their own appropriate set of input parameters and numerical solution methods (see Fig. 1). For example, in the dynamic case, the user can choose the amount of heating from neutrinos released through β-decay, turn on/off nuclear fission, alpha or neutron decay and specify the density profile of the expanding material. The nuclear input data that is provided along with r-Java uses nuclear masses from experimental sources as much as possible (Audi & Wapstra 1995), and where none was available we used 3. For the two nuclear inputs: the neutron separation energy (Sn) and the β-decay rates (λβ), we used 3. A freedom provided to the users of r-Java is the ability to modify any nuclear input in order to simulate the effect that such a change would have on abundance yields. There are many different outputs that can be generated by r-Java as well as the ability to animate how a given parameter or abundance changes during a simulation run. A unique module included in r-Java displays the periodic table and in real-time shows which elements are being created and destroyed during a simulation. Also included is a module to study nucleosynthesis in a Quark-Nova (Jaikumar et al. 2007), which will be discussed in detail in subsequent papers. The following subsections are dedicated to describing in detail the physics involved in each of the modules of r-Java and presenting corresponding simulation results.

2.1. Nuclear statistical equilibrium

As a prelude to more complex situations, we begin by developing a NSE code where the determining nuclear property is the binding energy. Consider an astrophysical plasma composed of photons, free neutrons and protons, and a mix of seed nuclei such as 56Fe, with temperature and density high enough that the nuclear reactions assembling nuclei from free neutrons and protons are much faster than the expansion timescale. For high enough temperatures T9 = T/(109) K  ≃ 6, the system would then reach NSE with forward and backward reaction rates equal. As long as such conditions prevail, the detailed balance equation holds (1)where μi is the chemical potential of the species i (with Zi, Ni, Ai the corresponding atomic number, neutron number, and mass number, respectively). Assuming all the nuclei in the gas follow Maxwell-Boltzmann statistics, the particle number densities are (e.g. Pathria 1977) (2)where B is the binding energy and g is the statistical weight. In addition, mass and charge conservation imply that Combining Eqs. (2) and (1) with the conservation Eqs. (3) and (4), we obtain the following equations for the neutron and proton chemical potentials: (5)For any given Ye (electron fraction), T (the temperature in Kelvin) and baryon mass density ρ (in g cm-3), we then solve this system by the Newton-Raphson method for μn and μp (Fig. 1), which we use to obtain the mass fractions of every nucleus through Eq. (8), where NA is the Avogadro number. The nuclei that dominate the abundance distribution in NSE depend only on ρ, T and Ye, once nuclear parameters B (binding energy) and statistical weights are known. High densities favour large nuclei due to the ρ-dependence while high temperatures favour light nuclei since the Planck distribution contains energetic photons which photodisintegrate heavy nuclei. We found the detailed balance method (Eq. (1)) numerically amenable to the Newton-Raphson method. However, we should caution that the success of the Newton-Raphson method depends sensitively on a proper choice for the initial value of μp, μn and the starting composition. We chose to use the NSE calculations of 5 as a test of our code. We have plotted our NSE results in Fig. 2. We also display the mass fractions as a function of Ye in Fig. 3. Comparing the former to the results of 1, and the latter with Fig. 3 of 5, we find very good agreement with these established results. In Fig. 4, we have also shown a new plot: the mass fractions obtained with our code as a function of density, for fixed Ye = 0.5 and T = 6.5 × 109 K.

thumbnail Fig. 2

Abundances of nucleons and selected nuclei in the NSE state with varying temperature, at ρ = 1 × 107 g cm-3 and Ye = 0.5. The neutron and proton abundances overlap for temperatures exceeding  ~7T9.

Open with DEXTER

thumbnail Fig. 3

Abundances of selected nuclei in the NSE state with a varying electron fraction, at ρ = 1 × 107 g cm-3 and T = 3.5T9.

Open with DEXTER

thumbnail Fig. 4

Abundances of nucleons and selected nuclei in the NSE state with varying density, at Ye = 0.5 and T = 6.5 T9.

Open with DEXTER

3. Reaction network and static case

In general, a whole variety of reactions can occur simultaneously, producing and destroying nucleus i, whose reaction induced density evolution is described by (e.g. Landau & Lifshitz 1958): (6)where one-body reactions, such as decays and photodisintegrations (e.g. rj = λjnj) and two body reactions ( where δjk is the Kronecker delta) dominate. The individual Ni s in Eq. (6) are positive or negative numbers specifying creation and annihilation. Then, the number abundances Yi and mass abundances Xi are given by In terms of Yi, the nuclear reaction network Eq. (6) becomes (Hix & Meyer 2006) (9)Applying this reaction network to the r-process, we will limit ourselves at the moment to the following interactions:

  • neutron capture;

  • photodisintegration;

  • β-decay.

For the static simulation, we assume a static, non-expanding chunk of neutron-rich material near drip density; the drip density is the density at which point it is no longer energetically favourable to add more neutrons to a nucleus. This is estimated to be  ≈ 4 × 1011 g/cc in the Baym-Pethick-Sutherland equation of state, which describes nuclear matter upto such densities (Baym et al. 1971b, herafter BPS)2. The chunk of neutron-rich material is bombarded by a large and constant flux of neutrons over a given period of time. Equation (9) is then (10)where  ⟨ σv ⟩  is the thermally averaged neutron capture cross section. Since we have a high neutron density (1022–1026 cm-3) and high temperatures (T ≳ 2 × 109 K), neutron captures and photodisintegrations occur much faster than β-decays. For such conditions we can assume that (n) ⇌ (γ,n) equilibrium is established within every isotopic chain. We also assume steady flow while the high density and temperature conditions remain. The steady flow approximation implies that the rate of inflow of a nucleus into the network equals the outflow – this simplifies the solution to the reaction network, enabling one to solve for the abundances in an algebraic manner (Cameron et al. 1983a), viz., (11)The steady flow approximation can be further simplified when the previously mentioned conditions (high neutron density and temperature) hold and (n) ⇌ (γ,n) equilibrium applies.

thumbnail Fig. 5

Final abundances in the static simulations after one second of neutron bombardment and at different physical conditions. Top left: nn = 1026 cm-3 and T9 = 2, top right: nn = 1028 cm-3 and T9 = 2, bottom left: nn = 1026 cm-3 and T9 = 4, bottom right: nn = 1028 cm-3 and T9 = 4. The solar abundance is shown in each case with vertical lines marking the location of the important peaks.

Open with DEXTER

The abundance in a particular chain is then primarily that of the nuclide that has slowest β-decay rate, this is the waiting-point nucleus. In the standard waiting-point approximation we can express the abundance ratio of two neighbouring isotopes within one isotopic chain in terms of only the neutron separation energy (Sn). In terms of the reduced reaction network (10), we then have (Qian 2003): (12)Further, we introduce the abundance of an isotopic chain, Y(Z), defined as the sum of every isotope abundance in that chain: Y(Z) =  ∑ AY(Z,A). Also, g(Z,A) in the above equation is the statistical weight of nucleus (Z,A). With Eq. (12), we can now express every Y(Z,A) in terms of only the corresponding Y(Z) and population coefficients defined through Y(Z,A) = P(Z,A)Y(Z) which depend only on nn, T and Sn. This allows us to dramatically reduce the number of equations from several thousand down to 136, the number of isotopic chains we have in our network. Now we study the evolution of Y(Z) through β-decay (Cameron et al. 1983a; Goriely & Arnould 1996): (13)This leads to (14)Introducing effective β-decay rates: , we obtain a system of coupled differential equations: (15)As the abundances Y(Z) can vary by many orders of magnitude within the integration interval, Eq. (15) is a very stiff system, so we use an implicit scheme, the Backward Euler Scheme, to treat the problem. To present the capabilities of r-Java in the static regime we varied the neutron density number and the temperature of the material to observe their impact on the r-process abundance (see Fig. 5).

Comparing the left and right figures in the top panel of Fig. 5 shows clearly that the third r-process peak at A = 195 is realized when the neutron number density is high (nn = 1028 cm-3) and the temperature is close to 2 × 109 K. In reality, in dense neutron-rich environments neutron number densities above 1032 cm-3 are possible, and we might expect large A nuclei to be copiously produced, but as we have not included fission in the static simulations, we have restrained our code at present to lower neutron densities (nn = 1026 − 1028 cm-3), so that nuclei do not build up without limit.

3.1. Static simulation results

From the top panels of Fig. 5, it can be seen that our static simulations produce a weak third r-process peak at A = 195 even when the neutron number density is relatively high (nn = 1028 cm-3). Higher neutron densities can produce a strong third peak (nn ~ 1032   cm-3), but the absence of fission presently limits our static code to moderate neutron densities. Absence of fission cycling also leads to a broad peak structure. The bottom panels of Fig. 5 show the effects of doubling the temperature, relative to the top panel, keeping neutron density the same. As temperature increases heavy nuclei tend to disappear since photodisintegrations of neutron-rich nuclei become highly effective freezing the nuclear flow at lower A. In this case, the peaks in the mass distribution appear at lower mass numbers compared to the corresponding usual r-process peaks. This was found to be the case even for longer periods of neutron exposure, possibly suggesting that photodisintegrations are dominating over neutron captures within an isotopic chain, violating the conditions of (n,γ) equilibrium. Such a “shift” of the peaks to incorrect locations in mass number has been noticed in earlier dynamical situations, such as neutron star mergers (e.g. Freiburghaus et al. 1999b), where it is attributed to the choice of Ye. In our static simulations Ye does not cause the shift. Rather, in our case, this effect is driven by the high temperature, which favours photodisintegrations and takes the r-process path closer to the valley of β-stability, thereby downplaying the role of waiting-point nuclei in determining the abundance peaks (see Goriely & Arnould 1996, for the high neutron density and low temperature, T9 < 2, regime). Extremely high neutron densities are required to restore the peak to its correct position.

thumbnail Fig. 6

Final abundances in the static simulations showing how the elements β-decay back to stability once bombardment of neutrons stops. The cases shown are one second and one hundred seconds following the end of neutron bombardment. The left panel is the corresponding Z versus N distribution while the right panel displays how the abundances versus charge number evolve in time.

Open with DEXTER

3.2. The β-decay module

thumbnail Fig. 7

Schematic representation of a single time step in our dynamic code. In this chart ρ0 is the initial crust density, Ye,0 is the initial electron fraction, Y(Z0) is the initial chemical abundance of the crust, T0 is the initial temperature, Z0 is the initial atomic number and n0 is the initial particle density of the crust. The neutron density (nn) is first found using the secant method and then reaction network is solved using the Crank-Nicholson method. After solving the reaction network evolution of the density (ρ(t)), heating due to β and neutron decay and nuclear fission is calculated. The neutron density is then re-evaluated and finally the maximal temperature change criterion is checked before moving on to the next time step.

Open with DEXTER

Although we have implemented the waiting-point approximation, r-Java also provides the option of evolving the elements by β-decay back to stability from the time bombardment of neutrons is abruptly terminated. Thus, one can isolate key β-decay rates for particular reaction networks and studies of freeze-out, which is known to be important, for example in explaining the rare-earth peak near europium (Surman et al. 1997). We use the Euler implicit scheme for each isotopic chain of elements with decay rate of each element exactly the ones we used to compute the effective rates defined in Eq. (15). Shown in Fig. 6 are final abundances from the static simulation (under same input conditions as in Fig. 5) β-decaying one second and one hundred seconds after the bombardment of neutrons has stopped. An added convenience in r-Java is the fact that the β-decay module can be run independent of nucleosynthesis calculations by specifying the initial composition as a vector. We added a modification of the bi-diagonal solver, so that this solver accepts a vector of variable size. Currently the β-decay module is limited to nuclei with mass between 50 and 300.

4. Dynamic case

The dynamic module of r-Java allows the user to input the initial temperature (T0), density (ρ0), electron fraction (Ye,0) and dominantly abundant element (Z0) of a chunk of material. The material then expands at a user specified expansion timescale (τ) for a given duration. In order for r-Java to accommodate a wide variety of ejection mechanisms we have allowed the density profile () to be freely modified by the user3. The relativistic expansion regime will be discussed in an upcoming paper, here we discuss the non-relativistic case. Users of r-Java are given the choice to include a simple fission model for dynamic simulations, in which the last species of the network can be specified. For our dynamic simulations we assume that the last species in the network is Z = 92, A = 276 which fissions into two fragments-one with Z = 40 and one with Z = 52. The code does not presently include other types of fission (neutron-induced, neutrino-induced, β-delayed) or β-delayed neutron emission. Further options given to the user of r-Java for dynamic simulations are the ability to specify the heat fraction deposited from β-decays, and allowing β-decay and neutron decay to be turned on or off.

Once the user provides the input parameters as above, the initial entropy/nucleon is determined (in kB = 1 units) from the initial density, electron fraction and temperature using (17)where μe = (3π2Ye   ρ)1/3 is the chemical potential of the fully degenerate relativistic electrons, mN = 939.1 MeV is the neutron mass and g0 = 2 (spins) is the statistical weight of the neutron.

Once the code begins to run, at each and every time step, we update the abundance of nucleus i if the fractional change dYi/Yi < 0.1, else we decrease the time step until this condition is satisfied. From Eq. (18), this implies that |ni/ni − ∂ρ/ρ| < 0.1. For non-relativistic expansion, both ni/ni and ∂ρ/ρ are small, but for relativistic expansion, ∂ρ/ρ can be large at early times. We can still ensure that the condition on the Yi’s is satisfied by choosing our time step small enough that ∂ρ/ρ < 0.1, provided ∂ρ/ρ dominates over ni/ni. This is equivalent to an upper limit on the entropy which can be converted to an upper limit on the initial temperature Tup (K) = , where n0 is the initial ejecta number density (cm-3). For example, taking Z0 = ZFe (iron-rich ejecta), we find that an initial density ρ0 ~ 1013 g/cc requires that the initial temperature satisfy T0 ≲ 1011 K. If the temperature input by the user violates this density relation r-Java asks for a new choice of initial temperature and shows the user the upper bound. Introducing heat into the system increases the entropy and can potentially reduce the maximum initial temperature for which our abundance criterion can be satisfied. However, heating due to nuclear transmutations can be ignored at early times in the case of fast relativistic expansion, since the expansion is much faster than typical β-decay timescales. Thus, the above estimate for the maximum initial temperature remains valid for fast expansions. It is automatically satisfied for the non-relativistic case, where ∂ρ/ρ is much smaller. Furthermore, once satisfied at the outset, the abundance criterion dYi/Yi < 0.1 is seen to be always satisfied for all future times. These features are naturally built into our code, as detailed by the algorithm flowchart seen in Fig. 7.

Equation (6) for the evolution of Y(Z) must now be modified due to the density evolution, which was neglected before in the static calculations. The total derivative of Y(Z) is (18)The first term, which describes changes in ejecta density, can be rewritten in term of Yi and then incorporated into our reaction network: (19)For the dynamic simulation, we used a Crank-Nicholson method to solve Eq. (19). As in the static case, we were able to solve the linear system of abundance equations to find the variations in abundances within a time step. These nuclear transmutations generate entropy and lead to heating, so to keep track of this, we inverted the entropy from these changes to compute the change in temperature. This procedure was implemented in a self-consistent manner, i.e, the new temperature was fed back to check the value of entropy. The heating of the material via β-decay is added to Eq. (17) as δs = δq/T and this whole expression is inverted to get the new temperature after each time step. We check that the temperature remains above 2 × 109 K at all times during r-process, so that the waiting-point approximation we employ is still valid.

To complete the time step, we also need to explicitly calculate the neutron density, which we can do using (n) ⇌ (γ,n) equilibrium, as long as temperatures and densities remain high enough. This ensures equilibrium within any isotopic chain, which will instantly redistribute Y(Z) between all isotopes in accordance with the neutron separation energy (Sn). Once again, we use the population coefficients which depend on the neutron density, so to get the equilibrium value of neutron density at each time step, we solve the implicit equation of mass conservation: (20)where  ∑ A is the mass of all nuclei except neutrons. This quantity depends on the neutron density through the population coefficients, which give the effective mass of every isotopic chain. We used a secant method to solve this equation and get the equilibrium value of nn. Then the system is completely described and the code can move to the next time step. We can run the code as long as the waiting-point approximation is valid, for nn > 1020 cm-3 and T > 2 × 109 K (see algorithm in Fig. 7).

4.1. Dynamic simulation results

thumbnail Fig. 8

The final abundances (solid red line) for dynamic simulations (one second duration) of nucleosynthesis from expansion of a chunk of pure iron situated in neutron-rich matter with Ye = 0.03. Solar abundances are shown as black crosses. For each panel the same initial temperature (4 × 109 K) and density (1 × 1011 g/cc) is used. The left column of panels use an expansion timescale of τ = 0.1 s and for the right column of panels the expansion timescale is τ = 0.001 s. Each row considers a different density profile: top ρ(t) ∝ (t/τ)-1, middle ρ(t) ∝ (t/τ)-2, bottom ρ(t) ∝ (t/τ)-3.

Open with DEXTER

thumbnail Fig. 9

The final Z versus N distribution is plotted for each dynamic simulation (after one second duration). In each panel the solid green line is at N = Z, and the large crosses mark the location of the closed neutron shells. For each panel the same initial temperature (4 × 109 K), density (1 × 1011 g/cc) and electron fraction (0.03) are used. The left column of panels use an expansion timescale of τ = 0.1 s and for the right column of panels, the expansion timescale is τ = 0.001 s. Each row considers a different density profile: top ρ(t) ∝ (t/τ)-1, middle ρ(t) ∝ (t/τ)-2, bottom ρ(t) ∝ (t/τ)-3.

Open with DEXTER

thumbnail Fig. 10

The evolution of temperature over the simulation run is plotted for each dynamic simulation (50% heating from β-decays is assumed). For each panel the same initial temperature (4 × 109 K), density (1 × 1011 g/cc) and electron fraction (0.03) are used. The left column of panels use an expansion timescale of τ = 0.1 s and for the right column of panels the expansion timescale is τ = 0.001 s. Each row considers a different density profile: top ρ(t) ∝ (t/τ)-1, middle ρ(t) ∝ (t/τ)-2, bottom ρ(t) ∝ (t/τ)-3.

Open with DEXTER

In Fig. 8, we explore the effect on final abundance of changing the density profile of the expanding material. To do this we considered the following general density profile: (21)For Figs. 810 in the top panels α = 1, the middle panels α = 2 and the bottom panels α = 3, where increasing values of α roughly simulate increasing degrees of freedom for expansion. The left column of these figures use an expansion timescale of τ = 0.001 s while for the right column τ = 0.1 s. We left the remaining input parameters the same for all panels (T0 = 4 × 109 K, Ye,0 = 0.03, Z0 = 26 and a simulation duration of one second). From the left column of Fig. 8 it can be seen that the density drops off too rapidly to produce elements heavier than A = 120 in all but the α = 1 case. This is because the expansion is too fast for neutron captures to occur for a sufficiently long period of time to make heavy elements. An opposite trend can be seen in the τ = 0.1 s cases where the heavy-element production increases with the degrees of freedom for expansion. This is because for slow expansions, the density, particularly for low degrees of freedom (α = 1) remains high for an extended period of time, blocking β-decays that are needed to move the flow up to higher mass numbers. Consequently, the third peak is absent or weak. For α = 2 or 3, the density fall off is rapid enough to compensate this effect, so a strong 3rd peak is seen. The effect of the β-decay blocking can be clearly seen in the Z vs. N plots shown in Fig. 9. In the top panel, right column of Fig. 9, it is clear that the slower expansion case allows elements to build up beyond the magic numbers, while for faster expansion (top-left panel of Fig. 9) or multiple degrees of freedom (right column of 2nd and 3rd panel) there exists clear steps at each magic number, revealing the imprint of the extra stability. It is interesting to note that the final abundances seen in the top left (τ = 0.001, α = 1) and bottom right (τ = 0.1, α = 3) panels of both Figs. 8 and 9 are remarkably similar. This similarity is due to the fact that after one second of expansion the effects of the different α and τ values will directly offset one another leaving the final density of the expanding chunk the same. The slight difference in final abundances of the two simulations are due to temperature effects. Finally, Fig. 10 displays the temperature evolution for each of our dynamical simulations. The left column of Fig. 10 exemplifies how in the short expansion timescale adiabatic cooling dominates the heating from β-decays, while in the longer expansion timescale cases (the right column of Fig. 10) β-decays have a strong heating effect early on, specially for α = 2, 3.

5. Discussion and conclusion

We have constructed and tested a nucleosynthesis code for the r-process under static and dynamic conditions. Called r-Java, it is made available online through an easy-to-use interface at http://quarknova.ucalgary.ca/. The user is given a number of input options as well as the ability to easily modify the nuclear data used in the code. In this paper, we have studied NSE, as well as both static and dynamic simulations.

Exploring static simulations for a range of neutron irradiation and temperatures, we found that neutron density  ~1028 cm-3 and above is required for a strong third r-process peak. With increasing temperature beyond T9 = 2, nuclei photodisintegrate extremely fast, before significant pile up can occur at the true waiting-point nuclei. This resulted in abundance peaks at smaller mass numbers. This feature is compensated by increasing the neutron density; however, the effect of fission will then play a role in determining the final abundances. We intend to include fission rates in our static simulations in the future, so that higher neutron densities can be reliably addressed. We have included a β-decay module that allows the user to track the flow of β-decays once neutron bombardment stops. This should be useful in studying freeze-out conditions and quasi-equilibrium reaction networks.

Using the dynamic module of r-Java we compared different density profiles for two expansion timescales. We found that at the shorter expansion timescale (τ = 0.001 s) only the slowest evolving density profile (ρ(t) ∝ ρ0(t/τ)-1) produced a significant abundance of elements above A = 120. For the longer expansion timescale (τ = 0.1 s) the heavy abundance yield increased with faster evolving density profiles, where the strongest A = 195 peak occurred in the simulation that considered ρ(t) ∝ ρ0(t/τ)-3. At present, all our calculations are carried out in the waiting-point approximation.

The large qualitative differences in the r-process pattern we observe as a function of dynamical expansion parameters are important for alternate sites of the r-process, such as neutron star mergers or spherically symmetric decompressing neutron star matter. In the case of mergers, tidal forces can lead to a “tube of toothpaste” quasi-1D ejection, whereas spherical decompressions resemble the 3D case. From Fig. 8, we see that fast expansion is required if heavy elements are made in mergers, while less energy (slower expansion) is required for spherical decompression. This is the case for example in Quark-Novae, which will be discussed in a follow-up paper. We also note the curious fact that the peak at A = 3, corresponding to tritium, is formed from the sequence of reactions comprised of neutron decay followed by formation of deuterium and then tritium. Since the half-life of tritium is relatively short (12.32 years), a large amount of its daughter element, He-3, would be formed and dispelled into the surrounding space. In the case of a dual-shock Quark-Nova (Leahy & Ouyed 2008; Ouyed et al. 2010), this can have important consequences, through spallation on supernova ejecta, for cosmic rays and anomalous abundance ratios in some supernovae (e.g., Cas A).

Our future efforts are aimed in three directions: (1) applications of r-Java to Neutron Star Mergers and the Quark-Nova scenario. The Quark-Nova is an ejection mechanism for neutron star material powered by the phase conversion from hadronic to deconfined quark matter occurring at the core of a neutron star (Ouyed et al. 2002; Keränen et al. 2005). The extreme energetics and high neutron-to-seed ratio makes the Quark-Nova an ideal candidate for an astrophysical r-process site (Jaikumar et al. 2007); (2) including the effects of freeze-out through a set of quasi-equilibrium reaction networks and other processes neglected in Eq. (13) such as β-delayed neutron emission, to obtain a more realistic pattern of final abundances, which can then be used for schematic chemical evolution studies; (3) based on 4, further numerical study of the dynamics of the nuclear-quark conversion front inside the neutron star to better constrain the Lorentz factor of the ejecta, and possibly develop a semi-analytic approach to constrain the parameter space, similar to what has been done for neutrino-driven winds from supernovae.


1

http://quarknova.ucalgary.ca/software/rJava/. If using this code for calculations for a paper, we request that you please acknowledge this website and cite the present paper.

2

This leaves only T, nn and Z0 (the dominant nucleus at that density) as inputs since Ye is prescribed.

3

The default density profile in the dynamic module of r-Java comes from studies of non-relativistic decompression by 2: (16)

Acknowledgments

R.O., M.K. and N.K. are supported by the Natural Sciences and Engineering Research Council of Canada. N.K. acknowledges support from Alberta Ingenuity and the Killam Trusts. P.J. acknowledges start-up funding from California State University Long Beach.

References

  1. Audi, G., & Wapstra, A. H. 1995, Nucl. Phys. A, 595, 409 [NASA ADS] [CrossRef] [Google Scholar]
  2. Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14 [NASA ADS] [CrossRef] [Google Scholar]
  4. Burbidge, G. R., Burbidge, W. A., Fowler, W. A., & Hoyle, F. 1957, Rev. Mod. Phys., 29, 547 [NASA ADS] [CrossRef] [Google Scholar]
  5. Cameron, A. G. W. 1957, Pub. Astr. Soc. Pacific, 69, 201 [NASA ADS] [CrossRef] [Google Scholar]
  6. Cameron, A. G. W., Cowan, J. J., Klapdor, H. V., et al. 1983a, Astrophys. Space Sci., 91, 221 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cameron, A. G. W., Cowan, J. J., & Truran, J. W. 1983b, Astrophys. Space Sci., 91, 235 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cardall, C. Y., & Fuller, G. M. 1997, ApJ, 486, L111 [NASA ADS] [CrossRef] [Google Scholar]
  9. Clifford, F., & Tayler, R. 1965, Mem. Roy. Astron. Soc., 69, 21 [NASA ADS] [Google Scholar]
  10. Freiburghaus, C., Rembges, F., Rauscher, T., et al. 1999a, ApJ, 516, 381 [NASA ADS] [CrossRef] [Google Scholar]
  11. Freiburghaus, C., Rosswog, S., & Thielemann, F.-K. 1999b, ApJ, 525, L121 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  12. Goriely, S., & Arnould, M. 1996, A&A, 312, 327 [NASA ADS] [Google Scholar]
  13. Goriely, S., Demetriou, P., Janka, H.-T., Pearson, J. M., & Samyn, M. 2004, Nucl. Phys. A, 758, 587 [NASA ADS] [CrossRef] [Google Scholar]
  14. Hix, W. R., & Meyer, M. S. 2006, Nucl. Phys. A, 777, 188 [NASA ADS] [CrossRef] [Google Scholar]
  15. Jaikumar, P., Meyer, B. S., Otsuki, K., & Ouyed, R. 2007, A&A, 471, 227 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Keränen, P., Ouyed, R., & Jaikumar, P. 2005, ApJ, 618, 485 [NASA ADS] [CrossRef] [Google Scholar]
  17. Landau, L. D., & Lifshitz, E. M. 1958, in Statistic. Physical (London: Pergamon Press), Ch. X, 5 [Google Scholar]
  18. Leahy, D., & Ouyed, R. 2008, MNRAS, 387, 1193 [NASA ADS] [CrossRef] [Google Scholar]
  19. Meyer, B. S. 1989, ApJ, 343, 254 [NASA ADS] [CrossRef] [Google Scholar]
  20. Meyer, B. S., & Brown, J. S. 1997, ApJS, 112, 199 [NASA ADS] [CrossRef] [Google Scholar]
  21. Möller, P., Nix, J. R., Myers, W. D., & Swiatecki, W. J. 1995, Atom. Data Nucl. Data Tables, 59, 185 [NASA ADS] [CrossRef] [Google Scholar]
  22. Möller, P., Nix, J. R., & Kratz, K.-L. 1997, Atom. Data Nucl. Data Tables, 66, 131 [NASA ADS] [CrossRef] [Google Scholar]
  23. Niebergal, B., Ouyed, R., & Jaikumar, P. 2010, Phys. Rev. C, 82, 062801 [NASA ADS] [CrossRef] [Google Scholar]
  24. Otsuki, K., Tagoshi, H., Kajino, T., & Wanajo, S. 2000, ApJ, 533, 424 [NASA ADS] [CrossRef] [Google Scholar]
  25. Ouyed, R., Dey, J., & Dey, M. 2002, A&A, 390, 39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Ouyed, R., Kostka, M., Koning, N., Leahy, D., & Steffen, W. 2010 [arXiv:1010.5530] [Google Scholar]
  27. Pathria, R. K. 1977, Stat. Mech. (Oxford: Pergamon Press) [Google Scholar]
  28. Qian, Y.-Z. 2000, ApJ, 534, L67 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  29. Qian, Y.-Z. 2003, Prog. Part. Nucl. Phys., 50, 153 [NASA ADS] [CrossRef] [Google Scholar]
  30. Qian, Y.-Z., & Woosley, S. E. 1996, ApJ, 471, 331 [NASA ADS] [CrossRef] [Google Scholar]
  31. Qian, Y.-Z., & Wasserburg, G. J. 2007, Phys. Rep., 442, 237 [NASA ADS] [CrossRef] [Google Scholar]
  32. Seitenzahl, I. R., Timmes, F. X., Mariu-Laflèche, A., et al. 2008, ApJ, 685, 129 [NASA ADS] [CrossRef] [Google Scholar]
  33. Sneden, C., Cowan, J. J., Lawler, J. E., et al. 2003, ApJ, 591, 936 [NASA ADS] [CrossRef] [Google Scholar]
  34. Sumiyoshi, K., Terasawa, M., Mathews, G. J., et al. 2001, ApJ, 562, 880 [NASA ADS] [CrossRef] [Google Scholar]
  35. Surman, R., Engel, J., Bennett, J. R., & Meyer, B. S. 1997, Phys. Rev. Lett., 79, 1809 [NASA ADS] [CrossRef] [Google Scholar]
  36. Surman, R., McLaughlin, G. C., & Hix, W. R. 2006, ApJ, 643, 1057 [NASA ADS] [CrossRef] [Google Scholar]
  37. Takahashi, K., Witti, J., & Janka, H.-T. 1994, A&A, 286, 857 [NASA ADS] [Google Scholar]
  38. Thompson, T. A., Burrows, A., & Meyer, B. S. 2001, ApJ, 562, 887 [NASA ADS] [CrossRef] [Google Scholar]
  39. Truran, J. W., Cowan, J. J., Pilachowski, C. A., & Sneden, C. 2002, PASP, 114, 1293 [NASA ADS] [CrossRef] [Google Scholar]
  40. Wanajo, S., Toshitaka, K., Mathews, G. J., & Otsuki, K. 2001, ApJ, 554, 578 [NASA ADS] [CrossRef] [Google Scholar]
  41. Wanajo, S., Tamamura, M., Naoki, I., et al. 2003, ApJ, 593, 968 [NASA ADS] [CrossRef] [Google Scholar]
  42. Woosley, S. E., Wilson, J. R., Mathews, G. J., Hoffman, R. D., & Meyer, B. S. 1994, ApJ, 433, 229 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

Schematic representation of r-Java’s capabilities and tools (see text for meaning of symbols). More details can be found in the manual available on the quark nova project website: http://quarknova.ucalgary.ca/.

Open with DEXTER
In the text
thumbnail Fig. 2

Abundances of nucleons and selected nuclei in the NSE state with varying temperature, at ρ = 1 × 107 g cm-3 and Ye = 0.5. The neutron and proton abundances overlap for temperatures exceeding  ~7T9.

Open with DEXTER
In the text
thumbnail Fig. 3

Abundances of selected nuclei in the NSE state with a varying electron fraction, at ρ = 1 × 107 g cm-3 and T = 3.5T9.

Open with DEXTER
In the text
thumbnail Fig. 4

Abundances of nucleons and selected nuclei in the NSE state with varying density, at Ye = 0.5 and T = 6.5 T9.

Open with DEXTER
In the text
thumbnail Fig. 5

Final abundances in the static simulations after one second of neutron bombardment and at different physical conditions. Top left: nn = 1026 cm-3 and T9 = 2, top right: nn = 1028 cm-3 and T9 = 2, bottom left: nn = 1026 cm-3 and T9 = 4, bottom right: nn = 1028 cm-3 and T9 = 4. The solar abundance is shown in each case with vertical lines marking the location of the important peaks.

Open with DEXTER
In the text
thumbnail Fig. 6

Final abundances in the static simulations showing how the elements β-decay back to stability once bombardment of neutrons stops. The cases shown are one second and one hundred seconds following the end of neutron bombardment. The left panel is the corresponding Z versus N distribution while the right panel displays how the abundances versus charge number evolve in time.

Open with DEXTER
In the text
thumbnail Fig. 7

Schematic representation of a single time step in our dynamic code. In this chart ρ0 is the initial crust density, Ye,0 is the initial electron fraction, Y(Z0) is the initial chemical abundance of the crust, T0 is the initial temperature, Z0 is the initial atomic number and n0 is the initial particle density of the crust. The neutron density (nn) is first found using the secant method and then reaction network is solved using the Crank-Nicholson method. After solving the reaction network evolution of the density (ρ(t)), heating due to β and neutron decay and nuclear fission is calculated. The neutron density is then re-evaluated and finally the maximal temperature change criterion is checked before moving on to the next time step.

Open with DEXTER
In the text
thumbnail Fig. 8

The final abundances (solid red line) for dynamic simulations (one second duration) of nucleosynthesis from expansion of a chunk of pure iron situated in neutron-rich matter with Ye = 0.03. Solar abundances are shown as black crosses. For each panel the same initial temperature (4 × 109 K) and density (1 × 1011 g/cc) is used. The left column of panels use an expansion timescale of τ = 0.1 s and for the right column of panels the expansion timescale is τ = 0.001 s. Each row considers a different density profile: top ρ(t) ∝ (t/τ)-1, middle ρ(t) ∝ (t/τ)-2, bottom ρ(t) ∝ (t/τ)-3.

Open with DEXTER
In the text
thumbnail Fig. 9

The final Z versus N distribution is plotted for each dynamic simulation (after one second duration). In each panel the solid green line is at N = Z, and the large crosses mark the location of the closed neutron shells. For each panel the same initial temperature (4 × 109 K), density (1 × 1011 g/cc) and electron fraction (0.03) are used. The left column of panels use an expansion timescale of τ = 0.1 s and for the right column of panels, the expansion timescale is τ = 0.001 s. Each row considers a different density profile: top ρ(t) ∝ (t/τ)-1, middle ρ(t) ∝ (t/τ)-2, bottom ρ(t) ∝ (t/τ)-3.

Open with DEXTER
In the text
thumbnail Fig. 10

The evolution of temperature over the simulation run is plotted for each dynamic simulation (50% heating from β-decays is assumed). For each panel the same initial temperature (4 × 109 K), density (1 × 1011 g/cc) and electron fraction (0.03) are used. The left column of panels use an expansion timescale of τ = 0.1 s and for the right column of panels the expansion timescale is τ = 0.001 s. Each row considers a different density profile: top ρ(t) ∝ (t/τ)-1, middle ρ(t) ∝ (t/τ)-2, bottom ρ(t) ∝ (t/τ)-3.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.