# R-Process Nucleosynthesis in MHD Jet Explosions of Core-Collapse Supernovae

Authors: Motoaki Saruwatari; Masa-aki Hashimoto; Ryohei Fukuda; Shin-ichiro Fujimoto
Journal: Journal of Astrophysics
 (2013)

---

## Abstract

We investigate ther-process nucleosynthesis during the magnetohydrodynamical (MHD) explosion of a supernova in a helium star of 3.3 M⊙, 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 r-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 r-process site owing to MHD jets in supernovae.


---

## Body

## 1. Introduction

Study of ther-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 Se80, Te130, and Pt195, in the abundance pattern for the solar system r-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 r-process [2, 3], explosion mechanism is not still completely resolved, where supernova explosions are originated from the gravitational collapse of massive stars of M≥8 M⊙ [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 M≥13 M⊙ 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 M⊙ stars to explode, the explosion energies are very weak [8–15]. Therefore, a plausible site/mechanism of the r-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 [16–19] 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 [20–22]. 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 r-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, 23–27]. 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 r-process during the supernova explosion of massive stars under the assumption of adiabatic explosion. They have shown the pattern of distributions of the r-elements of the solar system abundances. However, they have also found that the electron fractions (Ye) increase significantly enough to destroy the r-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 r-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 r-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 r-process. It is noted that Winteler et al. [32] have emphasized the possible importance of the early appearance of r-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.3M⊙ untill the final simulation time tf≃600 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 1.6×1014 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 1014–1016 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 r-process under adiabatic MHD explosion [30], we include the effects of neutrinos using a leakage scheme [35–40] with some modifications explained in Section 2. Finally, we investigate the possibility of the r-process in the MHD jets with use of our large nuclear reaction network. We find the region that produces the r-process elements having the particular distribution of low Ye.In Section2, our supernova models that include the neutrino effects are given, and we also explain the initial models and r-process networks. The results of the r-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]:
(1)dρdt+ρ∇·v=0,ρdvdt=-∇P-ρ∇Φ+14π(∇×B)×B,∂B∂t=∇×(v×B),∇2Φ=4πGρ,ρd(e/ρ)dt=-P∇·v+Q+-Q-,
where ρ is the density, v is the velocity, P is the pressure, B is the magnetic field, e is the internal energy density, and Q+ and Q- 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 (Lν) at a neutrino sphere can be estimated from the average ν-energy ϵ¯ν,esc (see (15) later) that escapes freely:
(2)Lν=∫V‍ϵ¯ν,escnντescdV,
where nν is the ν-number density and τesc is the escape time for a neutrino to reach the ν-sphere Rν, that is, obtained from the leakage scheme [38, 41] in terms of the ν-mean free path (λtot) defined by
(3)∫Rν∞1λtotdr=23.The mean free path of neutrinos is given as(4)λtot-1=ρYNmuσcoh+ρYpmuσsc,p+ρYnmuσsc,n+ρYnmuσab,
where mu is the atomic mass unit and Yn, Yp, and YN are the number fractions relative to baryons for neutrons, protons, and nuclei, respectively. The values of σsc,p, σsc,n, and σcoh 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 outsideRν. Since the minimum size of the mesh interval in our hydrodynamic calculations is 105 cm, we prescribe the time step so that for each time step of Δt, neutrinos run by cΔt, which is typically 104 cm in our simulations.The terms ofQ+ and Q- are calculated in the following [42, 43]. Outside the neutrino sphere, the two terms are, respectively,
(5)Q+=σabnnFν,Q-=Γpnp+Γnnn+Ree+Rγ,
where nn(np), Fν, and σab 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 (n+νe→p+e-), that is, the most important heating source and a function of average energy and density, respectively [44]. Neutrino production (emission) rates, Γp, Γn, and Ree, Rγ are explained below (see (8) and the equations below (18)). Inside the neutrino sphere, we calculate only the term Q- as follows:
(6)Q-=ϵ¯ν,escnντesc.For theν¯-absorption by free protons, we can get both Q+ and Q- by replacing physical quantities between ν and ν¯.


### 2.3. Neutrino Processes and Physical Inputs

The change in electron fractionYe is given by (e.g., Kotake et al. [39])
(7)dYedt=-ΓpYp-ΓNYN-cλνYn+cλν¯Yp,
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 Ye after around 200 ms measured from the bounce.The electron capture rate by a proton (p+e-→n+νe) with Qp=1.3 MeV is obtained from Epstein and Pethick [45]:
(8)Γp=12π3ℏGF2CV2(1+3a2)(h3c3)2Ip,(9)Ip=∫Qμe‍dEeEe2Eνmfe(1-fν).Herem=2, Q=Qp+μν, and fα is the Fermi-Dirac distribution of particles α(=e,ν):
(10)fα=11+e(ϵα-μα)/Tα,
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, Tν≠Tm. The fundamental constant of GF is the Fermi coupling, CV=0.97 the pseudovector coupling, and CA the axial vector coupling ones. We set the ratio |CA/CV|=1.27. The electron capture rate by a nucleus of the atomic number Z with the Q-value QN is
(11)ΓN=12712π3hGF2CA2(Z-20)(h3c3)2IN,
where IN is obtained from (9) with the values of m=2 and Q=QN+μν. Note that this capture process is assumed to be inhibited above the neutron number N=40 due to the effects of shell blocking [45–47]. 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(12)Eν,p=JpIp,Eν,N=JNIN,
where Jp and JN are obtained from Ip and IN with m=3, respectively [45, 48]. The average energy of neutrinos emitted by electron captures is written as follows:
(13)ϵ¯ν=Eν,pY˙p+Eν,NY˙NY˙p+Y˙N.This average energy is added to obtain theν-energy density at the next time step of simulations.Inside theν-sphere, Tν and μν are evaluated in terms of nν and the ν-energy density eν:
(14)nν=Tν3(ℏc)3∫‍4πϵν2fνdϵν,eν=Tν4(ℏc)3∫‍4πϵν3fνdϵν,
where the Fermi-Dirac distribution is assumed for neutrinos. Outside the ν-sphere, neutrino radiation can be approximated to be the black body one with μν=0. The average ν-energy is written as follows;
(15)ϵ¯ν,esc=F3F2Tν,sp,
where Tν,sp is the ν-temperature at the ν-sphere and the Fermi integrals, F2 and F3, are, respectively:
(16)F2=∫0∞‍x21+exdx,F3=∫0∞‍x31+exdx.We can approximate these terms by the analytic formula given by Epstein and Pethick [45]. Since Tν≥Tm inside the ν-sphere, the baryon energy density em increases due to the energy flow from neutrinos:
(17)demdt=cλνeν,
where λν is the mean free path of relevant neutrinos [44]. In this region, we replace Tν,sp by Tm to obtain ϵ¯ν,esc 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Γn can be applied to anti-neutrinos of the reaction, n+e+→p+ν¯e, with the substitutions
(18)μe⟶-μe,μν⟶μν¯,Tν⟶Tν¯.Other neutrino processes are as follows:(19)e++e-⟶νx+νx¯,γ⟶νx+νx¯,x=e,τ,μ.Here,x 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 νe or νe¯ by electron-positron pair annihilation is given by
(20)Ree(νe,νe¯)=(C1+C2)νe,νe¯36σ0c(mec2)2ϵe-·ϵe+×(1-fνe(ϵ))ee(1-fνe¯(ϵ))ee,
where σ0=1.76×10-44 cm2, and ϵe- or ϵe+ indicates electron/positron energy density. The weak interaction constants are (C1+C2)=(CV-CA)2+(CV+CA)2, and (1-fνe(ϵ))ee is the blocking factor in the neutrino phase space and approximately expressed by
(21)(1-fνe(ϵ))ee=(1+exp(-(12F4(ηe)F3(ηe)+12F4(-ηe)F3(-ηe)-ηνe)))-1,
where Fn(η) means fermi integral:
(22)Fn(η)=Tn+1(ℏc)n+1∫η∞‍xn1-ex-ηdx,
with η=μe/T.For the production ofνμ, νμ¯ and ντ, ντ¯, the corresponding rate is
(23)Ree(νx)=(C1+C2)νx,νx¯9σ0c(mec2)2×ϵe-·ϵe+(1-fνx(ϵ))ee2,
where (C1+C2)νx,νx¯=(CV-CA)2+(CV+CA-2)2.The rate of creation ofνe or νe¯ by the decay of transversal plasmons can be written with sufficient accuracy as
(24)Rγ(νe,νe¯)=π33α*CV2σ0c2(mec2)2T8(hc)6γ6exp(-γ)×(1+γ)(1-fνe(ϵ))γ(1-fνe¯(ϵ))γ,
and the corresponding rate for producing νx becomes
(25)Rγ(νx)=4π33α*(CV-1)2σ0c2(mec2)2T8(hc)6γ6×exp(-γ)(1+γ)(1-fνx)γ.The fine structure constant,α*=1/137.036, and γ=5.5565×10-2(1/3)(π2+3ηe2), and (1-fνx)γ is the blocking factor,
(26)(1-fνx)γ=[1+exp(-(1+12γ21+γ-ηνx))]-1.


### 2.4. Neutrinos outside the Neutrino Sphere

Neutrino number densitynν on each mesh can be calculated as follows:
(27)dnνdt=Γpnp+Γnnn+Ree+Rγ-nν,esc,
where nν,esc 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:
(28)nν,esc=nνβτesc,
where βτesc is the neutrino escape time scale and it is estimated as follows,
(29)τesc=max(ΔRc,3ΔR2π2cλtot),
where ΔR is the distance from each point to the neutrino sphere Rν. 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 β=23, 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:
(30)∂nν∂t+1r∂(nνrc)∂r=Γpnp+Γnnn+Ree+Rγ-nab,
where nab is the absorbed neutrino density calculated from the heating rate (Q+nν/eν). 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 Figure1, 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.Table 1
Initial parameters of precollapse models.
Model
T
/
|
W
| (%)
E
m
/
|
W
| (%)
X
0
*
Z
0
*
Ω
0 (s−1)
B
0 (G)
Model 1
0
0
0
0
0
0
Model 2
0.5
0.1
1
1
5.2
5.4
×
1
0
12
Model 3
0.5
0.1
0.5
1
7.9
1.0
×
1
0
13
Model 4
0.5
0.1
0.1
1
42.9
5.2
×
1
0
13
Model 5
1.5
0.1
0.1
1
72.9
5.2
×
1
0
13
Note:X0*=X0/108 cm and Z0*=Z0/108 cm. Models 1 to 4 have the same initial parameters as those adopted by Nishimura et al. [30]. We add a model 5 that has the largest value of a parameter of Ω0 among five models, where other parameters are the same as those of the model 4.(a) It shows the evolution of the neutrino luminosity after the collapse for a model 1. (b) It shows changes in electron fraction (Ye) when the central density becomes 1012, 1013, and 1014 g cm−3, respectively, and neutrino fraction (Yν) when the central density becomes 1013 and 1014 g cm−3, respectively. Our results could approximate those of Liebendörfer et al. [31].
(a)
(b)


### 2.5. Initial Models andr-Process Networks

In all computations, spherical coordinates(r,θ) are adopted. The computational region is set to be 0≤r≤4000 km and 0≤θ≤π/2, where the included mass in the precollapse models amounts to 1.42 M⊙. The first quadrant of the meridian section is covered with 300(r)×30(θ) 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 of0.449≤Ye≤0.49 from the center to the place of 1.3 M⊙ (r=2200 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 Bϕ as follows:
(31)Ω(X,Z)=Ω0X02X2+X02Z04Z4+Z04,Bϕ(X,Z)=B0X02X2+X02Z04Z4+Z04,
where X and Z are the distances from the rotational axis and the equatorial plane. X0 and Z0 are model parameters. Both Ω0 and B0 are the initial values at X=0 and Z=0. 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 Ω0 among five models. As shown in the next section, among the five models, model 5 can eject material with very low electron fraction Ye.Here, we note that Winteler et al. [32] have used a shellular rotation law as their initial model of 15 M⊙ [51]. The value of T/|W|=7.63×10-3 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 Em/|W|=2.63×10-8. 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 ther-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 r-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 Z=100. Included reactions are two-body ones, (n,γ), (p,γ), (α,γ), (p,n), (α,p), (α,n), plus three-body one, (3α,γ), 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.


## 2.1. MHD Equations

Ideal MHD equations are enumerated as follows [13, 39]:
(1)dρdt+ρ∇·v=0,ρdvdt=-∇P-ρ∇Φ+14π(∇×B)×B,∂B∂t=∇×(v×B),∇2Φ=4πGρ,ρd(e/ρ)dt=-P∇·v+Q+-Q-,
where ρ is the density, v is the velocity, P is the pressure, B is the magnetic field, e is the internal energy density, and Q+ and Q- 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 (Lν) at a neutrino sphere can be estimated from the average ν-energy ϵ¯ν,esc (see (15) later) that escapes freely:
(2)Lν=∫V‍ϵ¯ν,escnντescdV,
where nν is the ν-number density and τesc is the escape time for a neutrino to reach the ν-sphere Rν, that is, obtained from the leakage scheme [38, 41] in terms of the ν-mean free path (λtot) defined by
(3)∫Rν∞1λtotdr=23.The mean free path of neutrinos is given as(4)λtot-1=ρYNmuσcoh+ρYpmuσsc,p+ρYnmuσsc,n+ρYnmuσab,
where mu is the atomic mass unit and Yn, Yp, and YN are the number fractions relative to baryons for neutrons, protons, and nuclei, respectively. The values of σsc,p, σsc,n, and σcoh 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 outsideRν. Since the minimum size of the mesh interval in our hydrodynamic calculations is 105 cm, we prescribe the time step so that for each time step of Δt, neutrinos run by cΔt, which is typically 104 cm in our simulations.The terms ofQ+ and Q- are calculated in the following [42, 43]. Outside the neutrino sphere, the two terms are, respectively,
(5)Q+=σabnnFν,Q-=Γpnp+Γnnn+Ree+Rγ,
where nn(np), Fν, and σab 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 (n+νe→p+e-), that is, the most important heating source and a function of average energy and density, respectively [44]. Neutrino production (emission) rates, Γp, Γn, and Ree, Rγ are explained below (see (8) and the equations below (18)). Inside the neutrino sphere, we calculate only the term Q- as follows:
(6)Q-=ϵ¯ν,escnντesc.For theν¯-absorption by free protons, we can get both Q+ and Q- by replacing physical quantities between ν and ν¯.


## 2.3. Neutrino Processes and Physical Inputs

The change in electron fractionYe is given by (e.g., Kotake et al. [39])
(7)dYedt=-ΓpYp-ΓNYN-cλνYn+cλν¯Yp,
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 Ye after around 200 ms measured from the bounce.The electron capture rate by a proton (p+e-→n+νe) with Qp=1.3 MeV is obtained from Epstein and Pethick [45]:
(8)Γp=12π3ℏGF2CV2(1+3a2)(h3c3)2Ip,(9)Ip=∫Qμe‍dEeEe2Eνmfe(1-fν).Herem=2, Q=Qp+μν, and fα is the Fermi-Dirac distribution of particles α(=e,ν):
(10)fα=11+e(ϵα-μα)/Tα,
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, Tν≠Tm. The fundamental constant of GF is the Fermi coupling, CV=0.97 the pseudovector coupling, and CA the axial vector coupling ones. We set the ratio |CA/CV|=1.27. The electron capture rate by a nucleus of the atomic number Z with the Q-value QN is
(11)ΓN=12712π3hGF2CA2(Z-20)(h3c3)2IN,
where IN is obtained from (9) with the values of m=2 and Q=QN+μν. Note that this capture process is assumed to be inhibited above the neutron number N=40 due to the effects of shell blocking [45–47]. 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(12)Eν,p=JpIp,Eν,N=JNIN,
where Jp and JN are obtained from Ip and IN with m=3, respectively [45, 48]. The average energy of neutrinos emitted by electron captures is written as follows:
(13)ϵ¯ν=Eν,pY˙p+Eν,NY˙NY˙p+Y˙N.This average energy is added to obtain theν-energy density at the next time step of simulations.Inside theν-sphere, Tν and μν are evaluated in terms of nν and the ν-energy density eν:
(14)nν=Tν3(ℏc)3∫‍4πϵν2fνdϵν,eν=Tν4(ℏc)3∫‍4πϵν3fνdϵν,
where the Fermi-Dirac distribution is assumed for neutrinos. Outside the ν-sphere, neutrino radiation can be approximated to be the black body one with μν=0. The average ν-energy is written as follows;
(15)ϵ¯ν,esc=F3F2Tν,sp,
where Tν,sp is the ν-temperature at the ν-sphere and the Fermi integrals, F2 and F3, are, respectively:
(16)F2=∫0∞‍x21+exdx,F3=∫0∞‍x31+exdx.We can approximate these terms by the analytic formula given by Epstein and Pethick [45]. Since Tν≥Tm inside the ν-sphere, the baryon energy density em increases due to the energy flow from neutrinos:
(17)demdt=cλνeν,
where λν is the mean free path of relevant neutrinos [44]. In this region, we replace Tν,sp by Tm to obtain ϵ¯ν,esc 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Γn can be applied to anti-neutrinos of the reaction, n+e+→p+ν¯e, with the substitutions
(18)μe⟶-μe,μν⟶μν¯,Tν⟶Tν¯.Other neutrino processes are as follows:(19)e++e-⟶νx+νx¯,γ⟶νx+νx¯,x=e,τ,μ.Here,x 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 νe or νe¯ by electron-positron pair annihilation is given by
(20)Ree(νe,νe¯)=(C1+C2)νe,νe¯36σ0c(mec2)2ϵe-·ϵe+×(1-fνe(ϵ))ee(1-fνe¯(ϵ))ee,
where σ0=1.76×10-44 cm2, and ϵe- or ϵe+ indicates electron/positron energy density. The weak interaction constants are (C1+C2)=(CV-CA)2+(CV+CA)2, and (1-fνe(ϵ))ee is the blocking factor in the neutrino phase space and approximately expressed by
(21)(1-fνe(ϵ))ee=(1+exp(-(12F4(ηe)F3(ηe)+12F4(-ηe)F3(-ηe)-ηνe)))-1,
where Fn(η) means fermi integral:
(22)Fn(η)=Tn+1(ℏc)n+1∫η∞‍xn1-ex-ηdx,
with η=μe/T.For the production ofνμ, νμ¯ and ντ, ντ¯, the corresponding rate is
(23)Ree(νx)=(C1+C2)νx,νx¯9σ0c(mec2)2×ϵe-·ϵe+(1-fνx(ϵ))ee2,
where (C1+C2)νx,νx¯=(CV-CA)2+(CV+CA-2)2.The rate of creation ofνe or νe¯ by the decay of transversal plasmons can be written with sufficient accuracy as
(24)Rγ(νe,νe¯)=π33α*CV2σ0c2(mec2)2T8(hc)6γ6exp(-γ)×(1+γ)(1-fνe(ϵ))γ(1-fνe¯(ϵ))γ,
and the corresponding rate for producing νx becomes
(25)Rγ(νx)=4π33α*(CV-1)2σ0c2(mec2)2T8(hc)6γ6×exp(-γ)(1+γ)(1-fνx)γ.The fine structure constant,α*=1/137.036, and γ=5.5565×10-2(1/3)(π2+3ηe2), and (1-fνx)γ is the blocking factor,
(26)(1-fνx)γ=[1+exp(-(1+12γ21+γ-ηνx))]-1.


## 2.4. Neutrinos outside the Neutrino Sphere

Neutrino number densitynν on each mesh can be calculated as follows:
(27)dnνdt=Γpnp+Γnnn+Ree+Rγ-nν,esc,
where nν,esc 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:
(28)nν,esc=nνβτesc,
where βτesc is the neutrino escape time scale and it is estimated as follows,
(29)τesc=max(ΔRc,3ΔR2π2cλtot),
where ΔR is the distance from each point to the neutrino sphere Rν. 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 β=23, 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:
(30)∂nν∂t+1r∂(nνrc)∂r=Γpnp+Γnnn+Ree+Rγ-nab,
where nab is the absorbed neutrino density calculated from the heating rate (Q+nν/eν). 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 Figure1, 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.Table 1
Initial parameters of precollapse models.
Model
T
/
|
W
| (%)
E
m
/
|
W
| (%)
X
0
*
Z
0
*
Ω
0 (s−1)
B
0 (G)
Model 1
0
0
0
0
0
0
Model 2
0.5
0.1
1
1
5.2
5.4
×
1
0
12
Model 3
0.5
0.1
0.5
1
7.9
1.0
×
1
0
13
Model 4
0.5
0.1
0.1
1
42.9
5.2
×
1
0
13
Model 5
1.5
0.1
0.1
1
72.9
5.2
×
1
0
13
Note:X0*=X0/108 cm and Z0*=Z0/108 cm. Models 1 to 4 have the same initial parameters as those adopted by Nishimura et al. [30]. We add a model 5 that has the largest value of a parameter of Ω0 among five models, where other parameters are the same as those of the model 4.(a) It shows the evolution of the neutrino luminosity after the collapse for a model 1. (b) It shows changes in electron fraction (Ye) when the central density becomes 1012, 1013, and 1014 g cm−3, respectively, and neutrino fraction (Yν) when the central density becomes 1013 and 1014 g cm−3, respectively. Our results could approximate those of Liebendörfer et al. [31].
(a)
(b)


## 2.5. Initial Models andr-Process Networks

In all computations, spherical coordinates(r,θ) are adopted. The computational region is set to be 0≤r≤4000 km and 0≤θ≤π/2, where the included mass in the precollapse models amounts to 1.42 M⊙. The first quadrant of the meridian section is covered with 300(r)×30(θ) 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 of0.449≤Ye≤0.49 from the center to the place of 1.3 M⊙ (r=2200 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 Bϕ as follows:
(31)Ω(X,Z)=Ω0X02X2+X02Z04Z4+Z04,Bϕ(X,Z)=B0X02X2+X02Z04Z4+Z04,
where X and Z are the distances from the rotational axis and the equatorial plane. X0 and Z0 are model parameters. Both Ω0 and B0 are the initial values at X=0 and Z=0. 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 Ω0 among five models. As shown in the next section, among the five models, model 5 can eject material with very low electron fraction Ye.Here, we note that Winteler et al. [32] have used a shellular rotation law as their initial model of 15 M⊙ [51]. The value of T/|W|=7.63×10-3 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 Em/|W|=2.63×10-8. 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 ther-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 r-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 Z=100. Included reactions are two-body ones, (n,γ), (p,γ), (α,γ), (p,n), (α,p), (α,n), plus three-body one, (3α,γ), 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, andr-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 Section2. Our results of MHD calculations are summarized in Table 2, where Eexp is the explosion energy when the shock reaches the edge of the Fe-core and Mej 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 40° 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 1010 K around the original layers of the Si + Fe core, where the nuclear statistical equilibrium is realized.Table 2
Calculated quantities that are crucial in ther-process.
Model
t
b
t
f
T
/
|
W
|
f
E
m
/
|
W
|
f
E
exp
*
M
ej
/
M
⊙
M
r
ej
/
M
⊙
Model 1
111
283
0
0
0.023
—
—
Model 2
125
311
6.91
0.053
0.127
—
—
Model 3
129
329
8.74
0.116
0.164
—
—
Model 4
133
433
8.80
0.142
1.13
0.111
—
Model 5
180
624
15.3
0.339
0.484
0.022
5.90
×
1
0
-
3
Note:tb indicates the time (ms) at the bounce. The calculations are stopped at the time tf (ms). The ratios T/|W|f and Em/|W|f are expressed in %. Eexp*=Eexp/1051 ergs. Mej is the sum of the ejected tracer particles. Mrej is the ejected mass of the r-element for A≥63.Figure 2
Snapshots of tracer particles att=100 ms (top) and t=200 ms (bottom) after the bounce in models 4 (left) and 5 (right).In model 4, the equatorial region is ejected as shown in Figure2 having rather high value of Ye≃0.50 (Figure 3(a)). In model 5, materials are ejected with the jets along the polar regions, whose total angle is subtended over 20° from the axis (see Figure 2). The corresponding evolutions of Ye relevant to the r-process are shown in Figure 4. The lowest value of Ye<0.20 is found around the polar region as seen in Figure 3(b).Contour ofYe over the range 0.1–0.5 at the final stage of calculation in model 4 (a) and model 5 (b).
(a)
(b)Time evolution ofYe in model 4 (a) and model 5 (b). For model 4, the explosion resembles to spherical one and all Ye have increased above 0.4 due to the neutrino capture process.
(a)
(b)Figure5 shows the ejected mass against Ye in the range 0.05≤Ye≤0.50. In model 4, the ejecta with Ye>0.40 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, Ye 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 Ye, n+νe→p+e- becomes ineffective.Figure 5
Ejected mass as a function ofYe at around tf=600 ms (model 5). The thick lines indicate the results of model 4 (tf=433 ms). Two separate distributions in Ye are due to different neutrino irradiation just outside the neutrino sphere.Density contour over 107–1014 g cm−3 at t=200 ms after the bounce in model 4 (a) and model 5 (b). Model 5 indicates a jet in which density is rather low along the polar axis.
(a)
(b)We note that we distribute 9000 tracer particles on each mesh point. To check the change in distribution ofYe owing to that of tracer particles for model 5, (1) we scatter 15000 tracer particles inside the computational region and (2) nine thousand particles between 500 km and 2200 km. As a result, we have confirmed that deficiency of Ye between 0.2–0.3 is the same for all cases. For case (1), tiny amounts of Ye appears below Ye=0.2. Therefore, nucleosynthesis results qualitatively do not depend on the method how to distribute tracer particles.We calculate ther-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 1010 K to around 9×109 K. Then, the nuclear reaction network of the r-process has been operated till the temperature decreases to 2-3×109 K (t~600 ms) using the results of the MHD calculations. After that, network calculations are performed until T~107 K (t~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 (n+νe→p+e-, p+ν¯e→n+e+) 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 (t~tf) under the assumption of constant neutrino luminosity, since Lν does not change appreciably after t~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 (Figure6(a)). It is, however, impossible to reproduce the r-elements even up to the second peak of the solar r-process abundance pattern, because Ye of the ejected materials distributes in the range of high values of Ye≥0.4. 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 r-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.Nucleosynthesis calculation results using ETFSI (a) and FRDM (b). The third peak is reproduced for both mass formula. However, the slope of the abundance pattern of FRDM is too steep aroundA≃200 to compare with the solar pattern. The last peak of A>200 cannot be reproduced in both cases.
(a)
(b)


## 4. Summary and Discussion

We have barely shown a qualitative possibility of ther-process nucleosynthesis during MHD explosion in a massive star of 13 M⊙. 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 Ye 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 r-process in the previous study [30], we can show at least the elemental distributions of the r-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 Ye 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 r-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 r-process for low metallicities would give a clue about nuclear cosmochronology represented by 232Th which half life is as long as the age of the universe [52].We could conclude that supernova explosions of massive stars associated with ther-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 andr-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 r-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 [11–13]. 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 Ye 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 Ye matter from deep inside the region. In this point, MRE explosion with a realistic magnetic field configuration should be investigated [20–22] to see the effects on the r-process.


### 5.5. Simulations of Other Massive Star Models

We have adopted a 3.3M⊙ 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 (1014–1016 G). This might suggest the formation of magnetar [63–72]. 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.


## 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 [11–13]. 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 Ye 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 Ye matter from deep inside the region. In this point, MRE explosion with a realistic magnetic field configuration should be investigated [20–22] to see the effects on the r-process.


## 5.5. Simulations of Other Massive Star Models

We have adopted a 3.3M⊙ 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 (1014–1016 G). This might suggest the formation of magnetar [63–72]. 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.



---

*Source: 506146-2013-10-22.xml*