We investigate the -process nucleosynthesis during the magnetohydrodynamical (MHD) explosion of a supernova in a helium star of 3.3 , where effects of neutrinos are taken into account using the leakage scheme in the two-dimensional (2D) hydrodynamic code. Jet-like explosion due to the combined effects of differential rotation and magnetic field is able to erode the lower electron fraction matter from the inner layers. We find that the ejected material of low electron fraction responsible for the -process comes out from just outside the neutrino sphere deep inside the Fe-core. It is found that heavy element nucleosynthesis depends on the initial conditions of rotational and magnetic fields. In particular, the third peak of the distribution is significantly overproduced relative to the solar system abundances, which would indicate a possible -process site owing to MHD jets in supernovae.

1. Introduction

Study of the -process has been developed considerably keeping pace with the terrestrial experiments of nuclear physics far from the stability line of nuclides [1]. In particular, among the three peaks, which correspond to the elements of , , and , in the abundance pattern for the solar system -elements, the transition from the second to third peak elements has been stressed by nuclear physicists [2]. Although supernovae could be one of the astrophysical sites of the -process [2, 3], explosion mechanism is not still completely resolved, where supernova explosions are originated from the gravitational collapse of massive stars of [4, 5]. However it is unclear whether neutron-rich elements could be ejected or not during the shock wave propagation. As far as the one-dimensional calculations, almost all realistic numerical simulations concerning the collapse-driven supernovae of have failed to explode the outer layer above the Fe-core due to drooping of the energetic shock wave propagation [6, 7]. Although there exist calculations for 8 and 11  stars to explode, the explosion energies are very weak [815]. Therefore, a plausible site/mechanism of the -process has not yet been clarified.

On the other hand, models of magnetorotational explosion (MRE) for core-collapse supernovae have been presented as a supernova mechanism [1619] since both rapid rotations and/or strong magnetic fields could be resulted for neutron stars after the explosions. Furthermore, MRE with a realistic magnetic field configuration has been investigated [2022]. In their series of papers, it has been shown that magnetorotational instability plays a critical role concerning the explosion energies which would explain the explosion energies of Type II and Ib supernovae. However, changes in the electron fractions and/or the heavy element nucleosynthesis have not been discussed well. Therefore, it should be studied whether MHD explosions affect the -process even within a qualitative method.

Two-dimensional (2D) magnetohydrodynamical (MHD) calculations have been performed with the various initial parameters concerning rotation and magnetic field [13, 2327]. The ZEUS-2D code [28] has been modified to include an equation of state [29], electron captures with a simple scheme of neutrino (-) transport [23]. Adopting these achievements, Nishimura et al. [30] have performed 2D/MHD calculations to study possibilities of the -process during the supernova explosion of massive stars under the assumption of adiabatic explosion. They have shown the pattern of distributions of the -elements of the solar system abundances. However, they have also found that the electron fractions () increase significantly enough to destroy the -process elemental distributions if the neutrino capture processes are included, where the processes were obtained from the results of spherical explosion calculations with the realistic neutrino transport included [31]. The problem is remained whether the adopted method for neutrino captures can be legitimate or not; effects of neutrino transport have not been included at all, and instead Nishimura et al. [30] have used the profiles of densities and temperatures obtained from the adiabatic calculations. It should be done to check their results by including the neutrino transport and to study the explosion energies even under extreme parameters of rotation/magnetic fields. In the meanwhile, Winteler et al. [32] have shown the -process nucleosynthesis with the use of results of a magnetorotationally driven supernova simulation, where they have performed 3D calculations with a 3D spectral scheme for neutrinos. Although they have obtained enough -elements, their MHD simulation has finished at around 31 ms after the bounce. It would be useful to show results of two-dimensional MHD calculations with longer simulation time and/or higher resolution for calculations of the -process. It is noted that Winteler et al. [32] have emphasized the possible importance of the early appearance of -process matter in low metallicities which could be originated by MHD explosion of supernovae.

In the present paper, we give the calculational results of the MHD explosion for the He-star of 3.3  untill the final simulation time  ms. Although these results have already been briefly reported [33], we show the details of simulation procedure and their results to compensate the contents. For the MHD calculations, five models are adopted for the initial configuration of rotation and magnetic fields. In the initial magnetic field, we assume very strong magnetic fields which validity is not known. In the observations of magnetar with the magnetic field  G [34], they have suggested a rather massive progenitor mass from the age of all the early type stars. In the complete online magnetar catalog cited by Mori et al. [34], magnetic field of 26 magnetars ranges  G which have been obtained from the analysis of the decrease in the rotational period under the assumption of magnetic dipole braking in a vacuum. These observations encourage us to set the initial condition of a strong magnetic field. Contrary to the previous investigation of the -process under adiabatic MHD explosion [30], we include the effects of neutrinos using a leakage scheme [3540] with some modifications explained in Section 2. Finally, we investigate the possibility of the -process in the MHD jets with use of our large nuclear reaction network. We find the region that produces the -process elements having the particular distribution of low .

In Section 2, our supernova models that include the neutrino effects are given, and we also explain the initial models and -process networks. The results of the -process nucleosynthesis calculations are presented in Section 3. We summarize our results in Section 4, discuss remained problems, and propose future works in Section 5.

2. Supernova Models

2.1. MHD Equations

Ideal MHD equations are enumerated as follows [13, 39]: where is the density, is the velocity, is the pressure, is the magnetic field, is the internal energy density, and and are the neutrino () heating and cooling rates, respectively. The gravitational potential is solved from the Poisson solver [28].

2.2. Neutrino Leakage Scheme

Neutrino luminosity () at a neutrino sphere can be estimated from the average -energy (see (15) later) that escapes freely: where is the -number density and is the escape time for a neutrino to reach the -sphere , that is, obtained from the leakage scheme [38, 41] in terms of the -mean free path () defined by

The mean free path of neutrinos is given as where is the atomic mass unit and , , and are the number fractions relative to baryons for neutrons, protons, and nuclei, respectively. The values of , , and are the cross sections for scattering on protons, neutrons, and nuclei, respectively.

Contrary to the original leakage scheme, we do not adopt free stream approximation for neutrinos outside . Since the minimum size of the mesh interval in our hydrodynamic calculations is  cm, we prescribe the time step so that for each time step of , neutrinos run by , which is typically  cm in our simulations.

The terms of and are calculated in the following [42, 43]. Outside the neutrino sphere, the two terms are, respectively, where , , and are the number density of free neutrons (protons), energy flux of neutrinos at each point calculated by the equation of continuity, and the -absorption cross section by free neutrons (), that is, the most important heating source and a function of average energy and density, respectively [44]. Neutrino production (emission) rates, , , and , are explained below (see (8) and the equations below (18)). Inside the neutrino sphere, we calculate only the term as follows:

For the -absorption by free protons, we can get both and by replacing physical quantities between and .

2.3. Neutrino Processes and Physical Inputs

The change in electron fraction is given by (e.g., Kotake et al. [39]) where is the mean free path of electron neutrino and is that of antineutrino [44]. The last two terms in the right hand side, which manifest the effects of -radiation, play an essential role to change after around 200 ms measured from the bounce.

The electron capture rate by a proton () with  MeV is obtained from Epstein and Pethick [45]:

Here , , and is the Fermi-Dirac distribution of particles : where and in units of the Boltzmann constant are the single particle energy and the chemical potential, respectively. We note that outside the equilibrium region between neutrinos and baryons; significant thermal deviation comes out between temperatures of neutrinos and baryons, that is, . The fundamental constant of is the Fermi coupling, the pseudovector coupling, and the axial vector coupling ones. We set the ratio . The electron capture rate by a nucleus of the atomic number with the -value is where is obtained from (9) with the values of and . Note that this capture process is assumed to be inhibited above the neutron number due to the effects of shell blocking [4547]. We should note that the temperature effects and correlations are responsible for removing the inhibition [47].

Energies of emitted neutrinos by the individual electron captures are given by where and are obtained from and with , respectively [45, 48]. The average energy of neutrinos emitted by electron captures is written as follows:

This average energy is added to obtain the -energy density at the next time step of simulations.

Inside the -sphere, and are evaluated in terms of and the -energy density : where the Fermi-Dirac distribution is assumed for neutrinos. Outside the -sphere, neutrino radiation can be approximated to be the black body one with . The average -energy is written as follows; where is the -temperature at the -sphere and the Fermi integrals, and , are, respectively:

We can approximate these terms by the analytic formula given by Epstein and Pethick [45]. Since inside the -sphere, the baryon energy density increases due to the energy flow from neutrinos: where is the mean free path of relevant neutrinos [44]. In this region, we replace by to obtain in (15), because neutrino temperature decreases due to the diffusive effect around the sphere, which would underestimate the -energy. This is our first modification of the original leakage scheme.

The similar procedure for the positron capture rate can be applied to anti-neutrinos of the reaction, , with the substitutions

Other neutrino processes are as follows:

Here, means electron, , and -neutrinos. These processes are important for a late stage of the explosion, because the most neutrino luminosity for the late stage of the explosion comes from these processes. The rates of these processes have been obtained by Ruffert et al. [36, 37]. While electron neutrinos emitted by these processes contribute to neutrino cooling and heating, , -neutrinos do only neutrino cooling. The emission rate of or by electron-positron pair annihilation is given by where  cm2, and or indicates electron/positron energy density. The weak interaction constants are , and is the blocking factor in the neutrino phase space and approximately expressed by where means fermi integral: with .

For the production of , and , , the corresponding rate is where .

The rate of creation of or by the decay of transversal plasmons can be written with sufficient accuracy as and the corresponding rate for producing becomes

The fine structure constant, , and , and is the blocking factor,

2.4. Neutrinos outside the Neutrino Sphere

Neutrino number density on each mesh can be calculated as follows: where indicates number density of protons, neutrons, and escaping neutrino density estimated by escaping time scale (see (28)), respectively. For the second modification of the leakage scheme, the last term in the right hand side in (27) is calculated in the neutrino sphere as follows: where is the neutrino escape time scale and it is estimated as follows, where is the distance from each point to the neutrino sphere . The factor is introduced and it could take the value in the range 1–5 [49, 50]. In our case, we set the value to be , considering the isotropic diffusion of neutrinos from around the neutrino sphere. The escaping neutrino density is added to the neutrino density at the neutrino sphere. Outside the neutrino sphere, we consider streaming neutrinos (see the below equation of continuity) except for slightly absorbed neutrinos. Therefore, we calculate equation of continuity: where is the absorbed neutrino density calculated from the heating rate (). We should note that Kotake et al. [23, 39] set the neutrino fraction to be zero outside the -sphere. We solved the continuity equation outside the -sphere, which affects the location of the sphere. Furthermore, they utilized a postprocessing approach for the heating term. The heating term may change the dynamics. Due to the heating, jet formation may become preferable.

In Figure 1, we show the neutrino luminosities and electron fractions calculated by the method explained in this section which corresponds to a model of spherically symmetry (model 1 in Table 1). They are compared with the figures (Figures 8 and 1(d)) in Liebendörfer et al. [31]. Considering the simple scheme of the neutrino transport, our results could approximate more accurate ones which adopt the detailed neutrino transport scheme.

2.5. Initial Models and -Process Networks

In all computations, spherical coordinates are adopted. The computational region is set to be  km and , where the included mass in the precollapse models amounts to . The first quadrant of the meridian section is covered with mesh points (Fe-core plus some amounts of Silicon-rich layer).

To acquire information of mass elements, nine thousand tracer particles are distributed on each mesh point within the region of from the center to the place of ( km) at the beginning of the collapse. Five initial models are prepared as shown in Table 1. We adopt cylindrical properties of the angular velocity and the toroidal component of the magnetic field as follows: where and are the distances from the rotational axis and the equatorial plane. and are model parameters. Both and are the initial values at and . Here, we add a model 5 in addition to the initial models adopted by Nishimura et al. [30]. This model has the largest value of among five models. As shown in the next section, among the five models, model 5 can eject material with very low electron fraction .

Here, we note that Winteler et al. [32] have used a shellular rotation law as their initial model of 15  [51]. The value of may not be so different from those of our models. On the other hand, their initial distribution of the magnetic field is assumed to be purely poloidal field having . The difference of the distribution of the magnetic field could affect the formation of jet along the polar direction.

Although significant improvements could have been performed for the nuclear data of nucleosynthesis calculations, we utilize the same nuclear reaction networks that have been constructed for the -process calculations by Nishimura et al. [30]. This is because our purpose is to investigate some qualitative effects by including the neutrino transport on the -process. Let us briefly explain the nuclear data included in the networks. The networks have been extended toward the neutron-rich side till the neutron-drip line. Each network consists of about 4000 nuclear species up to . Included reactions are two-body ones, (), (), (), (), (), (), plus three-body one, (), and their inverses. Two kinds of the network, called FRDM and ETFSI, have been constructed. The mass formula of FRDM is constructed by the Nilsson-Struntinsky model considering effects of shell and microscopic part. ETFSI approach is a semiclassical approximation to the Hartree-Fock method in which the shell corrections are calculated with the integral version of the Strutinsky theorem. Reaction rates are constructed based on experimental data if available which are supplemented by theoretical data with inverse reaction rates and partition functions with use of FRDM or ETFSI.

3. Explosion Models, Distribution of Electron Fraction, and -Process Calculations

We investigate hydrodynamical stages of the collapse, bounce, and propagation of the shock wave with use of ZEUS-2D code using a simple neutrino transport scheme as shown in Section 2. Our results of MHD calculations are summarized in Table 2, where is the explosion energy when the shock reaches the edge of the Fe-core and is the mass summed over the ejected tracer particles. We note that the explosion does not occur for model 1 to model 3. While the jet-like explosion occurs along the equator (up to from the equator) in model 4, a collimated jet emerged from the rotational axis in model 5 (Figure 2). A protoneutron star remains after the jet-like explosion. During the explosion, temperature exceeds  K around the original layers of the Si + Fe core, where the nuclear statistical equilibrium is realized.

In model 4, the equatorial region is ejected as shown in Figure 2 having rather high value of (Figure 3(a)). In model 5, materials are ejected with the jets along the polar regions, whose total angle is subtended over from the axis (see Figure 2). The corresponding evolutions of relevant to the -process are shown in Figure 4. The lowest value of is found around the polar region as seen in Figure 3(b).

Figure 5 shows the ejected mass against in the range . In model 4, the ejecta with comes from the Si-rich layers along the equatorial region, which is attributed to the enhanced centrifugal force relative to the magnetic one. We recognize that as against the spherical explosion, decreases significantly for model 5, due to the collimated jet along the rotational axis. This is because neutrino luminosity is low by a factor of ten compared to that of spherical explosion shown in Figure 1(a). This can be seen in Figure 6, where the density along polar axis is low compared to the case of model 4. Therefore, the reaction responsible for the increase in , becomes ineffective.

We note that we distribute 9000 tracer particles on each mesh point. To check the change in distribution of owing to that of tracer particles for model 5, we scatter 15000 tracer particles inside the computational region and nine thousand particles between 500 km and 2200 km. As a result, we have confirmed that deficiency of between 0.2–0.3 is the same for all cases. For case , tiny amounts of appears below . Therefore, nucleosynthesis results qualitatively do not depend on the method how to distribute tracer particles.

We calculate the -process nucleosynthesis for the explosion model 5. Before the nucleosynthesis calculation, we have assumed abundances to be in nuclear statistical equilibrium state (NSE) as has been done [44, 52]. The NSE code is used just after the temperature drops  K to around  K. Then, the nuclear reaction network of the -process has been operated till the temperature decreases to  K (~600 ms) using the results of the MHD calculations. After that, network calculations are performed until ~107 K (~10 s) with the method in Flower and Hoyle [53]. We include neutrino captures in our nucleosynthesis network that has been applied for tracer particles. The capture rates (, ) that were not incorporated by Nishimura et al. [30] are obtained from neutrino luminosity, neutrino flux, and other hydrodynamical information as we have described in the previous sections. Concerning the nucleosynthesis, we can get actually this information after the end of NSE stage (~) under the assumption of constant neutrino luminosity, since does not change appreciably after ~200 ms. We note that a difference between the neutron and proton single-particle energies in a dense medium may change neutrino capture cross sections significantly [54]. As a consequence, it is shown that the luminosities of all neutrino flavors are reduced while the spectral differences between electron neutrinos and antineutrinos are increased [55]. These changes in weak interaction processes should be examined by including them for the hydrodynamical simulations.

In case of the rapid rotation and strong magnetic field (model 4), barely jet-like explosion is obtained in the direction of the equatorial region (Figure 6(a)). It is, however, impossible to reproduce the -elements even up to the second peak of the solar -process abundance pattern, because of the ejected materials distributes in the range of high values of . In the model 5, where we have adopted a special initial configuration of concentrated magnetic field with strong differential rotation, jet-like explosion emerges in the direction of the rotational axis (Figure 6(b)). The difference is that model 5 has larger value of the angular velocity compared to model 4 by a factor of 1.7. We compared the produced heavy elements with the solar -elements in Figure 7, where results for two different mass formulas are shown. Model 5 may present a site for reproducing only the third peak with the first and second peaks underproduced.

4. Summary and Discussion

We have barely shown a qualitative possibility of the -process nucleosynthesis during MHD explosion in a massive star of 13 . Initial models have been constructed changing the distributions of rotation and magnetic fields parametrically [30].

We include neutrino effects by using the leakage scheme. This scheme treats neutrino effects approximately, where we assume Fermi-Dirac distribution for neutrinos. Furthermore, we add modification about the physical process just outside the neutrino sphere (15) and timescale of neutrino drift (28) and (29) in addition to the original leakage scheme. Validity for 2D calculations cannot be assured for this scheme, because flow to the direction is not included. If jets are so strong that the corresponding density becomes low, increase in would be suppressed compared to the 1D case. Therefore, we can say that the discussion of Nishimura et al. [30] with use of an analytical formula of 1D results is inaccurate.

Other schemes have been developed for solving neutrino transport [56]. For example, IDSA (Isotropic Diffusion Source Approximation) solves Boltzmann equation approximately, and the result obtained by using this scheme is consistent with one dimensional simulations, where the Fermi-Dirac distribution has been found to be a good approximation [56].

We have calculated MHD simulations by varying two parameters in initial models (the distribution of rotation and magnetic field). The previous adiabatic simulations without neutrino effects have succeeded in explosion except for the model 3 [30]. However, in the present case, only two models with neutrino effects whose distribution of magnetic field and rotation concentrated in the central region result in explosion. Generally, for strong and concentrated rotation case with some magnetic field, a collimated jet-like explosion occurs. For relatively weak and concentrated rotation model, large region around from rotational axis to equatorial axis is blasted, which is assisted with the magnetic field effect. In the present study, ejected region is different from that of Nishimura et al. [30]. Model 5 shows that deep region is ejected in comparison with Nishimura’s model 4. It means that the produced composition depends crucially on rotation parameter. Model 5 may present an appropriate site for reproducing only the third peak relative to the first and second peaks insufficiently built up. Contrary to the negative conclusion against the possible -process in the previous study [30], we can show at least the elemental distributions of the -elements as far as the third peak of the solar pattern is concerned. This is due to the lower neutrino luminosities from 200 to 500 ms after the bounce. At the same time, we have shown the possibility for the lower materials to be ejected significantly if the neutrino transport mechanism works appropriately.

Observations of -ray bursts associated with supernovae are rare in the present observations [57]. Considering the relations between -ray bursts and the MHD explosions, the new site of the -process to produce significantly only the third peak of the solar abundance pattern should be also rare. The nuclear process to produce the abundances after the third peak would have some relations to our MHD jet model. Since -ray bursts should have continued after the formation of the first star, a new model beyond our jet model (e.g., Winteler et al. [32]) and their motivation of the -process for low metallicities would give a clue about nuclear cosmochronology represented by Th which half life is as long as the age of the universe [52].

We could conclude that supernova explosions of massive stars associated with the -process cannot be excluded under some assumed conditions; a progenitor has special distributions of rotational/magnetic fields inside the stellar core; simple neutrino transport scheme such as a leakage scheme can be applied.

5. Future Work

We propose that both the supernova mechanism and -process nucleosynthesis still remain to be some crucial problems. This might be one of the reason why whole consistency for the origin of elements has not been understood. We should include the following effects for the -process calculations.

5.1. Detailed Neutrino Transport Scheme

In the present study, we use a simple leakage scheme for description of neutrino transport. Leakage scheme can describe neutrino transport in very easy way compared to the Boltzmann equation solver. Neutrino effect, however, is very sensitive to dynamics [58], and our method may not be inadequate for the 2D calculations. We need to adopt a more detailed neutrino transport scheme such as that of IDSA which may be capable of applying to multi-dimensional simulations. Due to the difficulty of multi-dimensional simulations with Boltzmann equation solver, IDSA simulation will play a more important role [13, 56].

5.2. Neutrino Oscillation

We do not consider neutrino oscillation because transport with neutrino oscillation is difficult to handle. Since mean energies of and neutrinos are higher than those of electron neutrino, the effect may cause significant effects for neutrino heating [59]. Furthermore, it has been shown that the charged current weak interaction processes affect the luminosities of neutrino flavors which are related to neutrino oscillations [55]. These fundamental processes should be studied with hydrodynamical simulations.

5.3. SASI

Standing accretion shock instability (SASI) has been focused on hydrodynamic simulations [1113]. It is considered to help the shock strength increase in heating regions. While SASI is in a good sense for explosion mechanism, it may play minus role for heavy element nucleosynthesis. Although a strong convection may occur behind the shock front by SASI, the convection effect tends to average the distribution [60].

5.4. Distribution of Magnetic Field

We consider only toroidal magnetic field. If we input polar magnetic field as an initial parameter [61], even a weak rotation model may lead to succeed in explosion and draw up low matter from deep inside the region. In this point, MRE explosion with a realistic magnetic field configuration should be investigated [2022] to see the effects on the -process.

5.5. Simulations of Other Massive Star Models

We have adopted a 3.3  He-core model, because it has been convenient to set up an initial model. We need to simulate for higher stellar masses (e.g., Ono et al. [62]). Neutrino luminosity may rise for more massive models. Increasing neutrino luminosity could unfortunately cause neutrino capture process furthermore and heavier nuclei may not be synthesized.

5.6. Magnetar

We find that the produced protoneutron star has very strong magnetic field ( G). This might suggest the formation of magnetar [6372]. However, in the present calculation with a numerical scheme adopted, the MRI cannot be resolved [21]; it is hopeful to study whether there exists any relation between the formation of magnetar and heavy elemental synthesis.


This work has been supported in part by a Grant-in-Aid for Scientific Research (19104006, 21540272, 24540278, and 25400281) of the Ministry of Education, Culture, Sports, Science, and Technology of Japan. The authors would like to thank the reviewers for improving the draft.