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/00046361/201117008  
Published online  16 June 2011 
rJava: an rprocess code and graphical user interface for heavyelement nucleosynthesis
^{1}
Department of Physics and AstronomyUniversity of Calgary, 2500 University Drive NW, Calgary, AB T2N 1N4, Canada
email: rouyed@ucalgary.ca
^{2}
Départment de Physique, École Normale Supérieure de Cachan, 94230 Cachan, France
^{3}
Department of Physics and Astronomy, California State University Long Beach, 1250 Bellflower Blvd., Long Beach, California 90840, USA
^{4}
Institute of Mathematical Sciences, CIT Campus, Chennai 600113, India
Received: 3 April 2011
Accepted: 19 May 2011
We present rJava, an rprocess code for open use that performs rprocess nucleosynthesis calculations. Equipped with a simple graphical user interface, rJava is capable of carrying out nuclear statistical equilibrium (NSE), as well as static and dynamic rprocess calculations, for a wide range of input parameters. In this introductory paper, we present the motivation and details behind rJava and results from our static and dynamic simulations. Static simulations are explored for a range of neutron irradiation and temperatures. Dynamic simulations are studied with a parameterized expansion formula. Our code generates the resulting abundance pattern based on a general entropy expression that can be applied to both degenerate and nondegenerate matter, allowing us to track the rapid density and temperature evolution of the ejecta during the initial stages of ejecta expansion. At present, our calculations are limited to the waitingpoint approximation. We encourage the nuclear astrophysics community to provide feedback on the code and related documentation, which is available for download from the website of the QuarkNova Project: http://quarknova.ucalgary.ca/.
Key words: nuclear reactions, nucleosynthesis, abundances
© ESO, 2011
1. Introduction
The rapid neutron capture process (rprocess) 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 neutronrich astrophysical environments present ideal conditions for the rprocess to take place. Possible candidate sites discussed in the literature include the neutrinodriven, neutronrich wind from protoneutron 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 gammaray bursts (Surman et al. 2006) etc. Most importantly, abundance data on rprocess elements in metalpoor stars (Sneden et al. 2003) and certain radionuclides in meteorites (Qian & Wasserburg 2007) point toward the distinct possibility of multiple rprocess sites, so that more than one of these sites may contribute to the observed abundance pattern of the rprocess elements (Truran et al. 2002).
In 1985, an important work by Bethe & Wilson brought neutrinodriven winds from type II supernovae (SNe) to the forefront of the discussion about astrophysical sites for the rprocess. 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 selfconsistent astrophysical model that reproduces the observed abundance of rprocess elements and that at the same time is consistent with spectroscopic and meteoritic data. The most natural explanation is that the observed rprocess pattern follows from a superposition of neutron capture events with differing neutrontoseed 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 neutrontoseed ratio than type II SNe, but occur so infrequently that they fail to explain the early enrichment of rprocess elements relative to iron observed in metalpoor stars (Qian 2000). Clearly, further study is needed to better understand the rprocess conditions, so to this end we present a fairly general and flexible rprocess code that is transparent and freely available to the nuclear physics and astrophysics communities.
Fig. 1
Schematic representation of rJava’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 rprocess code came in part from recognizing the difficulty in having access to an openly shared crossplatform 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 rareisotope experiments such as FRIB later in this decade will make available new and improved nuclear data for reaction rates relevant to the rprocess. Besides this advantage, the rprocess is of wide interest to the astrophysics community, who might wish to explore different dynamical environments for the rprocess 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 neutronrich 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 realtime interpretation and analysis of simulations through data and graphs. We have made the rprocess application freely available in Javabased format on the quark nova project website^{1}, 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 rJava, 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 rJava and display results for NSE that are in agreement with previous benchmarks. In Sects. 3 and 4, we present the simplified reaction network (waitingpoint approximation) and simulation results for heavyelement 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 Rprocess code: rJava
Fundamental to our code is the waitingpoint approximation which sets lower limits on the temperature of 2 × 10^{9} K and neutron density of 10^{20} cm^{3} for all our simulations (Cameron et al. 1983b). Each of the modules of rJava; 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 rJava 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 (S_{n}) and the βdecay rates (λ^{β}), we used 3. A freedom provided to the users of rJava 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 rJava as well as the ability to animate how a given parameter or abundance changes during a simulation run. A unique module included in rJava displays the periodic table and in realtime shows which elements are being created and destroyed during a simulation. Also included is a module to study nucleosynthesis in a QuarkNova (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 rJava 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 ^{56}Fe, 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 T_{9} = T/(10^{9}) 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 Z_{i}, N_{i}, A_{i} the corresponding atomic number, neutron number, and mass number, respectively). Assuming all the nuclei in the gas follow MaxwellBoltzmann 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 Y_{e} (electron fraction), T (the temperature in Kelvin) and baryon mass density ρ (in g cm^{3}), we then solve this system by the NewtonRaphson method for μ_{n} and μ_{p} (Fig. 1), which we use to obtain the mass fractions of every nucleus through Eq. (8), where N_{A} is the Avogadro number. The nuclei that dominate the abundance distribution in NSE depend only on ρ, T and Y_{e}, 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 NewtonRaphson method. However, we should caution that the success of the NewtonRaphson 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 Y_{e} 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 Y_{e} = 0.5 and T = 6.5 × 10^{9} K.
Fig. 2
Abundances of nucleons and selected nuclei in the NSE state with varying temperature, at ρ = 1 × 10^{7} g cm^{3} and Y_{e} = 0.5. The neutron and proton abundances overlap for temperatures exceeding ~7T_{9}. 

Open with DEXTER 
Fig. 3
Abundances of selected nuclei in the NSE state with a varying electron fraction, at ρ = 1 × 10^{7} g cm^{3} and T = 3.5T_{9}. 

Open with DEXTER 
Fig. 4
Abundances of nucleons and selected nuclei in the NSE state with varying density, at Y_{e} = 0.5 and T = 6.5 T_{9}. 

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 onebody reactions, such as decays and photodisintegrations (e.g. r_{j} = λ_{j}n_{j}) and two body reactions ( where δ_{jk} is the Kronecker delta) dominate. The individual N^{i} s in Eq. (6) are positive or negative numbers specifying creation and annihilation. Then, the number abundances Y_{i} and mass abundances X_{i} are given by In terms of Y_{i}, the nuclear reaction network Eq. (6) becomes (Hix & Meyer 2006) (9)Applying this reaction network to the rprocess, we will limit ourselves at the moment to the following interactions:

neutron capture;

photodisintegration;

βdecay.
For the static simulation, we assume a static, nonexpanding chunk of neutronrich 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 × 10^{11} g/cc in the BaymPethickSutherland equation of state, which describes nuclear matter upto such densities (Baym et al. 1971b, herafter BPS)^{2}. The chunk of neutronrich 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 (10^{22}–10^{26} cm^{3}) and high temperatures (T ≳ 2 × 10^{9} 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.
Fig. 5
Final abundances in the static simulations after one second of neutron bombardment and at different physical conditions. Top left: n_{n} = 10^{26} cm^{3} and T_{9} = 2, top right: n_{n} = 10^{28} cm^{3} and T_{9} = 2, bottom left: n_{n} = 10^{26} cm^{3} and T_{9} = 4, bottom right: n_{n} = 10^{28} cm^{3} and T_{9} = 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 waitingpoint nucleus. In the standard waitingpoint approximation we can express the abundance ratio of two neighbouring isotopes within one isotopic chain in terms of only the neutron separation energy (S_{n}). 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) = ∑ _{A}Y(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 n_{n}, T and S_{n}. 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 rJava in the static regime we varied the neutron density number and the temperature of the material to observe their impact on the rprocess abundance (see Fig. 5).
Comparing the left and right figures in the top panel of Fig. 5 shows clearly that the third rprocess peak at A = 195 is realized when the neutron number density is high (n_{n} = 10^{28} cm^{3}) and the temperature is close to 2 × 10^{9} K. In reality, in dense neutronrich environments neutron number densities above 10^{32} 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 (n_{n} = 10^{26} − 10^{28} 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 rprocess peak at A = 195 even when the neutron number density is relatively high (n_{n} = 10^{28} cm^{3}). Higher neutron densities can produce a strong third peak (n_{n} ~ 10^{32} 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 neutronrich 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 rprocess 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 Y_{e}. In our static simulations Y_{e} does not cause the shift. Rather, in our case, this effect is driven by the high temperature, which favours photodisintegrations and takes the rprocess path closer to the valley of βstability, thereby downplaying the role of waitingpoint nuclei in determining the abundance peaks (see Goriely & Arnould 1996, for the high neutron density and low temperature, T_{9} < 2, regime). Extremely high neutron densities are required to restore the peak to its correct position.
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
Fig. 7
Schematic representation of a single time step in our dynamic code. In this chart ρ_{0} is the initial crust density, Y_{e,0} is the initial electron fraction, Y(Z_{0}) is the initial chemical abundance of the crust, T_{0} is the initial temperature, Z_{0} is the initial atomic number and n_{0} is the initial particle density of the crust. The neutron density (n_{n}) is first found using the secant method and then reaction network is solved using the CrankNicholson 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 reevaluated 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 waitingpoint approximation, rJava 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 freezeout, which is known to be important, for example in explaining the rareearth 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 rJava 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 bidiagonal 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 rJava allows the user to input the initial temperature (T_{0}), density (ρ_{0}), electron fraction (Y_{e,0}) and dominantly abundant element (Z_{0}) of a chunk of material. The material then expands at a user specified expansion timescale (τ) for a given duration. In order for rJava to accommodate a wide variety of ejection mechanisms we have allowed the density profile () to be freely modified by the user^{3}. The relativistic expansion regime will be discussed in an upcoming paper, here we discuss the nonrelativistic case. Users of rJava 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 fragmentsone with Z = 40 and one with Z = 52. The code does not presently include other types of fission (neutroninduced, neutrinoinduced, βdelayed) or βdelayed neutron emission. Further options given to the user of rJava 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 k_{B} = 1 units) from the initial density, electron fraction and temperature using (17)where μ_{e} = (3π^{2}Y_{e} ρ)^{1/3} is the chemical potential of the fully degenerate relativistic electrons, m_{N} = 939.1 MeV is the neutron mass and g_{0} = 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 dY_{i}/Y_{i} < 0.1, else we decrease the time step until this condition is satisfied. From Eq. (18), this implies that ∂n_{i}/n_{i} − ∂ρ/ρ < 0.1. For nonrelativistic expansion, both ∂n_{i}/n_{i} and ∂ρ/ρ are small, but for relativistic expansion, ∂ρ/ρ can be large at early times. We can still ensure that the condition on the Y_{i}’s is satisfied by choosing our time step small enough that ∂ρ/ρ < 0.1, provided ∂ρ/ρ dominates over ∂n_{i}/n_{i}. This is equivalent to an upper limit on the entropy which can be converted to an upper limit on the initial temperature T_{up} (K) = , where n_{0} is the initial ejecta number density (cm^{3}). For example, taking Z_{0} = Z_{Fe} (ironrich ejecta), we find that an initial density ρ_{0} ~ 10^{13} g/cc requires that the initial temperature satisfy T_{0} ≲ 10^{11} K. If the temperature input by the user violates this density relation rJava 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 nonrelativistic case, where ∂ρ/ρ is much smaller. Furthermore, once satisfied at the outset, the abundance criterion dY_{i}/Y_{i} < 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 Y_{i} and then incorporated into our reaction network: (19)For the dynamic simulation, we used a CrankNicholson 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 selfconsistent 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 × 10^{9} K at all times during rprocess, so that the waitingpoint 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 (S_{n}). 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 n_{n}. 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 waitingpoint approximation is valid, for n_{n} > 10^{20} cm^{3} and T > 2 × 10^{9} K (see algorithm in Fig. 7).
4.1. Dynamic simulation results
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 neutronrich matter with Y_{e} = 0.03. Solar abundances are shown as black crosses. For each panel the same initial temperature (4 × 10^{9} K) and density (1 × 10^{11} 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 
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 × 10^{9} K), density (1 × 10^{11} 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 
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 × 10^{9} K), density (1 × 10^{11} 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. 8–10 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 (T_{0} = 4 × 10^{9} K, Y_{e,0} = 0.03, Z_{0} = 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 heavyelement 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 (topleft 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 rprocess under static and dynamic conditions. Called rJava, it is made available online through an easytouse 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 ~10^{28} cm^{3} and above is required for a strong third rprocess peak. With increasing temperature beyond T_{9} = 2, nuclei photodisintegrate extremely fast, before significant pile up can occur at the true waitingpoint 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 freezeout conditions and quasiequilibrium reaction networks.
Using the dynamic module of rJava 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 waitingpoint approximation.
The large qualitative differences in the rprocess pattern we observe as a function of dynamical expansion parameters are important for alternate sites of the rprocess, 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” quasi1D 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 QuarkNovae, which will be discussed in a followup 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 halflife of tritium is relatively short (12.32 years), a large amount of its daughter element, He3, would be formed and dispelled into the surrounding space. In the case of a dualshock QuarkNova (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 rJava to Neutron Star Mergers and the QuarkNova scenario. The QuarkNova 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 neutrontoseed ratio makes the QuarkNova an ideal candidate for an astrophysical rprocess site (Jaikumar et al. 2007); (2) including the effects of freezeout through a set of quasiequilibrium 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 nuclearquark conversion front inside the neutron star to better constrain the Lorentz factor of the ejecta, and possibly develop a semianalytic approach to constrain the parameter space, similar to what has been done for neutrinodriven winds from supernovae.
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.
The default density profile in the dynamic module of rJava comes from studies of nonrelativistic 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 startup funding from California State University Long Beach.
References
 Audi, G., & Wapstra, A. H. 1995, Nucl. Phys. A, 595, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Burbidge, G. R., Burbidge, W. A., Fowler, W. A., & Hoyle, F. 1957, Rev. Mod. Phys., 29, 547 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, A. G. W. 1957, Pub. Astr. Soc. Pacific, 69, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, A. G. W., Cowan, J. J., Klapdor, H. V., et al. 1983a, Astrophys. Space Sci., 91, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, A. G. W., Cowan, J. J., & Truran, J. W. 1983b, Astrophys. Space Sci., 91, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Cardall, C. Y., & Fuller, G. M. 1997, ApJ, 486, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Clifford, F., & Tayler, R. 1965, Mem. Roy. Astron. Soc., 69, 21 [NASA ADS] [Google Scholar]
 Freiburghaus, C., Rembges, F., Rauscher, T., et al. 1999a, ApJ, 516, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Freiburghaus, C., Rosswog, S., & Thielemann, F.K. 1999b, ApJ, 525, L121 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Goriely, S., & Arnould, M. 1996, A&A, 312, 327 [NASA ADS] [Google Scholar]
 Goriely, S., Demetriou, P., Janka, H.T., Pearson, J. M., & Samyn, M. 2004, Nucl. Phys. A, 758, 587 [NASA ADS] [CrossRef] [Google Scholar]
 Hix, W. R., & Meyer, M. S. 2006, Nucl. Phys. A, 777, 188 [NASA ADS] [CrossRef] [Google Scholar]
 Jaikumar, P., Meyer, B. S., Otsuki, K., & Ouyed, R. 2007, A&A, 471, 227 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keränen, P., Ouyed, R., & Jaikumar, P. 2005, ApJ, 618, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1958, in Statistic. Physical (London: Pergamon Press), Ch. X, 5 [Google Scholar]
 Leahy, D., & Ouyed, R. 2008, MNRAS, 387, 1193 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, B. S. 1989, ApJ, 343, 254 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, B. S., & Brown, J. S. 1997, ApJS, 112, 199 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Möller, P., Nix, J. R., & Kratz, K.L. 1997, Atom. Data Nucl. Data Tables, 66, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Niebergal, B., Ouyed, R., & Jaikumar, P. 2010, Phys. Rev. C, 82, 062801 [NASA ADS] [CrossRef] [Google Scholar]
 Otsuki, K., Tagoshi, H., Kajino, T., & Wanajo, S. 2000, ApJ, 533, 424 [NASA ADS] [CrossRef] [Google Scholar]
 Ouyed, R., Dey, J., & Dey, M. 2002, A&A, 390, 39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ouyed, R., Kostka, M., Koning, N., Leahy, D., & Steffen, W. 2010 [arXiv:1010.5530] [Google Scholar]
 Pathria, R. K. 1977, Stat. Mech. (Oxford: Pergamon Press) [Google Scholar]
 Qian, Y.Z. 2000, ApJ, 534, L67 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Qian, Y.Z. 2003, Prog. Part. Nucl. Phys., 50, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, Y.Z., & Woosley, S. E. 1996, ApJ, 471, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, Y.Z., & Wasserburg, G. J. 2007, Phys. Rep., 442, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Seitenzahl, I. R., Timmes, F. X., MariuLaflèche, A., et al. 2008, ApJ, 685, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Sneden, C., Cowan, J. J., Lawler, J. E., et al. 2003, ApJ, 591, 936 [NASA ADS] [CrossRef] [Google Scholar]
 Sumiyoshi, K., Terasawa, M., Mathews, G. J., et al. 2001, ApJ, 562, 880 [NASA ADS] [CrossRef] [Google Scholar]
 Surman, R., Engel, J., Bennett, J. R., & Meyer, B. S. 1997, Phys. Rev. Lett., 79, 1809 [NASA ADS] [CrossRef] [Google Scholar]
 Surman, R., McLaughlin, G. C., & Hix, W. R. 2006, ApJ, 643, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, K., Witti, J., & Janka, H.T. 1994, A&A, 286, 857 [NASA ADS] [Google Scholar]
 Thompson, T. A., Burrows, A., & Meyer, B. S. 2001, ApJ, 562, 887 [NASA ADS] [CrossRef] [Google Scholar]
 Truran, J. W., Cowan, J. J., Pilachowski, C. A., & Sneden, C. 2002, PASP, 114, 1293 [NASA ADS] [CrossRef] [Google Scholar]
 Wanajo, S., Toshitaka, K., Mathews, G. J., & Otsuki, K. 2001, ApJ, 554, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Wanajo, S., Tamamura, M., Naoki, I., et al. 2003, ApJ, 593, 968 [NASA ADS] [CrossRef] [Google Scholar]
 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
Fig. 1
Schematic representation of rJava’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 
Fig. 2
Abundances of nucleons and selected nuclei in the NSE state with varying temperature, at ρ = 1 × 10^{7} g cm^{3} and Y_{e} = 0.5. The neutron and proton abundances overlap for temperatures exceeding ~7T_{9}. 

Open with DEXTER  
In the text 
Fig. 3
Abundances of selected nuclei in the NSE state with a varying electron fraction, at ρ = 1 × 10^{7} g cm^{3} and T = 3.5T_{9}. 

Open with DEXTER  
In the text 
Fig. 4
Abundances of nucleons and selected nuclei in the NSE state with varying density, at Y_{e} = 0.5 and T = 6.5 T_{9}. 

Open with DEXTER  
In the text 
Fig. 5
Final abundances in the static simulations after one second of neutron bombardment and at different physical conditions. Top left: n_{n} = 10^{26} cm^{3} and T_{9} = 2, top right: n_{n} = 10^{28} cm^{3} and T_{9} = 2, bottom left: n_{n} = 10^{26} cm^{3} and T_{9} = 4, bottom right: n_{n} = 10^{28} cm^{3} and T_{9} = 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 
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 
Fig. 7
Schematic representation of a single time step in our dynamic code. In this chart ρ_{0} is the initial crust density, Y_{e,0} is the initial electron fraction, Y(Z_{0}) is the initial chemical abundance of the crust, T_{0} is the initial temperature, Z_{0} is the initial atomic number and n_{0} is the initial particle density of the crust. The neutron density (n_{n}) is first found using the secant method and then reaction network is solved using the CrankNicholson 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 reevaluated and finally the maximal temperature change criterion is checked before moving on to the next time step. 

Open with DEXTER  
In the text 
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 neutronrich matter with Y_{e} = 0.03. Solar abundances are shown as black crosses. For each panel the same initial temperature (4 × 10^{9} K) and density (1 × 10^{11} 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 
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 × 10^{9} K), density (1 × 10^{11} 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 
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 × 10^{9} K), density (1 × 10^{11} 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.