Magnetorotational Collapse of Population III Stars
arXiv:0704.1945v1  [astro-ph]  16 Apr 2007
Magnetorotational Collapse of Population III Stars
Yudai Suwa1, Tomoya Takiwaki1, Kei Kotake2, and, Katsuhiko Sato1,3
1Department of Physics, School of Science, The University of Tokyo, Tokyo 113-0033
2National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588,Japan
3Research Center for the Early Universe, School of Science, the University of Tokyo,7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-0033, Japan
suwa@utap.phys.s.u-tokyo.ac.jp
(Received 2006 October 12; accepted 2007 April 11)
Abstract
We perform a series of two-dimensional magnetorotational core-collapse simula-
tions of Pop III stars. Changing the initial distributions of rotation and magnetic
fields prior to collapse in a parametric manner, we compute 19 models. By so do-
ing, we systematically investigate how rotation and magnetic fields affect the collapse
dynamics and explore how the properties of the black-hole formations and neutrino
emissions could be affected. As for the microphysics, we employ a realistic equation of
state and approximate the neutrino transfer by a multiflavour leakage scheme. With
these computations, we find that the jet-like explosions are obtained by the magneto-
driven shock waves if the initial magnetic field is as large as 1012G. We point out that
the black-hole masses at the formation decrease with the initial field strength, on the
other hand, increase with the initial rotation rates. As for the neutrino properties,
we point out that the degree of the differential rotation plays an important role to
determine which species of the neutrino luminosity is more dominant than the others.
Furthermore, we find that the stronger magnetic fields make the peak neutrino lumi-
nosities smaller, because the magnetic pressure acts to halt the collapse in the central
regions, leading to the suppression of the releasable gravitational binding energies.
Key words: stars: supernovae: general — black hole physics — neutrinos —
methods: numerical — magnetohydrodynamics: MHD
1. INTRODUCTION
Great attention has been paid to Population III, the first stars to form in the universe,
because they are related to many unsettled problems in cosmology and the stellar physics. The
Population III (Pop III) stars ionize and enrich the metallicity of the intergalactic medium and
thus provide important clues to the subsequent star formation history (for reviews, see e.g.,
http://arxiv.org/abs/0704.1945v1
Barkana & Loeb 2001; Bromm & Larson 2004; Glover 2005). Pop III stars are also important
for the understanding of the chemical evolution history. Recent discovery of hyper metal poor
stars such as HE 0107-5240 (Christlieb et al. 2002) and HE 1327-2326 (Frebel et al. 2005)
has given us good opportunities to investigate the nucleosynthesis in Pop III stars (Heger &
Woosley 2002; Umeda & Nomoto 2002; Umeda & Nomoto 2003; Daigne et al. 2004; Iwamoto
et al. 2005). Gamma-ray bursts at very high redshift are pointed out to be accompanied by the
gravitational collapse of Pop III stars (Schneider et al. 2002; Bromm & Loeb 2006). The Swift
satellite, which is now running 1, is expected to directly detect the Pop III stars accompanied
by the high-z gamma-ray bursts.
The evolutions of Pop III stars have been also studied for long. From their studies, Pop
III stars are predicted to have been predominantly very massive with M >∼ 100M⊙ (Nakamura
& Umemura 2001; Abel et al. 2002; Bromm et al. 2002, see references therein). Massive stars in
the range of 100M⊙ <∼M
∼ 260M⊙ encounter the electron-positron pair instability during their
evolution. This instability sets off explosive oxygen burning, and if the burning provides enough
energy to reverse the collapse, the stars are thought to become pair instability supernovae (Bond
et al. 1984; Fryer et al. 2001; Heger & Woosley 2002), whose detectability has been recently
reported (Scannapieco et al. 2005; Weinmann & Lilly 2005). More massive stars, which also
encounter pair-instability, are so tightly bound and the fusion of oxygen is unable to reverse
infall. Such stars are thought to collapse to black holes (BHs) finally (Bond et al. 1984; Fryer
et al. 2001), which we pay attention to in this paper.
So far there have been a few hydrodynamic simulations studying the gravitational-
collapse of the BH forming Pop III stars. In the two-dimensional, gray neutrino transport
simulations by Fryer et al. (2001), they investigated the collapse of a rotating Pop III star of
300M⊙, leading to the BH formation. They discussed the effects of rotation on the emitted
neutrino luminosities, gravitational waves, and furthermore, the possibility of such stars to be
the gamma-ray bursts. In their Newtonian study, the central BH was excised and treated as
an absorbing boundary after the formation. Although such simplification is not easy to be
validated, they followed the dynamics long after the formation of the BH and obtained many
findings. More recently, Nakazato et al. (2006) performed one-dimensional, but, general rel-
ativistic simulations in the range of 100 ∼ 10000M⊙, in which the state-of-the-art neutrino
physics are taken into account. Their detailed calculations revealed the properties of the emer-
gent neutrino spectrum, and based on that, they discussed the detectability of such neutrinos
as the supernova relic neutrino background (see also Ando & Sato 2004; Iocco et al. 2005).
They successfully saw the formation of the apparent horizon, however, the dynamics in the
later phases was not referred.
In this paper we study the magnetorotational collapse of Pop III stars by performing
the two-dimensional magnetohydrodynamic (MHD) simulations (see, also, Akiyama et al. 2003;
1 See http://swift.gsfc.nasa.gov
Kotake et al. 2004a; Kotake et al. 2004b; Takiwaki et al. 2004; Yamada & Sawai 2004; Ardeljan
et al. 2005; Sawai et al. 2005; Obergaulinger et al. 2006, for MHD computations of core-collapse
supernovae, and Kotake et al. 2006 for a review). As for the microphysics, we employ a realistic
equation of state based on the relativistic mean field theory and take into account the neutrino
cooling by a multiflavor leakage scheme, in which state-of-the-art reactions of neutrinos are
included. In our Newtonian simulations, the formation of the BHs is ascribed to a certain
condition, and after the formation, the central region is excised and treated as an absorbing
boundary in order to follow the dynamics later on. Since the distributions of rotation and
magnetic fields in the progenitors of Pop III stars are highly uncertain, we change them in
a parametric manner and systematically investigate how rotation and magnetic fields affect
the dynamics. We also explore how the natures of explosions, the properties of the BHs and
neutrino luminosities could be affected due to the incursion of the rotation and magnetic fields.
This paper is organized as follows. In §2, we describe the numerical methods and the
initial conditions. In §3, we present the results. We give a summary and discussion in §4.
2. METHOD
2.1. Basic Equations
The basic evolution equations are written as follows,
+ ρ∇ ·v = 0, (1)
=−∇P − ρ∇Φ+
(∇×B)×B, (2)
=−P∇ ·v−Lν , (3)
=∇× (v×B) , (4)
△Φ= 4πGρ, (5)
where ρ,P,v,e,Φ,B,Lν ,
, are the mass density, the gas pressure including the radiation pres-
sure from neutrino’s, the fluid velocity, the internal energy density, the gravitational potential,
the magnetic field, the neutrino cooling rate, and Lagrange derivative, respectively. In our 2D
calculations, axial symmetry and reflection symmetry across the equatorial plane are assumed.
Spherical coordinates (r, θ) are employed with logarithmic zoning in the radial direction and
regular zoning in θ. One quadrant of the meridian section is covered with 300 (r)× 30 (θ) mesh
points. The minimum and maximum mesh spacings are 2 km and 60 km, respectively. We also
calculated some models with 60 angular mesh points, however, any significant difference was
obtained. Therefore, we will report in the following the results obtained from the models with
30 angular mesh points. We employed the ZEUS-2D code (Stone & Norman 1992) as a base
and added major changes to include the microphysics. First we added an equation for elec-
tron fraction to treat electron captures and neutrino transport by the so-called leakage scheme
(Kotake et al. 2003). Furthermore, we extend the scheme to include all 6 species of neutrino
(νe, ν̄e, νX), which is indispensable for the computations of the Pop III stars. Here νX means
νµ, ν̄µ, ντ and ν̄τ . As for the reactions of νX , pair, photo, and plasma processes are included
using the rates by Itoh et al. (1989). The Lν , in Eq. (3) is the cooling rate of the relevant
neutrino reactions (see Takiwaki et al. 2007, for details). As for the equation of state, we have
incorporated the tabulated one based on relativistic mean field theory instead of the ideal gas
EOS assumed in the original code (Shen et al. 1998).
2.2. Initial Models and Boundary Condition
In this paper, we set the mass of the Pop III star to be 300M⊙. This is consistent with
the recent simulations of the star-formation phenomena in a metal free environment, providing
an initial mass function peaked at masses 100−300M⊙ (see, e.g., Nakamura & Umemura 2001).
We choose the value because we do not treat the nuclear-powered pair instability supernovae
(M<∼260M⊙) and, for convenience, for the comparison with the previous study, which employed
the same stellar mass (Fryer et al. 2001).
We start the collapse simulations of 180M⊙ core of the 300M⊙ star. The core, which
is the initial condition of our simulations, is produced in the following way. According to the
prescription in Bond et al. (1984), we set the polytropic index of the core to n= 3 and assume
that the core is isentropic of ∼ 10kB per nucleon (Fryer et al. 2001) with the constant electron
fraction of Ye = 0.5. We adjust central density to 5× 10
6 g cm−3, by which the temperature
of the central regions become high enough to photodisintegrate the iron (∼ 5× 109K), thus
initiating the collapse. Given central density, the distribution of electron fraction, and entropy,
we construct numerically the hydrostatic structures of the core.
Since we know little of the angular momentum distributions in the cores of Pop III stars
(see, however, Fryer et al. 2001), we add the following rotation profiles in a parametric manner
to the non-rotating core mentioned above. We assume the cylindrical rotation of the core and
change the degree of differential rotation in the following two ways.
1. As for the differential rotation models, we assume the following distribution of the initial
angular velocity,
Ω(X,Z) = Ω0
X2+X20
Z4+Z40
, (6)
where Ω is the angular velocity and Ω0 is the model constant. X and Z denote distance
from rotational axis and the equatorial plane, respectively. We adopt the value of param-
eters, X0 and Z0, as 2×10
8cm,2×109cm, respectively. Since the radius of the outer edge
of the core is taken to be as large as 3.5× 109cm, the above profile represents that the
cores rotate strongly differentially.
2. As for the rigid rotation models, the initial angular velocity is given by,
Ω(X,Z) = Ω0. (7)
As for the initial configuration of the magnetic fields, we assume that the field is nearly
uniform and parallel to the rotational axis in the core and dipolar outside (see Figure 1). For
the purpose, we consider the following effective vector potential,
Ar = Aθ = 0, (8)
r3+ r30
r sinθ, (9)
where Ar,θ,φ is the vector potential in the r,θ,φ direction, respectively, r is the radius, r0 is the
radius of the core, and B0 is the model constant. In this study, we adopt the value of r0 as
3.5× 109 cm. This vector potential can produce the uniform magnetic fields when r is small
compared with r0, and the dipole magnetic fields for vice versa. We set the outflow boundary
conditions for the magnetic fields at the outer boundary of the calculated regions. It is noted
that this is a far better way than the loop current method for constructing the dipole magnetic
fields (LeBlanc & Wilson 1970), because our method produces no divergence of the magnetic
fields near the loop current.
Fig. 1. The configuration of the initial magnetic fields. Note that B0 = 10
12G for this figure. The arrows
represent the vector of the poloidal magnetic fields. The contour shows the logarithm of the magnetic
pressure (: B2/8π).
Changing the initial rotational and magnetic energies by varying the values of Ω0 and B0,
we compute 19 models in this paper, namely, one spherical and 18 magnetorotational models.
In Table 1, we summarize the differences of the initial models. Note that the models are named
after this combination, with the first letters, B12, B11, B10, indicating the strength of initial
magnetic field, the following letter, TW1, TW2, TW4 indicating the initial T/|W | and final
capital letter D or R representing the initial rotational law (D: Differential rotation, R: Rigid
rotation). Note that T/|W | represents the ratio of the rotational to the gravitational energy.
Table 1. Models and Parameters.∗
T/|W |
B0 1% 2% 4%
1010G B10TW1{D,R} B10TW2{D,R} B10TW4{D,R}
1011G B11TW1{D,R} B11TW2{D,R} B11TW4{D,R}
1012G B12TW1{D,R} B12TW2{D,R} B12TW4{D,R}
∗ This table shows the name of the models. In the table they are labeled by
the strength of the initial magnetic field and rotation. T/|W | represents
the ratio of the rotational to the gravitational energy. B0 represents the
strength of the initial magnetic field.
In this paper we assume a BH is formed when the condition
6Gm(r)
>r is satisfied, where
c,G,m(r) are the speed of light, the gravitational constant and the mass coordinate, respectively.
This condition means that we assume that fluids cannot escape from the inner region below the
radius of the marginally stable orbit of a Schwarzschild BH. When this condition is satisfied,
we excise the region inside the radius calculated and then treat it as an absorbing boundary.
Afterwards, we enlarge the boundary of the excised region to take into account the growth of
the mass infalling into the central region. Although it is not accurate at all to refer the central
region as the BH, we cling to the simplification in this paper in order to follow and see the
dynamics later on.
3. RESULT
3.1. Spherical Collapse
First of all, we briefly describe the hydrodynamic features of spherical collapse as a
baseline for the MHD models mentioned later. Note in the following that by “massive stars”,
we mean the stars of ≈ O(10)M⊙ with the initial composition of the solar metallicity, which
are considered to explode as supernovae at their ends of the evolution (Heger et al. 2003).
As in the case of massive stars, the gravitational collapse is triggered by the electron-
capture reactions and the photodisintegration of iron nuclei. On the other hand, the gravita-
tional contraction is stopped not by the nuclear forces as in the case of massive stars but by
the (gradient of) thermal pressure. This is because the progenitor of Pop III stars has high
entropy, i.e. high temperature. We call this bounce as “thermal bounce” for convenience. The
evolution of density, temperature, entropy and radial velocity around the thermal bounce are
shown in Figure 2. Unlike the case of massive stars, no outgoing shock propagates outward
after the thermal bounce. At the bounce, the size of the inner core, which is 200 km in radius
and 6M⊙ in the masscoordinate, grows gradually due to the mass accretion. As seen from the
figure, the materials in the accreting shock regions obtain higher entropy and temperature than
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
-1×1010
-8×109
-6×109
-4×109
-2×109
 2×109
 0  500  1000  1500  2000
Radius [km]
-1×1010
-8×109
-6×109
-4×109
-2×109
 2×109
 0  500  1000  1500  2000
Radius [km]
-1×1010
-8×109
-6×109
-4×109
-2×109
 2×109
 0  500  1000  1500  2000
Radius [km]
-1×1010
-8×109
-6×109
-4×109
-2×109
 2×109
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
 0  500  1000  1500  2000
Radius [km]
Fig. 2. The evolutions of density, temperature, entropy and radial velocities in the spherical model. Solid
line is for -37 ms from bounce, dashed line is for -1 ms, dotted line is for 19 ms and dashed dotted line is
for 32 ms, respectively.
the ones in the case of massive stars.
The higher temperatures are good for producing a large amount of µ- and τ - neutrinos
through the pair annihilation of electrons and positrons. This also makes different features
of the neutrino emissions from the case of massive stars, in which the electron-neutrino (νe)
luminosity dominates over those of the other species near the epoch of core bounce. As shown
in the top panel of Figure 3, the total luminosity of νX (νµ, ντ , ν̄µ and ν̄τ ) begins to dominate
over the total luminosity of electron neutrinos and anti-electron neutrinos at 25 msec after
bounce (see the first intersection of the lines in the figure).
At 87 msec after the bounce, the core is so heavy that it promptly collapses to a BH. In
Figure 4, we show the evolution of the BH calculated by the procedure described in §2.2. The
mass of the BH is initially 20M⊙, rapidly increases to 35 M⊙ because the rest of the dense inner
core falls into the BH soon after the formation. The growth rate of the mass is slowed down
afterwards when the quasi-steady accretion flow to the BH is established (Fig 4). The rapid
decrease of the neutrino luminosity ∼100 msec after bounce (top panel of Figure 3) corresponds
to the epoch when the neutrinospheres are swallowed into the BH. Note in the bottom panel,
-0.2 -0.1  0  0.1  0.2  0.3  0.4  0.5
Time After Bounce [s]
-0.2 -0.1  0  0.1  0.2  0.3  0.4  0.5
Time After Bounce [s]
Fig. 3. Upper panel: Time evolutions of neutrino luminosity in the spherical model. The time is measured
from the thermal bounce. Solid line represents the total luminosity of electron neutrinos and anti-electron
neutrinos. Dashed line represents the total luminosity of νX (νµ, ντ , ν̄µ and ν̄τ ) neutrino. Before the
thermal bounce, the luminosity of electron + anti-electron neutrinos dominate that of νX luminosity,
though, after the bounce, it reverses. At ∼ 0.1 second after the bounce, the luminosities drastically
decrease due to the BH formation. Lower panel: Time-integrated neutrino luminosities. Solid line and
Dashed line are total energy emitted by electron neutrinos and anti-electron neutrinos and X neutrinos,
respectively.
-0.2  0  0.2  0.4  0.6  0.8
Time [s]
Spherical
Fig. 4. Evolution of the BH mass in the spherical model. The time is measured from the BH formation.
The black filled circle indicates the epoch of the black-hole formation.
the total emitted energy are calculated by
dtLν , which represents the energy carried out from
the core by neutrinos. Again from the quantity, it is shown that νX deprives dominantly of the
gravitational energy of the core than νe and ν̄e.
3.2. Rotational Collapse
3.2.1. Effect of Differential Rotation
Now we move on to discuss the features in the rotational core-collapse. The deviation of
the dynamics from the spherical collapse comes from the initial rotation rates and the degree
of the differential rotation initially imposed. To see the effects of the differential rotation on
the collapse-dynamics, we first take models of B10TW1D (differential rotation) and B10TW1R
(rigid rotation) as examples and mention the difference of them. The effects of the initial
rotation rates are discussed later in §3.2.2.
We first describe the collapse of model B10TW1D. As in the case of spherical collapse, the
rotating core experiences the collapse due to the neutrino emission and the photodisintegration,
but the difference appears at the time of the thermal bounce. Due to the pressure support
supplied by the centrifugal force, model B10TW1D bounces at the pole at the epoch 17 msec
later than that of the spherical collapse. The time evolutions after bounce is presented in Figure
5. It is shown that the materials of the inner core oscillate about 20 msec after bounce (see
from the top left down to the bottom), and then the shock wave begins to propagate along the
rotational axis (see from the top right down to the bottom). This jet-like shock wave finally
stalls at Z ∼ 2×108 cm, where Z is the distance from the center along the rotational axis. It is
noted that the shock wave formed at bounce does not stall in the strongly magnetized models
as discussed in §3.3. In this weakly magnetized model, the stellar mantle just collapses to the
central region after the shock-stall, and then leading to the formation of the BH. In this model,
we follow the hydrodynamics until more than 99 % of the materials outside collapse to the BH
(typically 2 sec after bounce).
Fig. 5. Entropy profiles of differential rotation model of B10TW1D 50 (left top), 63 (left middle), 73 (left
bottom), 87 (right top), 113 (right middle), and 127 (right bottom) ms after bounce, respectively. The
color coded contour shows the logarithm of entropy (kB) per nucleon and arrows represent the velocity
fields.
Model B10TW1R thermally bounces rather isotropically in the center, not like model
B10TW1D. This is because the central regions have less angular momentum in comparison with
the differentially rotating model of B10TW1D. In Figure 6, the time evolutions of entropy after
bounce are shown. Unlike B10TW1D (Figure 5), B10TW1R directly collapses to form the BH
without producing the outgoing shock waves. This is because the central part has less pressure
support from the centrifugal force due to the uniform rotation profile initially imposed. On the
other hand, the model B10TW1R has more angular momentum than that of model B10TW1D
in the outer part of the core. This leads to the suppression of the accretion rates of the infalling
matter to the inner core. As a result, the core of model B10TW1R oscillates in a longer period
than that of model B10TW1D because the dynamical timescale, which is proportional to ρ−1/2,
becomes longer due to the smaller density there.
Fig. 6. Same as Figure 5 but for the rigidly rotating model of B10TW1R at 56 (left top), 76 (left bottom),
106 (right top), and 134 (right bottom) ms after bounce, respectively. Color coded contour shows the
logarithm of entropy (kB) per nucleon and arrows represent the velocity fields.
Next we compare the masses of the BH at the formation and the subsequent growth
between the two models (see Figure 7). The initial mass of BH of models B10TW1R and
B10TW1D are 40 and 70 M⊙, respectively (see the black filled circles in the Figure). Both
of them are larger than that of the spherical collapse model (∼ 20 M⊙). As mentioned, the
reason that the earlier formation of less massive BH of model B10TW1R is that the model
has smaller centrifugal forces in the central regions than model B10TW1D. On the other hand,
reflecting the smaller mass accretion rates to the BH, the growth rate of BH’s mass of the model
B10TW1R is smaller than that of the model B10TW1D (compare the slopes of the lines in the
Figure after the BH formation).
The luminosity of neutrinos and the total leaked energy of the model B10TW1D are
shown in the left panel of Figure 8. It is found that the luminosity of µ and τ neutrinos (νX) do
not overwhelm that of the electron neutrinos even after the bounce unlike the spherical collapse
and that most energy are emitted by electron neutrinos (bottom panel). This is because the
rotation suppresses the compression of the core, which lowers the temperature in the central
regions than that of the spherical model. It should be noted that the energy production rates
by the pair annihilation processes sharply depend on the temperature. The neutrino features
of B10TW1R are found to be intermediate between the model B10TW1D and the spherical
model (see right panels).
If the initial rotation rates of the above two models become larger, the bounce occurs
-0.2  0  0.2  0.4  0.6  0.8
Time[s]
B10TW1D
B10TW1R
Fig. 7. Evolution of the masses of BH for the rotating models. Solid and dashed lines are for models
B10TW1D and B10TW1R, respectively. The black circle indicates the epoch of the black-hole formation.
Note that the time is measured from the epoch of the BH formation in the spherical model.
more later due to the stronger centrifugal forces. In addition, the interval of the core oscilla-
tions becomes longer. Except for such differences, the hydrodynamic features before the BH
formation are mainly determined by the degree of the differential rotation as mentioned above
and no qualitative changes are found as the initial rotation rates become larger.
3.2.2. Effects of Rotation on the BH mass and Neutrino emission
In this section, we proceed to describe how the initial rotation rate and the degree of
the differential rotation affect the growth of the BH masses and the neutrino emissions.
The effects of rotation on the initial masses of the BHs for the almost purely rotating
models, labeled by B10, are shown in Figure 9. As seen, larger the initial rotation rate becomes,
the heavier BH is found to be produced. This tendency is independent of the degree of the
differential rotation. This is simply because rapid rotation tends to halt the infall of the matter
to the center, thus heavier masses are required to fulfill the condition of BH formation. It
is furthermore found that the initial mass is larger for the differential rotation models than
the rigid rotation models. This is regardless of the initial rotation rates. This is due to the
smaller angular momentum of the rigid rotation models in the central regions than that of the
differential rotation models as mentioned.
In Figure 10, the growth of the BH mass for the corresponding models is shown. It is
found that the epoch of the formation is delayed as the initial rotation rates become larger
regardless of the degree of the differential rotation. As for the growth rates of BH’s mass, it
is found that they are almost the same for the differential rotation models regardless of the
initial rotation rates (see the left panel of Figure 10). This is because the outer part of the core
has little angular momentum due to the strong differential rotation imposed, and thus falls to
the center in the similar way. On the other hand, the initial rotation rates affect the evolution
of BHs in the rigid rotation models (see the right panel of Figure 10). As the initial rotation
-0.2 -0.1  0  0.1  0.2  0.3  0.4  0.5
Time After Thermal Bounce [s]
-0.2 -0.1  0  0.1  0.2  0.3  0.4  0.5
Time After Thermal Bounce [s]
-0.2 -0.1  0  0.1  0.2  0.3  0.4  0.5
Time After Thermal Bounce [s]
-0.2 -0.1  0  0.1  0.2  0.3  0.4  0.5
Time After Thermal Bounce [s]
Fig. 8. Same as Figure 3 but for models B10TW1D (left) and B10TW1R (right), respectively.
 0  1  2  3  4  5
T/|W| [%]
Differential Rotation
Rigid Rotation
Fig. 9. Effects of initial rotation rate and the degree of differential rotation on the initial mass of the BH.
In this figure, the sequence of the models labeled by “B10”, which are almost purely rotating model, is
chosen. Note that model B10TW4R is absent because this model does not produce the BH during the
simulation time.
rates become larger, the growth rates of the BHs become smaller due to the larger angular
momentum imposed initially.
-0.2  0  0.2  0.4  0.6  0.8
Time[s]
B10TW1D
B10TW2D
B10TW4D
Spherical
-0.2  0  0.2  0.4  0.6  0.8
Time[s]
B10TW1R
B10TW2R
Spherical
Fig. 10. Time evolution of the BH mass for the almost purely rotating models labeled by B10. Left
panel: Solid, dashed, dotted, and dashed-dotted lines, are for models B10TW1D, B10TW2D, B10TW4D,
and the spherical model, respectively. Right panel: Solid, dashed, dotted lines are for models B10TW1R,
B10TW2R, and the spherical model, respectively. The time is measured from the thermal bounce of each
model. The black filled circles of each panel represent the epoch of the BH formation.
BH’s mass at the formation affects the total energy emitted by the neutrinos because
the neutrinos in the region of BH cannot escape to the outside of the core afterwards. The
total energy emitted by neutrinos are shown in Figure 11. One can see the general trend in the
figure that the emitted energy rapidly rises and then becomes constant. The transition to the
constant phase corresponds to the formation of the BH.
Also in this case, differential rotation models have similar features after the formation
of the BHs (solid line of right panel of Figure 11). It is interesting that the model B10TW2R
(dashed line of right panel of Figure 11) by contrast has different behaviors of total emitted
energy. This is because this model produces the stable accretion disk around the central BH. As
a result, the materials of the disk accretes only slowly to the BH, and thus can emit neutrinos
for a longer time. The small accretion rate of this model is also prominent as seen in the Figure
It is interesting to note that only about 10% of the gravitational energy of the core can be
carried away by neutrinos even in the most rapidly rotating model considered here (B10TW4D).
On the other hand, it is well known that neutrinos carry away 99 % of the gravitational energy
of the protoneutron stars in case of the massive stars. The discrepancy stems obviously from
the fact that most part of the inner core is absorbed to the BH in case of the Pop III stars.
3.3. Magnetorotational Collapse
In this section, we present the results of the MHD models. First of all, we mention the
magnetohydrodynamic (MHD) features in section 3.3.1, then discuss the MHD effects on the
 2×1054
 4×1054
 6×1054
 8×1054
 1×1055
-0.4 -0.2  0  0.2  0.4
Time After Bounce [s]
B10TW1D
B10TW2D
B10TW4D
Spherical
 2×1054
 4×1054
 6×1054
 8×1054
 1×1055
-0.4 -0.2  0  0.2  0.4
Time After Bounce [s]
B10TW1R
B10TW2R
Spherical
Fig. 11. Time-integrated neutrino leakage energy. Left panel: Solid line means B10TW1D, dashed line
means B10TW2D, dotted line means B10TW4D and dash-dotted line means Spherical model. Right panel:
Solid line means B10TW1R, dashed line means B10TW2R and dotted line means Spherical model. The
time is measured from the thermal bounce of each model.
BH mass and neutrino emissions in section 3.3.2.
3.3.1. MHD Feature
Amongst the computed models, we find that the models with the strongest magnetic
field (B = 1012G) can only produce the jet-like shock waves along the rotational axis, which
can propagate outside of the core without shock-stall. First of all, we mention the properties
of such models taking model B12TW1D as an example.
The collapse dynamics before bounce is almost the same as the corresponding weak
magnetic field model of B10TW1D. This is because the amplified magnetic fields by the com-
pression and the field-wrapping are, of course larger than the weaker field model, but still much
smaller than the matter pressure in the central regions. After the bounce, the toroidal mag-
netic fields produced by the wrapping, provide the additional pressure support, thus acting to
push the infalling matter as jetlike outflow rather than rotates along the magnetic field. The
jet is launched when the magnetic pressure overcomes the local ram pressure of the accreting
matter. This feature is different from another jet driving mechanism, the magneto-centrifugal
acceleration (Blandford & Payne 1982).
The MHD features of model B12TW1D after bounce are presented in Figure 12. From
the right panels, it is shown that the regions behind the jet-like shock wave (Z ≥ 1.5×109 cm)
become dilute with the density of ρ ∼ 105 g cm−3 and have very high entropy s ∼ 102kB.
The bottom panel shows the jet is driven by magnetic pressure because the plasma beta (≡
gas pressure / magnetic pressure) of the region of jet is much smaller than unity. As the
jet propagates in the core, a newborn BH is produced (see the white circles of right panels of
Figure 12). The mass of the BH is initially 57.9M⊙, which is smaller than the one of B10TW1D
(70.4M⊙). The reason of the difference is mentioned in §3.3.2.
The properties of jet of model B12TW1D are shown in Figure 13. There are profiles of
Fig. 12. Time evolution of shock waves of the strongest magnetized model of B12TW1D. The top panel of
the figure shows the logarithm of entropy (kB) per nucleon, the middle panel of the figure shows logarithm
of density (g cm−3), and bottom panel shows the logarithm of plasma beta. All left figures are at 119 ms
from bounce, and the right figures are at 305 ms. The white circles of the right panels represent the BHs.
density, radial velocity, magnetic field, and pressure at 104 ms after bounce. The density of
matter in the jet region is ∼ 107 g cm−3. The speed of shock front is as large as 40 % of speed of
light, which is mildly relativistic. It is easily seen that the toroidal magnetic field overwhelms
the poloidal component behind the shock front. In the inner region, the poloidal magnetic field
is larger due to compression. The magnetic pressure overwhelms gas pressure throughout the
jet region as already depicted in Figure 12.
Now we move on to discuss how rotation affects the dynamics while fixing the initial
field strength. Figure 14 shows the properties of models B12TW1{D,R} and B12TW4{D,R}
when the jet-like shock wave reaches to 1×109 cm. As clearly seen, the main difference between
T/|W | = 1 % and 4 % is the degree of the collimation of the shock wave. As the initial rotation
rates become large in the differential rotation models, the compression of the magnetic fields
is hindered, thus leading to the suppression of the hoop stress of magnetic fields in the central
regions (Takiwaki et al. 2004). As a result, the collimation of the shock wave becomes less
 0  1×109  2×109
Radius [cm]
-4×1010
-2×1010
 2×1010
 0  1×109  2×109
Radius [cm]
 0  1×109  2×109
Radius [cm]
poloidal
toroidal
 0  1×109  2×109
Radius [cm]
magnetic pressure
gas pressure
Fig. 13. Various physical quantities around the rotational axis at 104 ms after bounce for model
B12TW1D. Density (left top), radial velocity (right top), absolute value of magnetic field (left bottom),
and pressure are shown. In the left bottom panel, the solid line and dashed line represent poloidal com-
ponent and toroidal component, respectively. In the right bottom panel, the dashed line represents gas
pressure and solid line represents magnetic pressure.
(compare the top panels). In contrast, the difference of the degree of collimation between
models B12TW1R and B12TW4R is smaller than B12TW1D and B12TW4D. This is because
the materials in the inner region of rigid rotating models rotate more slowly than the differential
models. Thus the degree of the collimation of the shocks depends weekly on the initial T/|W |.
Models with the weaker initial magnetic fields do not produce the jet-like explosion
except for model B11TW4R. These models collapse to BHs before formation of jets because
of weak magnetic pressure. After forming a BH, especially differentially rotating model, rest
parts of star rotate slowly so that the magnetic pressure does not grow up. Thus, when BH
is formed before jet rises, rest of core only collapse to BH and entire the core are absorbed by
3.3.2. MHD Effects on the BH Mass and Neutrino Emission
The MHD effects on the initial masses of BHs are shown in Tables 2 and 3. As seen,
the initial BH mass gets smaller when the initial magnetic field becomes stronger. The angular
Fig. 14. Profiles of the shock propagation for models B12TW1D (top left), B12TW4D (top right),
B12TW1R (bottom left), and B12TW4R (bottom right), respectively. They show the color coded contour
plots of logarithm of entropy (kB) per nucleon. Various profiles are found by changing the strength of the
initial magnetic field and rotation.
momentum transport by the magnetic fields is an important agent to affect the BH mass.
This feature is seen in Figure 15, which represents the distribution of mean specific angular
momentum of models B12TW1D and B10TW1D. The central region of B12TW1D has smaller
angular momentum than B10TW1D due to angular momentum transport by magnetic fields.
The peak of model B12TW1D represents the position of shock front on the equatorial plane.
The transport of the angular momentum makes centrifugal force of central region smaller and
enhances the collapse. This leads the BH mass smaller. This tendency is more prominent for
the rigid rotation models (compare Tables 2 and 3) because the rotation of the central region
is slower than differential rotation models and contraction of core is more significant, which
leads the amplification of magnetic field and larger angular momentum transport. Figure 16
shows the relation between BH mass and angular momentum at BH formation. The angular
momentum of BH gets larger with its mass at the time of BH formation. This is because the
matter with large angular momentum cannot collapse due to centrifugal force and requires large
Table 2. Initial Mass of Black Holes for Differentially Rotating Models [M⊙].
T/|W |
B0 1% 2% 4%
1010G 70.4 87.3 106.6
1011G 70.4 87.3 106.6
1012G 57.9 75.8 96.6
Table 3. Initial Mass of Black Holes for Rigidly Rotating Models [M⊙].
T/|W |
B0 1% 2% 4%
1010G 40.5 75.8 —
1011G 38.6 38.6 —
1012G 15.3 15.1 15.1
 0  20  40  60  80  100  120
Mass Coordinate [M
B12TW1D
B10TW1D
Fig. 15. Mean specific angular momentum over the shells as a function of the mass coordinate just before
BH formation. The solid line and dotted line represent model B12TW1D and B10TW1D, respectively.
amount of accreting matter to collapse BH, as already mentioned in §3.2.
Next, we discuss the MHD effects of neutrino emissions. Figure 17 shows the peak
neutrino luminosities as a function of the initial rotation rates. It is shown that the magnetic
fields make the peak luminosities smaller, when fixing the initial degree of the differential
rotation (compare B10D and B12D, and B10R and B12R). It is noted that the Pop III stars
have gentle slope of the density prior to core-collapse, so that materials of the outer region
have a great deal of the total gravitational energy of the iron core. Thus the stronger magnetic
pressures, which prevent the accretion, make the liberating gravitational energy of the accreting
matter smaller, and thus, results in the suppression of the peak luminosities. In case of the rigid
rotation, the stronger centrifugal forces in the outer regions, lead to the stronger suppression
of the releasable gravitational energy than in the case of the differential rotation. As a result,
 10  20  30  40  50  60  70  80  90  100  110
BH Mass [M
Fig. 16. Relation between BH mass and angular momentum at the BH formation.
the peak luminosities for the rigidly rotating models decreases more steeply with the initial
rotation rates than the ones for the differentially rotating models (compare B10R and B10D).
 2×1055
 4×1055
 6×1055
 0  1  2  3  4  5
T/|W| [%]
Fig. 17. Effects of rotation and magnetic fields on the peak luminosity of neutrinos.
4. SUMMARY AND DISCUSSION
We studied the magnetorotational core-collapse of Pop III stars by performing the two-
dimensional magnetohydrodynamic simulations. Since the distributions of rotation and mag-
netic fields in the progenitors of Pop III stars are highly uncertain, we changed them in a
parametric manner and systematically investigated how rotation and magnetic fields affect the
dynamics of Pop III stars. In addition, we explored how rotation and magnetic fields affect
the formation of the BHs and the neutrino emissions. In the current Newtonian simulations,
the BH formation was ascribed to a certain condition, and after the formation, the central
region was excised and treated as an absorbing boundary. As for the microphysics, we took
into account the neutrino cooling of 6 species by a leakage scheme with a realistic equation of
state. With these computations, we have obtained the following results,
1. In the spherical model, the gravitational contraction is stopped by the gradient of the ther-
mal pressure, not by the nuclear forces as in the case of the massive stars ≈O(10)M⊙ with
the initial composition of the solar metallicity, because the progenitor of Pop III stars has
high entropy, i.e. high temperature initially. Such high temperature also makes different
features of the neutrino emissions from the case of the massive stars. The luminosity of
µ- and τ - neutrinos dominates over that of the electron neutrinos after core bounce. Thus
the gravitational energy of the core is carried away dominantly by µ- and τ - neutrinos.
2. As the initial rotation rates of the core become larger, it is found that the epoch of the BH
formation is delayed later and that the initial masses of the BHs become larger. Fixing the
initial rotational energy, the BH masses at the formation become larger as the degree of
the differential rotation becomes stronger. As the initial degree of the differential rotation
becomes larger, the electron neutrino luminosity is found to be more dominant over that
of µ and τ neutrinos after core bounce, because the pair creation processes of µ and τ ,
which sharply depend on the temperature, are more suppressed.
3. We find that the jet-like explosions can be produced even in Pop III stars if the magnetic
field is as large as 1012G prior to core-collapse. This jet-like shock wave is completely
magneto-driven.
4. Jet-like shocks in the stronger magnetic field models are themselves found to make the
initial mass of the BH smaller. The angular momentum transport by magnetic fields is
found to be an important agent to make the initial mass of BHs smaller because the
transport of the angular momentum enhances the collapse of the central regions. As a
result, it is found that the initial BH masses for the most strongly magnetized models
are found to become smallest when fixing the initial rotation rates. As for the neutrino
luminosities, we point out that the stronger magnetic fields make the peak luminosities
smaller, because they can halt the collapse of the materials.
Here we shall refer to the limitations of this study. First of all, we mimicked the neutrino
transfer by the leakage scheme. Although the scheme is a radical simplification, we checked
that we could reproduce, at least, the qualitative features of the neutrino luminosities. The
supremacy of νX neutrinos’ luminosity in the spherical collapse of Pop III stars obtained in
this study is consistent with the foregoing studies by Fryer et al. (2001) and Nakazato et al.
(2006), in which the more elaborate neutrino transport schemes were employed. Furthermore,
the effects of rotation on the emergent neutrino luminosities are consistent with Fryer & Heger
(2000), in which one model of the rotational collapse of the massive stars was investigated.
Secondly, the simulations were done with the Newtonian approximation and we defined the
BH formation by the marginally stable orbit of a Schwarzschild BH. This treatment is totally
inaccurate because the core rotates so rapidly that fully general relativistic (MHD) simulations
with the appropriate implementations of the microphysics are necessary, however, are still too
computationally prohibitive and beyond our scope of this paper. Remembering these caveats,
this calculation is nothing but a demonstration showing how the combinations of rotation and
magnetic fields could produce the variety of the dynamics, and the important consequences in
the properties of the neutrino emissions, and these outcomes, of course, should be re-examined
by the more sophisticated simulations.
In this study, we followed the dynamics till ∼ 1 sec after the formation of the BHs
and saw the shock break-out from the cores in the strongly magnetized models. But if we
were to follow the dynamics in the much later phase in the weaker magnetized models, the
magneto-driven outflows might be produced, due to the long-term field-wrapping and/or the
development of the so-called magnetorotational instability, as demonstrated in the study of
collapsar (see, e.g., Proga et al. 2003; Fujimoto et al. 2006). But it should be noted that the
dynamical phases considered here and the other ones are apparently different (and thus they
are complimentary). In the latter studies, the central BHs with a rotationally supported disk
around are treated as an initial condition for the computations. Our core-collapse simulations
presented here showed that the outer region rotates much slowly than the Keplerian one and
most of them directly collapses to BH. Thus the amplification of the magnetic field in the disks,
which needs rapid rotation, might not be so efficient as previously demonstrated. To clarify it,
we are now preparing for the long-term simulations, in which the final states obtained here are
taken as an initial condition.
Then we discuss the validity of the initial strength of the magnetic fields assumed in this
study. For the purpose, we estimate the strength of the magnetic field just before collapse with
Eq. (13) of Maki & Susa (2004), in which the thermal history of the primordial collapsing clouds
was calculated in order to investigate the coupling of the magnetic field with the primordial
gas. For example, Bini∼10
−7 G and nH,ini∼10
3 cm−3, which are the values they employed, lead
B ∼ 1011 G if the magnetic flux is conserved during the contraction and the clouds collapse to
106 g cm−3. Although the above parameters chosen are slightly optimistic, the magnetic fields
assumed in this study may not be so unrealistic.
We pointed out that the total neutrino energy emitted from rotation models increases
several times than the one of the spherical collapse model. However, the detection of such
neutrinos as the diffusive backgrounds might be difficult because the Pop III stars are too
distant (see Iocco et al. 2005). Alternatively, the detection of gravitational waves from Pop III
stars as the backgrounds seems more likely (Buonanno et al. 2005; Sandick et al. 2006) by the
currently planning air-borne laser interferometers such as LISA2, DECIGO (Seto et al. 2001)
and BBO (Ungarelli et al. 2005), and needs further investigation.
We found that the Pop III stars are able to produce jet-like explosions with mass ejections
when the central cores are strongly magnetized. This may be important with respect to its
relevance to the nucleosynthesis in such objects (Ohkubo et al. 2006). This is also an interesting
topic to be investigated as a sequel of this paper.
2 http://lisa.jpl.nasa.gov
This study was supported in part by the Japan Society for Promotion of Science (JSPS)
Research Fellowships (T.T.), Grants-in-Aid for the Scientific Research from the Ministry of
Education, Science and Culture of Japan (No.S14102004, No.14079202, No.17540267, and No.
1840044). Numerical computations were in part carried out on VPP5000 at the Center for
Computational Astrophysics, CfCA, of the National Astronomical Observatory of Japan.
References
Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93
Akiyama, S., Wheeler, J. C., Meier, D. L., & Lichtenstadt, I. 2003, ApJ, 584, 954
Ando, S. & Sato, K. 2004, New Journal of Physics, 6, 170
Ardeljan, N. V., Bisnovatyi-Kogan, G. S., & Moiseenko, S. G. 2005, MNRAS, 359, 333
Barkana, R. & Loeb, A. 2001, Phys. Rep., 349, 125
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825
Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23
Bromm, V. & Larson, R. B. 2004, ARA&A, 42, 79
Bromm, V. & Loeb, A. 2006, ApJ, 642, 382
Buonanno, A., Sigl, G., Raffelt, G. G., Janka, H.-T., & Müller, E. 2005, Phys. Rev. D, 72, 084001
Christlieb, N., et al. 2002, Nature, 419, 904
Daigne, F., Olive, K. A., Vangioni, E., Silk, J., & Audouze, J. 2004, ApJ, 617, 693
Frebel, A., et al. 2005, Nature, 434, 871
Fryer, C. L. & Heger, A. 2000, ApJ, 541, 1033
Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372
Fujimoto, S.-i., Kotake, K., Yamada, S., Hashimoto, M.-a., & Sato, K. 2006, ApJ, 644, 1040
Glover, S. 2005, Space Sci. Rev., 117, 445
Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288
Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
Iocco, F., Mangano, G., Miele, G., Raffelt, G. G., & Serpico, P. D. 2005, Astroparticle Physics, 23,
Itoh, N., Adachi, T., Nakagawa, M., Kohyama, Y., & Munakata, H. 1989, ApJ, 339, 354
Iwamoto, N., Umeda, H., Tominaga, N., Nomoto, K., & Maeda, K. 2005, Science, 309, 451
Kotake, K., Sato, K., & Takahashi, K. 2006, Rep. Prog. Phys., 69, 971
Kotake, K., Sawai, H., Yamada, S., & Sato, K. 2004a, ApJ, 608, 391
Kotake, K., Yamada, S., & Sato, K. 2003, Phys. Rev. D, 68, 044023
Kotake, K., Yamada, S., Sato, K., Sumiyoshi, K., Ono, H., & Suzuki, H. 2004b, Phys. Rev. D, 69,
124004
LeBlanc, J. M. & Wilson, J. R. 1970, ApJ, 161, 541
Maki, H. & Susa, H. 2004, ApJ, 609, 467
Nakamura, F. & Umemura, M. 2001, ApJ, 548, 19
Nakazato, K., Sumiyoshi, K., & Yamada, S. 2006, ApJ, 645, 519
Obergaulinger, M., Aloy, M. A., & Müller, E. 2006, A&A, 450, 1107
Ohkubo, T., Umeda, H., Maeda, K., Nomoto, K., Suzuki, T., Tsuruta, S., & Rees, M. J. 2006, ApJ,
645, 1352
Proga, D., MacFadyen, A. I., Armitage, P. J., & Begelman, M. C. 2003, ApJL, 599, L5
Sandick, P., Olive, K. A., Daigne, F., & Vangioni, E. 2006, Phys. Rev. D, 73, 104024
Sawai, H., Kotake, K., & Yamada, S. 2005, ApJ, 631, 446
Scannapieco, E., Madau, P., Woosley, S., Heger, A., & Ferrara, A. 2005, ApJ, 633, 1031
Schneider, R., Guetta, D., & Ferrara, A. 2002, MNRAS, 334, 173
Seto, N., Kawamura, S., & Nakamura, T. 2001, Phys. Rev. Lett., 87, 221103
Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 1998, Nucl. Phys. A, 637, 435
Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753
Takiwaki, T., Kotake, K., Nagataki, S., & Sato, K. 2004, ApJ, 616, 1086
Takiwaki, T., Kotake, K., Yamada, S., & Sato, K. 2007, in preparation
Umeda, H. & Nomoto, K. 2002, ApJ, 565, 385
—. 2003, Nature, 422, 871
Ungarelli, C., Corasaniti, P., Mercer, R., & Vecchio, A. 2005, Class. Quant. Grav., 22, S955
Weinmann, S. M. & Lilly, S. J. 2005, ApJ, 624, 526
Yamada, S. & Sawai, H. 2004, ApJ, 608, 907
