arXiv:0704.1853v1  [astro-ph]  14 Apr 2007
PASJ: Publ. Astron. Soc. Japan , 1–??, 〈publication date〉
c© 2019. Astronomical Society of Japan.
Dissipation of Magnetic Flux in Primordial Star Formation:
From Run-away Phase to Mass Accretion Phase
Hideki Maki
Department of Physics, Rikkyo University, Nishi-Ikebukuro, Tokyo 171-8501
hide-mk@jcom.home.ne.jp
Hajime Susa
Department of Physics, Konan University, Okamoto, Kobe 658-8501
susa@center.konan-u.ac.jp
(Received 〈reception date〉; accepted 〈acception date〉)
Abstract
We investigate the dissipation of magnetic flux in primordial star-forming clouds throughout their
collapse including the run-away collapse phase as well as the accretion phase. We solve the energy equation
and the non-equilibrium chemical reactions in the collapsing gas, in order to obtain the radial distribution
of the ionized fraction during the collapse. As a result, we find the ionized fraction is high enough for the
magnetic field to couple with the gas throughout the evolution of the cloud. This result suggests that the
jet formation from protostars as well as the activation of magneto-rotational instability in the accretion
disk are enabled in the presence of the cosmological seed magnetic flux proposed by Langer, Puget &
Aghanim (2003).
Key words: accretion, accretion disks — diffusion — early universe — stars: formation — stars:
magnetic fields
1. INTRODUCTION
Typical mass or initial mass function of population III stars are fundamental parameters that have great impacts
on subsequent structure formation of the universe. Those stars are expected to be as massive as 100−1000M⊙, which
ionize/dissociate the surrounding media. As a result, star formation activities in the neighbourhood of the stars are
highly regulated (e.g. Susa (2007) and the references are therein).
There are two chief reasons that the population III stars are expected to be very massive. First one comes from
numerical studies by several authors (Nakamura & Umemura 1999; Bromm, Coppi & Larson 1999; Abel, Bryan &
Norman 2000; Nakamura & Umemura 2001; Bromm, Coppi & Larson 2002). Those studies reveal that the prestellar
cores formed as fragments of primordial gas are as massive as ∼ 103M⊙ − 104M⊙. Second ground is the very high
accretion rates of those stars, which is as high as 10−3− 10−2M⊙yr−1. These two facts are the direct consequence of
relatively high temperature (∼ 1000K) of primordial gas because of inefficient cooling by H2.
On the other hand, abundance pattern of the hyper metal-poor stars seems to be more consistent with that of
“faint” supernovae as remnants of less massive (∼ 25M⊙) population III stars(Christlieb et al. 2004; Frebel et al.
2005; Iwamoto et al. 2005). Recent theoretical studies suggest the possibilities for the formation of such less massive
population III stars. If the primordial gas was once ionized, enhanced fraction of H2 causes more efficient cooling and
HD formation. Since HD molecules have lower excitation energy than that of H2, gas can be cooled below 100K. As
a result, the mass of the fragments could be smaller than 100M⊙(O’Shea, Abel, Whalen & Norman 2005; Johnson &
Bromm 2006). O’Shea & Norman (2007) also demonstrates that such low mass population III stars cloud be formed
directly from the cosmological density fluctuations.
Those new ideas are quite interesting and promising, however, we still do not know the actual mass of population III
stars formed from ∼ 103M⊙ − 104M⊙ prestellar cores, which are commonly found in cosmological simulations. Since
they are only 104− 106cm−3, we have to follow the subsequent evolution of the collapsing cloud. Further evolution of
prestellar cores is first investigated by Omukai & Nishi (1998), and they found that the collapse proceeds in a run-away
fashion and converges to Larson-Penston type similarity solution (Larson 1969; Penston 1969; Suto & Silk 1988) with
polytropic index γ ≃ 1.09. They also find the mass accretion rate is very large compared to the present-day forming
stars, although spherical symmetry is assumed in their radiation hydrodynamic simulations. Recently, Yoshida et al.
(2006) have performed cosmological simulations in which the run-away collapsing core is traced up to a very dense
regime (nH ≫ 1010cm−3), taking the Sobolev type line transfer approximations. They also find basically consistent
results in 3D cosmological simulations with previous 1D results by Omukai & Nishi (1998).
However, most of the mass of protostar is accumulated in the accretion phase, which has not been traced especially in
http://arxiv.org/abs/0704.1853v1
2 Maki & Susa [Vol. ,
multi-dimensional simulations. Recent numerical simulations performed by Saigo, Matsumoto & Umemura (2005) show
that disks would be formed at the center of collapsing primordial clouds, which might result in the disks surrounding
the protostars, or binaries. Thus, the actual mass of a population III star should depends on the mechanism of the
angular momentum transport in the accretion disk. There are a few possibilities of the angular momentum transport
mechanism such as gravitational torque by the nonaxisymmetric structures in the accretion disk, the interaction
among the fragments of the disk (Stone et al. 2000; Bodenheimer et al. 2000), and the turbulent viscosity triggered
by Magneto-Rotational Instability (MRI) (Hawley & Balbus 1992; Sano, Inutsuka & Miyama 1998; Sano & Inutsuka
2001). As for MRI induced turbulence, the strength of magnetic field in the disk is the key quantity to activate
the instability(Tan & Mckee 2004; Tan & Blackman 2004). It is also worth noting that recent MHD simulations by
Machida et al. (2006) suggests the possibility of bipolar jets from proto-population III stars, which also could suppress
the mass accretion onto the central core. The formation of jets also requires the presence of some level of magnetic
field. Therefore, it is quite important to assess the magnetic field strength brought into the accretion disk from initial
weak cosmological seed field (e.g. Langer, Puget & Aghanim (2003)).
Maki & Susa (2004) investigated this issue by solving detailed chemical reaction rate equations coupled with energy
equation in run-away collapsing core. They found that magnetic field is always frozen to the collapsing core in case
the strength of initial magnetic field is smaller than 10−5(nH/10
3 cm−3)0.55 G. This is comparable to the maximal
strength of magnetic field which allow the clouds to collapse. Therefore they conclude magnetic field is always frozen
to the collapsing cloud in the run-away collapsing core. They also evaluate the minimal strength of magnetic field
which is required to activate MRI assuming magnetic field is frozen to the gas not only in the run-away phase, but
also in the accretion phase.
This assumption is based upon the argument that the temperature of the accretion flow would rise faster than
that of the core, because of the shock heating. As a result, the ionization degree is expected to be higher than those
in run-away phase, which guarantee the gas to be frozen to the magnetic field lines. However, such heating is only
important in the final phase of accretion, where the flow settle onto the accretion disk. In order to investigate the
coupling of matter and magnetic field in the accretion phase, we need to calculate the actual thermal and chemical
evolution of accretion flow. Aside from the issue on ionization degree, the ambipolar diffusion velocity is proportional
to the square inverse of density (vambB ∝ ρ−2). As a result, the diffusion velocity increases rapidly in the mass accretion
phase, since the density at a fixed radius decreases as the accretion proceeds. Thus, we need a detailed treatment to
test the coupling between the gas and magnetic field in accretion phase.
In this paper, we investigate the dissipation of magnetic field in collapsing primordial gas cloud from run-away
phase to accretion phase, by solving detailed thermal and chemical rate equations. In the next section, we describe
the formulations employed. In §3, results of our calculations are shown. The formation of jets and activation of MRI
is discussed in §4. Final section is devoted to summary.
2. METHOD OF CALCULATIONS
In order to evaluate the coupling of the gas and the magnetic field, we need to assess the amount of ions and electrons
in collapsing primordial cloud. The collapse of the cloud is expected to proceed in run-away fashion (run-away phase)
in the beginning, followed by the mass accretion phase after the formation of rotationally supported disk at the center.
In this paper, we follow the chemical and thermal evolution of the materials to form star, throughout the two phases.
We assume spherically symmetric collapse of the progenitor gas, although we expect the formation of rotationally
supported disk in the very dense regime (nH ≫ 1010cm−3). It is also assumed that the magnetic field is so weak that
the dynamics of the collapsing gas is not affected by the magnetic force.
2.1. Run-away Collapse Phase
The run-away collapse phase is traced by solving 1-dimensional hydrodynamics. The equation of continuity is
= 4πr2ρ, (1)
whereas the equation of motion is
=−4πr2 ∂p
, (2)
where m is the mass within radius r, ρ is the density of the cloud, u is the velocity and p is the pressure. Since
the equation of state of collapsing primordial gas in run-away phase is known to be approximated by polytrope with
γ = 1.09 (Omukai & Nishi 1998), we employ following equation:
p=Kργ , γ = 1.09 (3)
We solve above set of equations by spherically symmetric Lagrangian hydrodynamics code developed by ourselves,
following the Piecewise Parabolic Method (PPM) described in Colella & Woodward (1984).
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 3
2.2. Mass Accretion Phase
The dynamics of accretion phase could be approximated by simple free-fall, since the flow velocity is supersonic.
Basic equations of free-falling accretion flow consists of the equation of continuity and the equation of motion:
4πρ(t0,r0)r
0dr0 = 4πρ(t,r)r
2(t)dr, (4)
d2r (t)
=−GM0
r2(t)
, (5)
where r, t and t0 represent the position of fluid element, time, and the time when the mass accretion starts, respectively.
M0 denotes the mass within r0 at t= t0:
4πr′0ρ0dr
0. (6)
The initial radius r0 is defined as r0 ≡ r(t0). Therefore, the the solution of above equation of motion is given as
r = r(t; t0,r0). The velocity of the fluid element is given by
u(t; t0,r0) =−
2E0 +2
, (7)
where E0 is the total energy defined as E0 ≡ u2(t0)/2−GM0/r0. Using equation (4), the density at (t,r) is given as
ρ(t,r) =
(∂r/∂r0)t
. (8)
(∂r/∂r0)t is the partial differentiation by r0 keeping t fixed. An explicit expression of this term is given in appendix
1. In order to clarify the validity of the free-fall approximation, it is compared with the similarity solution with
γ = 1.09 (Suto & Silk 1988) in accretion phase. The initial condition of the accretion flow is set as the final state of
the run-away collapse with γ = 1.09. Figure 1 illustrates the density and velocity distributions of several snapshots.
It is obvious that the free-fall approximation can describe the matter distribution in accretion phase very well for
polytropic gas. In other words, the thermal evolution has little effects on the dynamics of the accretion flow. Thus,
the dynamics can be approximated by free-fall formula, irrespective of the internal energy equation.
Fig. 1. Time evolution of density/velocity profiles in accretion flow are shown. Dotted lines and marks represent the result with
free-fall approximation, whereas the solid lines denote the results from the similarity solution in accretion phase with γ = 1.09.
In the accretion flow, the effective polytropic index is basically unknown. Thus, we have to solve the following
energy equation for each fluid element:
4 Maki & Susa [Vol. ,
=−p d
−L(net). (9)
In this expression, p denotes the pressure, which is related to the internal energy ε by the equation of state:
p= (γad − 1)ρε, (10)
where γad is the adiabatic exponent, L(net) denotes the net energy loss rate per unit mass.
We take into account the rovibrational line cooling by H2(Lline), the continuum radiation from the gas (Lcont), and
cooling and heating associated with chemical reactions (Ldiss, GH− , GH+
and G3body). L(net) is the sum of all these
contributions.
The cooling rate by the rovibrational transition of hydrogen molecules are assessed by the fitting formula given by
Galli & Palla (1998) in case the lines are optically thin. In the optically thick regime, we employ the escape probability
formalism, as described in equations (12)-(15) in Omukai (2000), except that we use actual velocity gradient in our
calculations.
The cooling rate owing to the continuum radiation is given by
Lcont = 4σκgasT 4 (11)
where σ denotes the Stefan-Boltzmann constant, κgas is the Planck mean opacity of gas (Lenzuni, Chernoff & Salpeter
1991), which includes bound-free absorption, free-free absorption, photodissociation, Rayleigh scattering, and collision-
induced absorption. This formula is valid in case the accreting gas is optically thin for continuum radiation. We
confirmed that the optical depth is smaller than unity throughout the accretion phase.
Cooling due to the latent heat of H2 dissociation is given by
Ldiss = 4.48
eV s−1 g−1, (12)
where nH is the number density of hydrogen nuclei, and (dyH2/dt)− is the dissociation rate of H2.
On the other hand, the gas is heated when a hydrogen molecule is formed. The heating rate per unit mass is given
GH− =
1+ncr/nH
eV s−1 g−1, (13)
1+ncr/nH
eV s−1 g−1, (14)
G3body =
3body
1+ncr/nH
eV s−1 g−1, (15)
where (dyH2/dt)H− , (dyH2/dt)H+
and (dyH2/dt)3body are the formation rates of H2 by H
− process, H+2 process and
three body reactions (Hollenbach & McKee 1979). ncr is the critical density defined as
ncr =
T 1/2
1.6yH exp
+1.4yH exp
1.2× 104
T +1200
cm−3. (16)
2.3. Chemical Reactions
Since the dissipation of the magnetic flux, as well as the temperature of the gas, strongly depends on the chemical
abundances, we have to solve the non-equilibrium chemical reaction rate equations coupled with equation of motion
and energy equation described in previous subsections 2.1 and 2.2.
The evolution of the fraction of species i is followed by solving the equations,
nHklmylym+
nHklmnylymyn,
(i= 1,2,3, · · · ,24), (17)
where yi ≡ ni/nH is the fraction of species i, klm [cm3 s−1] and klmn [cm6 s−1] are the reaction rate coefficients with
respect to two-body processes and three-body processes, respectively. In our calculation, we include 24 species: e−,
H+, H, H−, H2, H
2 , H
3 , D, D
+, D−, HD, HD+, H2D
+, He, He+, He++, HeH+, Li, Li+, Li++, Li3+, Li−, LiH,
and LiH+. We employ the latest reaction rate coefficients in the following papers, Galli & Palla (1998), Omukai
(2000), Stancil, Lepp & Dalgarno (1998), Flower (2002) and Lepp, Stancil & Dalgarno (2002). As for the radiative
recombination, we use the rate coefficients based on Spitzer (1978). These reactions are the same as our previous
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 5
paper(Maki & Susa 2004). Here we stress the importance to include above rare elements such as Li, because the
coupling of magnetic field with gas can be attained by very low fractional abundance of electrons and ions (Maki &
Susa 2004).
2.4. Drift velocity of magnetic field
Magnetic field is dissipated from star-forming gas via ohmic loss and ambipolar diffusion. We assess the drift
velocity vBx of the field lines due to these two processes, which is compared to the accretion velocity of gas. We
evaluate the drift velocity following the formulation by Nakano & Umebayashi (1986). There are two important
quantities which characterize these diffusion processes. They are τν and ων which denote the viscous damping time of
the relative velocity of charged particle ν to the neutral particles, and the cyclotron frequency of the charged particle
ν, respectively. Then, τν is expressed as
µνnnνnn〈σv〉νn
, (18)
where µνn is the reduced mass, nν , nn, and ρν are, the mean number density for the charged particle ν, the neutral
particle n, and the mass density of charged particle ν, respectively. The averaged momentum-transfer rate coefficient
for a particle ν colliding with a neutral particle is expressed by 〈σv〉νn . We use the empirical formulae for the
momentum-transfer rate coefficients (Kamaya & Nishi 2000; Sano et al. 2000).
According to Nakano & Umebayashi (1986), the drift velocity is given by
vBx =
(j ×B)x, (19)
where
A=A21 +A
2, (20)
ρντνω
1+ τ2νω
, (21)
1+ τ2νω
, (22)
B is the mean magnetic field in the primordial cloud, the suffix x means x direction component in local Cartesian
coordinates where the z direction is taken as the direction of B. We replace (1/c)(j×B)x in equation (19) with the
mean magnetic force B2/4πr, where B is the mean field strength in the cloud, r is the radius of the cloud at which
we assess the drift velocity.
2.5. Initial Conditions
We consider primordial star-forming gas clouds that formed in the mini-halos with ∼ 106 M⊙. We set uniform
and spherical gas cloud with nH = 10
3 cm−3, T0 = 250 K, and M = 10
4 M⊙. Such clouds are commonly found
3D cosmological simulations (Bromm, Coppi & Larson 1999; Abel, Bryan & Norman 2000; Abel, Bryan & Norman
2002; Bromm, Coppi & Larson 2002; Yoshida et al. 2003; O’Shea & Norman 2006). We also use the cosmological
abundance given by Galli & Palla (1998), as the initial fraction of the chemical compositions. The initial magnetic
field strength is expected to be very weak, however, its magnitude is still under discussion. Thus, we regard the field
strength as a free parameter of the calculations.
2.6. Switching from Run-away Phase to Accretion Phase
The run-away phase and the subsequent accretion phase have to be treated in different formulation from each other
as described in subsections 2.1 and 2.2. Thus, we have to switch the scheme from the method in run-away phase to
that in accretion phase. In the light of physical arguments, two phases should be switched when the central run-away
collapse is stopped due to thermal pressure or centrifugal force. In our calculation, we switch the scheme when the
sonic point of the accreting flow go inside a certain small radius. We set this radius to be rsw = 700 R⊙, which is
comparable to the disk size when the accreted mass is comparable to ∼ 1M⊙(Tan & Mckee 2004). It is also worth
noting that the central density of the core at switching is ∼ 1016 cm−3.
3. RESULTS
3.1. Runaway Collapse Phase
3.1.1. The distribution of the chemical abundances
Two snapshots of yi distribution during the run-away collapse phase are shown in Figure 2. The left panel shows
the distribution at the time when the central density nH,c reaches 10
12 cm−3, whereas the right panel represents the
6 Maki & Susa [Vol. ,
Fig. 2. The distribution of the main species, e−, H+, H, H2, H
, Li, and Li+ are plotted. The vertical axis denotes the
fractional abundance yi of the above species, and the horizontal axis denotes the radius from the core center. Left panel represents
the snapshot when the central density satisfies nH,c = 10
12 cm−3, whereas right panel shows the results at nH,c = 10
16 cm−3.
Fig. 3. Drift velocity vBx as a function of the radius r and the field strength B at the central density . Contour map of
vBx/u is shown on r − B plane for nH,c = 10
16 cm−3 Solid curves represent the constant loci along which log(vBx/u) equal
to the values labeled on the curve. The dashed curve represents the critical field strength Bcr given by the equation (23).
distribution at nH,c = 10
16 cm−3. The horizontal axis denotes the radius from the center of the primordial gas cloud
r [cm], and the vertical axis is the fractional abundances yi of each species. Since the central density increases as the
collapse proceeds, right panel represents the later epoch than that of the left.
Roughly speaking, in both of the panels, fraction of electrons (ye) in the inner region are smaller than that in
the outer part. These results could be simply understood since the recombination process in the inner dense region
proceeds faster than in the outer envelope. We also find the opposite behavior in the inner most region of right panel,
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 7
at which ye increases as r decreases, because of the collisional ionization. It is also worth noting that Li becomes
the main provider of electrons around r ∼ 1014cm at nH,c = 1016 cm−3. Because of the combined effects of these two
(collisional ionization & presence of Li), ye never gets lower than 10
−11 as far as we consider nH,c < 10
16 cm−3 in
run-away phase.
3.1.2. Drift velocities in run-away collapse phase
In Figure 3, we show the ratio of the drift velocity vBx to the infall velocity u when nH,c = 10
16 cm−3 is satisfied.
The vertical axis shows the magnetic field strength, whereas the horizontal axis shows the distance from the cloud
center. The solid curves in Figure 3 show the contours of log10(vBx/u). The magnetic fields are dissipated in the
region where log10(vBx/u)> 0, in contrast, the fields are frozen to the gas in the region where log10(vBx/u)< 0. The
dashed curve in Figure 3 is the critical field strength Bcr that is defined by the equation
GM(r)ρ(r)
. (23)
Note that since we are interested in the collapsing gas cloud, the magnetic force needs to be weaker than the gravi-
tational force. Hence, our calculations are valid in the region where the field strength in the cloud satisfies B < Bcr.
We find clearly from Figure 3 that the frozen-in condition vBx/u < 1 is almost always satisfied as long as B is less
than the critical field strength Bcr. We also find the basically same results for other snapshots, i.e. the drift velocity
is smaller than the infall velocity anytime and anywhere if B <Bcr.
In addition, if B < Bcr is satisfied at some initial time and position (t0, r0), this condition also holds at some later
epoch (t,r(t)). This statement is proved as follows: Combining the magnetic flux conservation equation
2πr0dr0B(t0,r0) = 2πr(t)drB(t,r(t)) (24)
and the mass conservation law given in equation (4), we obtain
B(t,r(t)) =B(t0,r0)
ρ(t,r(t))r(t)
ρ(t0,r0)r0
, (25)
Using equations (23) and (25) we have
B(t,r(t))
Bcr(t,r(t))
B(t0,r0)
Bcr(t0,r0)
ρ(t,r(t))r3 (t)
ρ(t0,r0)r
The second term in right side is always less than or equals to unity in the run-away collapsing cloud, since ρ ∝ r−2
in the envelope, and ρ ∝ r−3 in the core. Thus, if B(t0,r0)/Bcr(t0,r0) < 1 is satisfied, B is less than Bcr throughout
the collapse. Remark that we assume the flux conservation as equation (25), which gives the maximal field strength.
Therefore, above inequality (B <Bcr) also holds even if the magnetic flux is dissipated.
In summary, it is concluded that the magnetic fields are always frozen to the whole cloud in the course of run-away
collapse phase, if the magnetic force is much smaller than gravitational force at the beginning of the collapse.
We also remark that the slight increase due to collisional ionization at nH,c =10
16 cm−3 implies that ye will become
larger as the central density increases up to nH,c> 10
16 cm−3, since the temperature gets higher as the density increases
for γ=1.09 polytropic gas. As a result, the magnetic field is expected to be frozen to the gas even for nH,c> 10
16 cm−3
in accretion phase.
3.2. Accretion Phase
3.2.1. Evolution of gas and chemical species in accretion phase
The evolutionary sequences of number density nH, temperature T , velocity u, and electron fraction ye in the accretion
phase are illustrated in Figure 4.
Upper left panel illustrates the evolution of density profile. The solid curves labeled as 0-4 corresponds to the epoch
at which the central accreted mass equals to 8.87× 10−2 M⊙,1.58 M⊙,12.7 M⊙,50.2 M⊙,100 M⊙, respectively. It
is clear that the density for r → 0 decreases as the collapse proceeds, which is basically the same feature found in
similarity solution of polytropic cloud in accretion phase.
Because of the high accretion velocity at later snapshots (lower left), the gas is heated efficiently by adiabatic
compression. As a result, temperature gets higher as the collapse proceeds, except at the beginning of the accretion
phase, when the gas is cooled efficiently due to the enhanced fraction of H2 molecules.
The ionization degree (electron fraction ye, lower right) also increases as the collapse proceeds (Figures.4 and 5),
due to the higher temperature for later snapshots. As a result, ye never gets lower than 10
−11. We also find that Li
is not so important as was in the final phase of run-away phase, since the hydrogen win back the position of chief
provider of electrons in the accretion phase( Figure 5).
8 Maki & Susa [Vol. ,
Fig. 4. The evolutionary sequences of the primordial gas cloud in accretion phase are shown. Four panels (a)-(d) show the spatial
distribution of following physical variables as functions of radius: (a) number density nH, (b) temperature T , (c) velocity u, and
(d) electron fraction ye. Five time sequences (0-4) are plotted. Accreted central mass Mc(t) is used as a clock. Corresponding
mass at the stages are: 0:Mc(t) = 8.87× 10
−2 M⊙, 1:Mc(t) = 1.58M⊙, 2:Mc(t) = 12.7M⊙, 3:Mc(t) = 50.2M⊙, 4:Mc(t) = 100M⊙.
3.2.2. Drift velocities in accretion phase
As briefly discussed in section 1, the ambipolar diffusion velocity is proportional to the inverse square of gas density.
Its dependence on various physical quantities are described as vambBx ∝ y−1e ρ−2B2r−1 (see equations (19)-(22)). On the
other hand, the accretion velocity scales as u ∝ (Mc/r)1/2 for r → 0. Since we consider the accretion phase starting
from the final phase of run-away collapse with γ = 1.09, the density profile and the central mass in accretion phase
also depend on γ. According to Suto & Silk (1988), we have
ρ(t,r) ∝ t(2−3γ)/2r−3/2 for r → 0 (27)
Mc(t)∝ t4−3γ (28)
Thus, the limiting behaviour of the infall velocity is described as
u(t,r)∝ t(4−3γ)/2r−1/2 for r→ 0 (29)
Combining above set of equations, we have the dependence of the ratio vambBx /u on t:
vambBx
∝ y−1e B2t(9γ−8)/2r5/2 (30)
This equation indicates that the ratio keeps growing in accretion phase for fixed r,B and ye, since 9γ− 8 is positive
for γ = 1.09. In reality, however, ye increases rapidly as the accretion proceeds (see Figures.4 and 5), and it offsets the
increase of vambBx /u.
Figures 6 and 7 illustrate the contours of log10(vBx/u) on r−B plane for two epochs. Two figures correspond to
the snapshots when the central accreted mass satisfies Mc = 12.7M⊙ and Mc = 100M⊙, respectively. The notations
are same as Figure 3. It is clear that the drift velocity is always smaller than the accretion velocity in two snapshots,
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 9
Fig. 5. The radial distributions of the fractional abundances of the main species, e, H+, H, H2, H
−, H+
, Li, and Li+
in the accretion phase. Two panels correspond to the snapshots at Mc(t) = 12.7 M⊙(left), and Mc(t) = 100 M⊙(right).
Fig. 6. The same as Fig.3, except that the contours are plotted in the accretion phase characterized at Mc =12.7M⊙.
as long as B <Bcr. In fact, we find that this is true all through the calculations. Besides, the equation (26) holds also
in the accretion phase. Considering that the density profile in the accretion phase is ρ∝ r−1.5∼−2, B is less than Bcr
throughout the collapse, if B(t0,r0)/Bcr(t0,r0) < 1, as discussed in §3.1.2. In other words, magnetic force is always
negligible if it can be ignored in the beginning of the collapse.
Thus, it is concluded that the dissipation of magnetic flux is negligible throughout the mass accretion phase, as well
as the run-away collapse phase.
4. DISCUSSION
We confirmed that the magnetic field is frozen to the star-forming primordial gas cloud even in the accretion phase.
Here we discuss the possibility of jet formation and activation of MRI considering the magnetic field strength brought
10 Maki & Susa [Vol. ,
Fig. 7. The same as Fig.6, except Mc = 100M⊙.
into the accretion disk surrounding the protostar.
The magnetic field strength brought into the accretion disk is assessed under the frozen-in condition as follows:
Bdisk =B0
rdisk
. (31)
Here B0 is the initial field strength, rdisk denotes the disk radius, whereas R describes the initial radius within which
includes the total mass of the disk-star star system M∗disk:
R= (3M∗disk/4πρ0)
1/3, (32)
where ρ0 represents the initial density of the cloud. rdisk can be evaluated by equation (15) in Tan & Mckee (2004):
rdisk ≃ 66.4AU
M∗disk
where fKep denotes the ratio of the rotation velocity to the Kepler velocity of accreting matter, which is found to be
∼ 0.5 in Abel, Bryan & Norman (2002). Thus, we have the magnetic field strength in the disk as follows:
Bdisk ≃ 7.5× 10−10G
3.7× 10−16G
103cm−3
)−2/3
M∗disk
)−40/21
Several possibilities to generate cosmological seed magnetic field have been proposed so far. Most of the mechanisms
predict BIGM <∼ 10
−19G(Widrow 2002), except the magnetic field generated by radiation transfer effects of powerful
ionizing sources such as quasars or first stars (Langer, Puget & Aghanim 2003). They suggests the possibility to
generate coherent magnetic field with BIGM ∼ 10−11G. Since the magnetic field is frozen to the primordial gas at low
densities (nH < 10
3cm−3), B0 at nH = 10
3cm−3 can be evaluated as
B0 = 3.7× 10−16G
10−19G
103cm−3
, (35)
Recently, 3-dimensional MHD simulations on primordial star formation have been performed by Machida et al.
(2006), assuming ideal MHD condition is always satisfied in the collapsing gas. In fact, our present results guarantee
this hypothesis. They found that the protostellar jet is driven in primordial environment if B0 >∼ 10
−9G at nH =
103cm−3. Comparing this condition with equation (35), it is concluded that 1) jets could be driven in first star
forming clouds if the seed field is generated by the mechanism proposed by Langer, Puget & Aghanim (2003), 2)
whereas the other mechanisms cannot generate the seed field enough to drive the jets.
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 11
MRI can be activated in the accretion disk in case the magnetic field in the disk is larger than a critical value (Tan
& Blackman 2004):
Bdisk >∼ 1.1× 10
M∗disk
10 M⊙
)1/4(
104 K
)−3/4(
)1/2(
ρdisk
5× 10−10 g cm−3
)1/2(
600 R⊙
)−3/4
, (36)
This threshold is assessed by the confrontation between the growth rate of MRI and the ohmic dissipation rate.
Combining equations (34)- (36), we find that MRI is activated in case the seed field satisfy
BIGM >∼ 1.5× 10
M∗disk
)2.155
104 K
)−3/4(
)1/2(
ρdisk
5× 10−10 g cm−3
)1/2(
600 R⊙
)−3/4
. (37)
Therefore, MRI is driven only if the seed field generation mechanism by the transfer effects of ionizing radiation works.
Based upon these arguments, we emphasize that the mechanism proposed by Langer, Puget & Aghanim (2003) should
be scrutinized since their results still based upon the argument of order-estimation.
5. SUMMARY
In this paper, we investigate the dissipation of magnetic flux in star-forming primordial gas cloud. We solve non-
equilibrium chemical reaction equations, coupled with thermal and dynamical evolution of the collapsing cloud all
through the run-away phase as well as the mass accretion phase. Thus, we obtain the detailed evolution of ionized
fraction of the gas, which enables us to assess the coupling between gas and magnetic field.
As a result, we find that the magnetic field is basically frozen to the gas anywhere in collapsing star-forming
primordial clouds at any time. Based upon this result, we find the cosmological seed magnetic field generated by most
of the mechanisms proposed so far is not sufficient to form jets as well as to activate MRI in the star-forming cloud.
Only one mechanism proposed by Langer, Puget & Aghanim (2003) is able to create sufficient field strength.
We thank Kazu Omukai for stimulating discussions. Noriaki Shibazaki and Ken Ohsuga are acknowledged for
continuous encouragement. The analysis has been made with computational facilities at Rikkyo University. This work
was supported in part by Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Grants-in-Aid,
Specially Promoted Research 16002003 and Young Scientists (B) 17740110.
Appendix 1. Density distribution of free-falling matter
In this Appendix, we derive the explicit formula for
which gives the density distribution of free-falling
matter in equation (8).
Equation of motion (5) is integrated as
ṙ2 =
−E0, E0 =
u20, (E0 > 0), (A1)
where u0 is the velocity of a fluid element at (t0,r0) where r0 denotes the position of the element at some initial time
t0. A solution of this equation is given as
r(α,α0) =
(1− cosα), (A2)
t(α,α0) = t0 −
(2E0)3/2
[(α− sinα)− (α0 − sinα0)] , (A3)
where α is the so-called development angle, and α0 is its value at t= t0. α0 also satisfies following relation
(1− cosα0). (A4)
Using above relations, (∂r/∂r0)t is derived as
M0,E0
, (A5)
where
12 Maki & Susa [Vol. ,
(1− cosα) = r
, (A6)
=−GM0
(1− cosα) =− r
, (A7)
M0,E0
sinα, (A8)
sinα0
, (A9)
= 4πr20ρ0, (A10)
, (A11)
(t− t0)
− 3(t− t0)
. (A12)
References
Abel, T., Bryan, G. L., & Norman, M. L. 2000, ApJ, 540, 39
Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295,
Bodenheimer, P., Burkert, A., Klein, R., & Boss, A. 2000,
Protostars and Planets IV 2000, eds. V. Mannings, A. P.
Boss, and S. S. Russel, U. of Arizona Press, 675
Bromm, V., Coppi, P. S., & Larson, R. B. 1999, ApJ, 527, L5
Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23
Christlieb, N. et al. 2004, ApJ, 603, 708
Colella, P., & Woodward, P.R. 1984, Jour. Comp. Phys., 54,
Flower, D. R. 2002, MNRAS, 333, 763
Frebel, A. et al. 2005, Nature, 434, 871
Galli, D., & Palla, F. 1998, A&A, 335, 403
Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595
Hollenbachm D., & McKee,C.F. 1979,ApJS, 41, 555
Iwamoto, N. et al. 2005, Science, 309, 451
Johnson, J.L. & Bromm, V. 2006, MNRAS, 366, 247
Kamaya, H., & Nishi, R. 2000, ApJ, 543, 257
Langer, M., Puget, J., & Aghanim, N. 2003, Phys. Rev. D, 67,
43505
Larson, R. 1969, MNRAS, 145, 271
Lenzuni, P., Chernoff, D. F., & Salpeter, E. E. 1991, ApJS, 76,
Lepp, S., Stancil, P. C., & Dalgarno, A. 2002, J. Phys. B, 35,
Levermore,C.D & Pomraning,G.C. 1981, ApJ, 248, L321
Machida, M. N., Omukai, K., Matsumoto, T., & Inutsuka, S.
2006, ApJL, 647, 1
Maki, H., & Susa, H. 2004, ApJ, 609, 467
Nakamura, F., & Umemura, M. 1999, ApJ, 515, 239
Nakamura, F., & Umemura, M. 2001, ApJ, 548, 19
Nakano, T., & Umebayashi, T. 1986, MNRAS, 218, 663
Omukai, K. 2000, ApJ, 534, 809
Omukai, K., & Nishi, R. 1998, ApJ, 508, 141
O’Shea, B. W., Abel, T., Whalen, D., Norman, M. L. 2005,
ApJ, 628, 5L
O’Shea, B. W. & Norman, M. L. 2006, ApJ, 648, 31
O’Shea, B. W. & Norman, M. L. 2007, ApJ, 654, 66
Penston, M. V. 1969, MNRAS, 144, 425
Saigo, K, Matsumoto, T & Umemura, M. 2005, ApJ, 615, 65
Sano, T., Inutsuka, S., & Miyama, S. M. 1998, ApJ, 506, L57
Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000,
ApJ, 543, 486
Sano, T., & Inutsuka, S. 2001, ApJ, 561, L179
Spitzer, L. 1978, Physical Processes in the Interstellar Medium,
(New York: Wiley Interscience) p.107
Stancil, P. C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1
Stone, J., Gammie, C., Balbus, S., & Hawley, J. 2000,
Protostars and Planets IV 2000, eds. V. Mannings, A. P.
Boss, and S. S. Russel, U. of Arizona Press, 589
Susa, H. 2007, ArXiv Astrophysics e-prints, arXiv:astro-
ph/0701172
Suto, Y., & Silk, J. 1988, ApJ, 326, 527
Tan, J. C., & Blackman, E. G. 2004, ApJ, 603, 401
Tan, J. C., & McKee, C. F. 2004, ApJ, 603, 383
Widrow, L. M. 2002, Rev. Mod. Phys., 74, 775
Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003,
ApJ, 592, 645
Yoshida, N., Omukai, K., Hernquist, L. & Abel, T. 2006, ApJ,
652, 6

  We investigate the dissipation of magnetic flux in primordial star-forming
clouds throughout their collapse including the run-away collapse phase as well
as the accretion phase. We solve the energy equation and the non-equilibrium
chemical reactions in the collapsing gas, in order to obtain the radial
distribution of the ionized fraction during the collapse. As a result, we find
the ionized fraction is high enough for the magnetic field to couple with the
gas throughout the evolution of the cloud. This result suggests that the jet
formation from protostars as well as the activation of magneto-rotational
instability in the accretion disk are enabled in the presence of the
cosmological seed magnetic flux proposed by Langer et al.(2003).

arXiv:0704.1853v1  [astro-ph]  14 Apr 2007
PASJ: Publ. Astron. Soc. Japan , 1–??, 〈publication date〉
c© 2019. Astronomical Society of Japan.
Dissipation of Magnetic Flux in Primordial Star Formation:
From Run-away Phase to Mass Accretion Phase
Hideki Maki
Department of Physics, Rikkyo University, Nishi-Ikebukuro, Tokyo 171-8501
hide-mk@jcom.home.ne.jp
Hajime Susa
Department of Physics, Konan University, Okamoto, Kobe 658-8501
susa@center.konan-u.ac.jp
(Received 〈reception date〉; accepted 〈acception date〉)
Abstract
We investigate the dissipation of magnetic flux in primordial star-forming clouds throughout their
collapse including the run-away collapse phase as well as the accretion phase. We solve the energy equation
and the non-equilibrium chemical reactions in the collapsing gas, in order to obtain the radial distribution
of the ionized fraction during the collapse. As a result, we find the ionized fraction is high enough for the
magnetic field to couple with the gas throughout the evolution of the cloud. This result suggests that the
jet formation from protostars as well as the activation of magneto-rotational instability in the accretion
disk are enabled in the presence of the cosmological seed magnetic flux proposed by Langer, Puget &
Aghanim (2003).
Key words: accretion, accretion disks — diffusion — early universe — stars: formation — stars:
magnetic fields
1. INTRODUCTION
Typical mass or initial mass function of population III stars are fundamental parameters that have great impacts
on subsequent structure formation of the universe. Those stars are expected to be as massive as 100−1000M⊙, which
ionize/dissociate the surrounding media. As a result, star formation activities in the neighbourhood of the stars are
highly regulated (e.g. Susa (2007) and the references are therein).
There are two chief reasons that the population III stars are expected to be very massive. First one comes from
numerical studies by several authors (Nakamura & Umemura 1999; Bromm, Coppi & Larson 1999; Abel, Bryan &
Norman 2000; Nakamura & Umemura 2001; Bromm, Coppi & Larson 2002). Those studies reveal that the prestellar
cores formed as fragments of primordial gas are as massive as ∼ 103M⊙ − 104M⊙. Second ground is the very high
accretion rates of those stars, which is as high as 10−3− 10−2M⊙yr−1. These two facts are the direct consequence of
relatively high temperature (∼ 1000K) of primordial gas because of inefficient cooling by H2.
On the other hand, abundance pattern of the hyper metal-poor stars seems to be more consistent with that of
“faint” supernovae as remnants of less massive (∼ 25M⊙) population III stars(Christlieb et al. 2004; Frebel et al.
2005; Iwamoto et al. 2005). Recent theoretical studies suggest the possibilities for the formation of such less massive
population III stars. If the primordial gas was once ionized, enhanced fraction of H2 causes more efficient cooling and
HD formation. Since HD molecules have lower excitation energy than that of H2, gas can be cooled below 100K. As
a result, the mass of the fragments could be smaller than 100M⊙(O’Shea, Abel, Whalen & Norman 2005; Johnson &
Bromm 2006). O’Shea & Norman (2007) also demonstrates that such low mass population III stars cloud be formed
directly from the cosmological density fluctuations.
Those new ideas are quite interesting and promising, however, we still do not know the actual mass of population III
stars formed from ∼ 103M⊙ − 104M⊙ prestellar cores, which are commonly found in cosmological simulations. Since
they are only 104− 106cm−3, we have to follow the subsequent evolution of the collapsing cloud. Further evolution of
prestellar cores is first investigated by Omukai & Nishi (1998), and they found that the collapse proceeds in a run-away
fashion and converges to Larson-Penston type similarity solution (Larson 1969; Penston 1969; Suto & Silk 1988) with
polytropic index γ ≃ 1.09. They also find the mass accretion rate is very large compared to the present-day forming
stars, although spherical symmetry is assumed in their radiation hydrodynamic simulations. Recently, Yoshida et al.
(2006) have performed cosmological simulations in which the run-away collapsing core is traced up to a very dense
regime (nH ≫ 1010cm−3), taking the Sobolev type line transfer approximations. They also find basically consistent
results in 3D cosmological simulations with previous 1D results by Omukai & Nishi (1998).
However, most of the mass of protostar is accumulated in the accretion phase, which has not been traced especially in
http://arxiv.org/abs/0704.1853v1
2 Maki & Susa [Vol. ,
multi-dimensional simulations. Recent numerical simulations performed by Saigo, Matsumoto & Umemura (2005) show
that disks would be formed at the center of collapsing primordial clouds, which might result in the disks surrounding
the protostars, or binaries. Thus, the actual mass of a population III star should depends on the mechanism of the
angular momentum transport in the accretion disk. There are a few possibilities of the angular momentum transport
mechanism such as gravitational torque by the nonaxisymmetric structures in the accretion disk, the interaction
among the fragments of the disk (Stone et al. 2000; Bodenheimer et al. 2000), and the turbulent viscosity triggered
by Magneto-Rotational Instability (MRI) (Hawley & Balbus 1992; Sano, Inutsuka & Miyama 1998; Sano & Inutsuka
2001). As for MRI induced turbulence, the strength of magnetic field in the disk is the key quantity to activate
the instability(Tan & Mckee 2004; Tan & Blackman 2004). It is also worth noting that recent MHD simulations by
Machida et al. (2006) suggests the possibility of bipolar jets from proto-population III stars, which also could suppress
the mass accretion onto the central core. The formation of jets also requires the presence of some level of magnetic
field. Therefore, it is quite important to assess the magnetic field strength brought into the accretion disk from initial
weak cosmological seed field (e.g. Langer, Puget & Aghanim (2003)).
Maki & Susa (2004) investigated this issue by solving detailed chemical reaction rate equations coupled with energy
equation in run-away collapsing core. They found that magnetic field is always frozen to the collapsing core in case
the strength of initial magnetic field is smaller than 10−5(nH/10
3 cm−3)0.55 G. This is comparable to the maximal
strength of magnetic field which allow the clouds to collapse. Therefore they conclude magnetic field is always frozen
to the collapsing cloud in the run-away collapsing core. They also evaluate the minimal strength of magnetic field
which is required to activate MRI assuming magnetic field is frozen to the gas not only in the run-away phase, but
also in the accretion phase.
This assumption is based upon the argument that the temperature of the accretion flow would rise faster than
that of the core, because of the shock heating. As a result, the ionization degree is expected to be higher than those
in run-away phase, which guarantee the gas to be frozen to the magnetic field lines. However, such heating is only
important in the final phase of accretion, where the flow settle onto the accretion disk. In order to investigate the
coupling of matter and magnetic field in the accretion phase, we need to calculate the actual thermal and chemical
evolution of accretion flow. Aside from the issue on ionization degree, the ambipolar diffusion velocity is proportional
to the square inverse of density (vambB ∝ ρ−2). As a result, the diffusion velocity increases rapidly in the mass accretion
phase, since the density at a fixed radius decreases as the accretion proceeds. Thus, we need a detailed treatment to
test the coupling between the gas and magnetic field in accretion phase.
In this paper, we investigate the dissipation of magnetic field in collapsing primordial gas cloud from run-away
phase to accretion phase, by solving detailed thermal and chemical rate equations. In the next section, we describe
the formulations employed. In §3, results of our calculations are shown. The formation of jets and activation of MRI
is discussed in §4. Final section is devoted to summary.
2. METHOD OF CALCULATIONS
In order to evaluate the coupling of the gas and the magnetic field, we need to assess the amount of ions and electrons
in collapsing primordial cloud. The collapse of the cloud is expected to proceed in run-away fashion (run-away phase)
in the beginning, followed by the mass accretion phase after the formation of rotationally supported disk at the center.
In this paper, we follow the chemical and thermal evolution of the materials to form star, throughout the two phases.
We assume spherically symmetric collapse of the progenitor gas, although we expect the formation of rotationally
supported disk in the very dense regime (nH ≫ 1010cm−3). It is also assumed that the magnetic field is so weak that
the dynamics of the collapsing gas is not affected by the magnetic force.
2.1. Run-away Collapse Phase
The run-away collapse phase is traced by solving 1-dimensional hydrodynamics. The equation of continuity is
= 4πr2ρ, (1)
whereas the equation of motion is
=−4πr2 ∂p
, (2)
where m is the mass within radius r, ρ is the density of the cloud, u is the velocity and p is the pressure. Since
the equation of state of collapsing primordial gas in run-away phase is known to be approximated by polytrope with
γ = 1.09 (Omukai & Nishi 1998), we employ following equation:
p=Kργ , γ = 1.09 (3)
We solve above set of equations by spherically symmetric Lagrangian hydrodynamics code developed by ourselves,
following the Piecewise Parabolic Method (PPM) described in Colella & Woodward (1984).
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 3
2.2. Mass Accretion Phase
The dynamics of accretion phase could be approximated by simple free-fall, since the flow velocity is supersonic.
Basic equations of free-falling accretion flow consists of the equation of continuity and the equation of motion:
4πρ(t0,r0)r
0dr0 = 4πρ(t,r)r
2(t)dr, (4)
d2r (t)
=−GM0
r2(t)
, (5)
where r, t and t0 represent the position of fluid element, time, and the time when the mass accretion starts, respectively.
M0 denotes the mass within r0 at t= t0:
4πr′0ρ0dr
0. (6)
The initial radius r0 is defined as r0 ≡ r(t0). Therefore, the the solution of above equation of motion is given as
r = r(t; t0,r0). The velocity of the fluid element is given by
u(t; t0,r0) =−
2E0 +2
, (7)
where E0 is the total energy defined as E0 ≡ u2(t0)/2−GM0/r0. Using equation (4), the density at (t,r) is given as
ρ(t,r) =
(∂r/∂r0)t
. (8)
(∂r/∂r0)t is the partial differentiation by r0 keeping t fixed. An explicit expression of this term is given in appendix
1. In order to clarify the validity of the free-fall approximation, it is compared with the similarity solution with
γ = 1.09 (Suto & Silk 1988) in accretion phase. The initial condition of the accretion flow is set as the final state of
the run-away collapse with γ = 1.09. Figure 1 illustrates the density and velocity distributions of several snapshots.
It is obvious that the free-fall approximation can describe the matter distribution in accretion phase very well for
polytropic gas. In other words, the thermal evolution has little effects on the dynamics of the accretion flow. Thus,
the dynamics can be approximated by free-fall formula, irrespective of the internal energy equation.
Fig. 1. Time evolution of density/velocity profiles in accretion flow are shown. Dotted lines and marks represent the result with
free-fall approximation, whereas the solid lines denote the results from the similarity solution in accretion phase with γ = 1.09.
In the accretion flow, the effective polytropic index is basically unknown. Thus, we have to solve the following
energy equation for each fluid element:
4 Maki & Susa [Vol. ,
=−p d
−L(net). (9)
In this expression, p denotes the pressure, which is related to the internal energy ε by the equation of state:
p= (γad − 1)ρε, (10)
where γad is the adiabatic exponent, L(net) denotes the net energy loss rate per unit mass.
We take into account the rovibrational line cooling by H2(Lline), the continuum radiation from the gas (Lcont), and
cooling and heating associated with chemical reactions (Ldiss, GH− , GH+
and G3body). L(net) is the sum of all these
contributions.
The cooling rate by the rovibrational transition of hydrogen molecules are assessed by the fitting formula given by
Galli & Palla (1998) in case the lines are optically thin. In the optically thick regime, we employ the escape probability
formalism, as described in equations (12)-(15) in Omukai (2000), except that we use actual velocity gradient in our
calculations.
The cooling rate owing to the continuum radiation is given by
Lcont = 4σκgasT 4 (11)
where σ denotes the Stefan-Boltzmann constant, κgas is the Planck mean opacity of gas (Lenzuni, Chernoff & Salpeter
1991), which includes bound-free absorption, free-free absorption, photodissociation, Rayleigh scattering, and collision-
induced absorption. This formula is valid in case the accreting gas is optically thin for continuum radiation. We
confirmed that the optical depth is smaller than unity throughout the accretion phase.
Cooling due to the latent heat of H2 dissociation is given by
Ldiss = 4.48
eV s−1 g−1, (12)
where nH is the number density of hydrogen nuclei, and (dyH2/dt)− is the dissociation rate of H2.
On the other hand, the gas is heated when a hydrogen molecule is formed. The heating rate per unit mass is given
GH− =
1+ncr/nH
eV s−1 g−1, (13)
1+ncr/nH
eV s−1 g−1, (14)
G3body =
3body
1+ncr/nH
eV s−1 g−1, (15)
where (dyH2/dt)H− , (dyH2/dt)H+
and (dyH2/dt)3body are the formation rates of H2 by H
− process, H+2 process and
three body reactions (Hollenbach & McKee 1979). ncr is the critical density defined as
ncr =
T 1/2
1.6yH exp
+1.4yH exp
1.2× 104
T +1200
cm−3. (16)
2.3. Chemical Reactions
Since the dissipation of the magnetic flux, as well as the temperature of the gas, strongly depends on the chemical
abundances, we have to solve the non-equilibrium chemical reaction rate equations coupled with equation of motion
and energy equation described in previous subsections 2.1 and 2.2.
The evolution of the fraction of species i is followed by solving the equations,
nHklmylym+
nHklmnylymyn,
(i= 1,2,3, · · · ,24), (17)
where yi ≡ ni/nH is the fraction of species i, klm [cm3 s−1] and klmn [cm6 s−1] are the reaction rate coefficients with
respect to two-body processes and three-body processes, respectively. In our calculation, we include 24 species: e−,
H+, H, H−, H2, H
2 , H
3 , D, D
+, D−, HD, HD+, H2D
+, He, He+, He++, HeH+, Li, Li+, Li++, Li3+, Li−, LiH,
and LiH+. We employ the latest reaction rate coefficients in the following papers, Galli & Palla (1998), Omukai
(2000), Stancil, Lepp & Dalgarno (1998), Flower (2002) and Lepp, Stancil & Dalgarno (2002). As for the radiative
recombination, we use the rate coefficients based on Spitzer (1978). These reactions are the same as our previous
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 5
paper(Maki & Susa 2004). Here we stress the importance to include above rare elements such as Li, because the
coupling of magnetic field with gas can be attained by very low fractional abundance of electrons and ions (Maki &
Susa 2004).
2.4. Drift velocity of magnetic field
Magnetic field is dissipated from star-forming gas via ohmic loss and ambipolar diffusion. We assess the drift
velocity vBx of the field lines due to these two processes, which is compared to the accretion velocity of gas. We
evaluate the drift velocity following the formulation by Nakano & Umebayashi (1986). There are two important
quantities which characterize these diffusion processes. They are τν and ων which denote the viscous damping time of
the relative velocity of charged particle ν to the neutral particles, and the cyclotron frequency of the charged particle
ν, respectively. Then, τν is expressed as
µνnnνnn〈σv〉νn
, (18)
where µνn is the reduced mass, nν , nn, and ρν are, the mean number density for the charged particle ν, the neutral
particle n, and the mass density of charged particle ν, respectively. The averaged momentum-transfer rate coefficient
for a particle ν colliding with a neutral particle is expressed by 〈σv〉νn . We use the empirical formulae for the
momentum-transfer rate coefficients (Kamaya & Nishi 2000; Sano et al. 2000).
According to Nakano & Umebayashi (1986), the drift velocity is given by
vBx =
(j ×B)x, (19)
where
A=A21 +A
2, (20)
ρντνω
1+ τ2νω
, (21)
1+ τ2νω
, (22)
B is the mean magnetic field in the primordial cloud, the suffix x means x direction component in local Cartesian
coordinates where the z direction is taken as the direction of B. We replace (1/c)(j×B)x in equation (19) with the
mean magnetic force B2/4πr, where B is the mean field strength in the cloud, r is the radius of the cloud at which
we assess the drift velocity.
2.5. Initial Conditions
We consider primordial star-forming gas clouds that formed in the mini-halos with ∼ 106 M⊙. We set uniform
and spherical gas cloud with nH = 10
3 cm−3, T0 = 250 K, and M = 10
4 M⊙. Such clouds are commonly found
3D cosmological simulations (Bromm, Coppi & Larson 1999; Abel, Bryan & Norman 2000; Abel, Bryan & Norman
2002; Bromm, Coppi & Larson 2002; Yoshida et al. 2003; O’Shea & Norman 2006). We also use the cosmological
abundance given by Galli & Palla (1998), as the initial fraction of the chemical compositions. The initial magnetic
field strength is expected to be very weak, however, its magnitude is still under discussion. Thus, we regard the field
strength as a free parameter of the calculations.
2.6. Switching from Run-away Phase to Accretion Phase
The run-away phase and the subsequent accretion phase have to be treated in different formulation from each other
as described in subsections 2.1 and 2.2. Thus, we have to switch the scheme from the method in run-away phase to
that in accretion phase. In the light of physical arguments, two phases should be switched when the central run-away
collapse is stopped due to thermal pressure or centrifugal force. In our calculation, we switch the scheme when the
sonic point of the accreting flow go inside a certain small radius. We set this radius to be rsw = 700 R⊙, which is
comparable to the disk size when the accreted mass is comparable to ∼ 1M⊙(Tan & Mckee 2004). It is also worth
noting that the central density of the core at switching is ∼ 1016 cm−3.
3. RESULTS
3.1. Runaway Collapse Phase
3.1.1. The distribution of the chemical abundances
Two snapshots of yi distribution during the run-away collapse phase are shown in Figure 2. The left panel shows
the distribution at the time when the central density nH,c reaches 10
12 cm−3, whereas the right panel represents the
6 Maki & Susa [Vol. ,
Fig. 2. The distribution of the main species, e−, H+, H, H2, H
, Li, and Li+ are plotted. The vertical axis denotes the
fractional abundance yi of the above species, and the horizontal axis denotes the radius from the core center. Left panel represents
the snapshot when the central density satisfies nH,c = 10
12 cm−3, whereas right panel shows the results at nH,c = 10
16 cm−3.
Fig. 3. Drift velocity vBx as a function of the radius r and the field strength B at the central density . Contour map of
vBx/u is shown on r − B plane for nH,c = 10
16 cm−3 Solid curves represent the constant loci along which log(vBx/u) equal
to the values labeled on the curve. The dashed curve represents the critical field strength Bcr given by the equation (23).
distribution at nH,c = 10
16 cm−3. The horizontal axis denotes the radius from the center of the primordial gas cloud
r [cm], and the vertical axis is the fractional abundances yi of each species. Since the central density increases as the
collapse proceeds, right panel represents the later epoch than that of the left.
Roughly speaking, in both of the panels, fraction of electrons (ye) in the inner region are smaller than that in
the outer part. These results could be simply understood since the recombination process in the inner dense region
proceeds faster than in the outer envelope. We also find the opposite behavior in the inner most region of right panel,
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 7
at which ye increases as r decreases, because of the collisional ionization. It is also worth noting that Li becomes
the main provider of electrons around r ∼ 1014cm at nH,c = 1016 cm−3. Because of the combined effects of these two
(collisional ionization & presence of Li), ye never gets lower than 10
−11 as far as we consider nH,c < 10
16 cm−3 in
run-away phase.
3.1.2. Drift velocities in run-away collapse phase
In Figure 3, we show the ratio of the drift velocity vBx to the infall velocity u when nH,c = 10
16 cm−3 is satisfied.
The vertical axis shows the magnetic field strength, whereas the horizontal axis shows the distance from the cloud
center. The solid curves in Figure 3 show the contours of log10(vBx/u). The magnetic fields are dissipated in the
region where log10(vBx/u)> 0, in contrast, the fields are frozen to the gas in the region where log10(vBx/u)< 0. The
dashed curve in Figure 3 is the critical field strength Bcr that is defined by the equation
GM(r)ρ(r)
. (23)
Note that since we are interested in the collapsing gas cloud, the magnetic force needs to be weaker than the gravi-
tational force. Hence, our calculations are valid in the region where the field strength in the cloud satisfies B < Bcr.
We find clearly from Figure 3 that the frozen-in condition vBx/u < 1 is almost always satisfied as long as B is less
than the critical field strength Bcr. We also find the basically same results for other snapshots, i.e. the drift velocity
is smaller than the infall velocity anytime and anywhere if B <Bcr.
In addition, if B < Bcr is satisfied at some initial time and position (t0, r0), this condition also holds at some later
epoch (t,r(t)). This statement is proved as follows: Combining the magnetic flux conservation equation
2πr0dr0B(t0,r0) = 2πr(t)drB(t,r(t)) (24)
and the mass conservation law given in equation (4), we obtain
B(t,r(t)) =B(t0,r0)
ρ(t,r(t))r(t)
ρ(t0,r0)r0
, (25)
Using equations (23) and (25) we have
B(t,r(t))
Bcr(t,r(t))
B(t0,r0)
Bcr(t0,r0)
ρ(t,r(t))r3 (t)
ρ(t0,r0)r
The second term in right side is always less than or equals to unity in the run-away collapsing cloud, since ρ ∝ r−2
in the envelope, and ρ ∝ r−3 in the core. Thus, if B(t0,r0)/Bcr(t0,r0) < 1 is satisfied, B is less than Bcr throughout
the collapse. Remark that we assume the flux conservation as equation (25), which gives the maximal field strength.
Therefore, above inequality (B <Bcr) also holds even if the magnetic flux is dissipated.
In summary, it is concluded that the magnetic fields are always frozen to the whole cloud in the course of run-away
collapse phase, if the magnetic force is much smaller than gravitational force at the beginning of the collapse.
We also remark that the slight increase due to collisional ionization at nH,c =10
16 cm−3 implies that ye will become
larger as the central density increases up to nH,c> 10
16 cm−3, since the temperature gets higher as the density increases
for γ=1.09 polytropic gas. As a result, the magnetic field is expected to be frozen to the gas even for nH,c> 10
16 cm−3
in accretion phase.
3.2. Accretion Phase
3.2.1. Evolution of gas and chemical species in accretion phase
The evolutionary sequences of number density nH, temperature T , velocity u, and electron fraction ye in the accretion
phase are illustrated in Figure 4.
Upper left panel illustrates the evolution of density profile. The solid curves labeled as 0-4 corresponds to the epoch
at which the central accreted mass equals to 8.87× 10−2 M⊙,1.58 M⊙,12.7 M⊙,50.2 M⊙,100 M⊙, respectively. It
is clear that the density for r → 0 decreases as the collapse proceeds, which is basically the same feature found in
similarity solution of polytropic cloud in accretion phase.
Because of the high accretion velocity at later snapshots (lower left), the gas is heated efficiently by adiabatic
compression. As a result, temperature gets higher as the collapse proceeds, except at the beginning of the accretion
phase, when the gas is cooled efficiently due to the enhanced fraction of H2 molecules.
The ionization degree (electron fraction ye, lower right) also increases as the collapse proceeds (Figures.4 and 5),
due to the higher temperature for later snapshots. As a result, ye never gets lower than 10
−11. We also find that Li
is not so important as was in the final phase of run-away phase, since the hydrogen win back the position of chief
provider of electrons in the accretion phase( Figure 5).
8 Maki & Susa [Vol. ,
Fig. 4. The evolutionary sequences of the primordial gas cloud in accretion phase are shown. Four panels (a)-(d) show the spatial
distribution of following physical variables as functions of radius: (a) number density nH, (b) temperature T , (c) velocity u, and
(d) electron fraction ye. Five time sequences (0-4) are plotted. Accreted central mass Mc(t) is used as a clock. Corresponding
mass at the stages are: 0:Mc(t) = 8.87× 10
−2 M⊙, 1:Mc(t) = 1.58M⊙, 2:Mc(t) = 12.7M⊙, 3:Mc(t) = 50.2M⊙, 4:Mc(t) = 100M⊙.
3.2.2. Drift velocities in accretion phase
As briefly discussed in section 1, the ambipolar diffusion velocity is proportional to the inverse square of gas density.
Its dependence on various physical quantities are described as vambBx ∝ y−1e ρ−2B2r−1 (see equations (19)-(22)). On the
other hand, the accretion velocity scales as u ∝ (Mc/r)1/2 for r → 0. Since we consider the accretion phase starting
from the final phase of run-away collapse with γ = 1.09, the density profile and the central mass in accretion phase
also depend on γ. According to Suto & Silk (1988), we have
ρ(t,r) ∝ t(2−3γ)/2r−3/2 for r → 0 (27)
Mc(t)∝ t4−3γ (28)
Thus, the limiting behaviour of the infall velocity is described as
u(t,r)∝ t(4−3γ)/2r−1/2 for r→ 0 (29)
Combining above set of equations, we have the dependence of the ratio vambBx /u on t:
vambBx
∝ y−1e B2t(9γ−8)/2r5/2 (30)
This equation indicates that the ratio keeps growing in accretion phase for fixed r,B and ye, since 9γ− 8 is positive
for γ = 1.09. In reality, however, ye increases rapidly as the accretion proceeds (see Figures.4 and 5), and it offsets the
increase of vambBx /u.
Figures 6 and 7 illustrate the contours of log10(vBx/u) on r−B plane for two epochs. Two figures correspond to
the snapshots when the central accreted mass satisfies Mc = 12.7M⊙ and Mc = 100M⊙, respectively. The notations
are same as Figure 3. It is clear that the drift velocity is always smaller than the accretion velocity in two snapshots,
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 9
Fig. 5. The radial distributions of the fractional abundances of the main species, e, H+, H, H2, H
−, H+
, Li, and Li+
in the accretion phase. Two panels correspond to the snapshots at Mc(t) = 12.7 M⊙(left), and Mc(t) = 100 M⊙(right).
Fig. 6. The same as Fig.3, except that the contours are plotted in the accretion phase characterized at Mc =12.7M⊙.
as long as B <Bcr. In fact, we find that this is true all through the calculations. Besides, the equation (26) holds also
in the accretion phase. Considering that the density profile in the accretion phase is ρ∝ r−1.5∼−2, B is less than Bcr
throughout the collapse, if B(t0,r0)/Bcr(t0,r0) < 1, as discussed in §3.1.2. In other words, magnetic force is always
negligible if it can be ignored in the beginning of the collapse.
Thus, it is concluded that the dissipation of magnetic flux is negligible throughout the mass accretion phase, as well
as the run-away collapse phase.
4. DISCUSSION
We confirmed that the magnetic field is frozen to the star-forming primordial gas cloud even in the accretion phase.
Here we discuss the possibility of jet formation and activation of MRI considering the magnetic field strength brought
10 Maki & Susa [Vol. ,
Fig. 7. The same as Fig.6, except Mc = 100M⊙.
into the accretion disk surrounding the protostar.
The magnetic field strength brought into the accretion disk is assessed under the frozen-in condition as follows:
Bdisk =B0
rdisk
. (31)
Here B0 is the initial field strength, rdisk denotes the disk radius, whereas R describes the initial radius within which
includes the total mass of the disk-star star system M∗disk:
R= (3M∗disk/4πρ0)
1/3, (32)
where ρ0 represents the initial density of the cloud. rdisk can be evaluated by equation (15) in Tan & Mckee (2004):
rdisk ≃ 66.4AU
M∗disk
where fKep denotes the ratio of the rotation velocity to the Kepler velocity of accreting matter, which is found to be
∼ 0.5 in Abel, Bryan & Norman (2002). Thus, we have the magnetic field strength in the disk as follows:
Bdisk ≃ 7.5× 10−10G
3.7× 10−16G
103cm−3
)−2/3
M∗disk
)−40/21
Several possibilities to generate cosmological seed magnetic field have been proposed so far. Most of the mechanisms
predict BIGM <∼ 10
−19G(Widrow 2002), except the magnetic field generated by radiation transfer effects of powerful
ionizing sources such as quasars or first stars (Langer, Puget & Aghanim 2003). They suggests the possibility to
generate coherent magnetic field with BIGM ∼ 10−11G. Since the magnetic field is frozen to the primordial gas at low
densities (nH < 10
3cm−3), B0 at nH = 10
3cm−3 can be evaluated as
B0 = 3.7× 10−16G
10−19G
103cm−3
, (35)
Recently, 3-dimensional MHD simulations on primordial star formation have been performed by Machida et al.
(2006), assuming ideal MHD condition is always satisfied in the collapsing gas. In fact, our present results guarantee
this hypothesis. They found that the protostellar jet is driven in primordial environment if B0 >∼ 10
−9G at nH =
103cm−3. Comparing this condition with equation (35), it is concluded that 1) jets could be driven in first star
forming clouds if the seed field is generated by the mechanism proposed by Langer, Puget & Aghanim (2003), 2)
whereas the other mechanisms cannot generate the seed field enough to drive the jets.
No. ] Dissipation of Magnetic Flux in Primordial Star Formation 11
MRI can be activated in the accretion disk in case the magnetic field in the disk is larger than a critical value (Tan
& Blackman 2004):
Bdisk >∼ 1.1× 10
M∗disk
10 M⊙
)1/4(
104 K
)−3/4(
)1/2(
ρdisk
5× 10−10 g cm−3
)1/2(
600 R⊙
)−3/4
, (36)
This threshold is assessed by the confrontation between the growth rate of MRI and the ohmic dissipation rate.
Combining equations (34)- (36), we find that MRI is activated in case the seed field satisfy
BIGM >∼ 1.5× 10
M∗disk
)2.155
104 K
)−3/4(
)1/2(
ρdisk
5× 10−10 g cm−3
)1/2(
600 R⊙
)−3/4
. (37)
Therefore, MRI is driven only if the seed field generation mechanism by the transfer effects of ionizing radiation works.
Based upon these arguments, we emphasize that the mechanism proposed by Langer, Puget & Aghanim (2003) should
be scrutinized since their results still based upon the argument of order-estimation.
5. SUMMARY
In this paper, we investigate the dissipation of magnetic flux in star-forming primordial gas cloud. We solve non-
equilibrium chemical reaction equations, coupled with thermal and dynamical evolution of the collapsing cloud all
through the run-away phase as well as the mass accretion phase. Thus, we obtain the detailed evolution of ionized
fraction of the gas, which enables us to assess the coupling between gas and magnetic field.
As a result, we find that the magnetic field is basically frozen to the gas anywhere in collapsing star-forming
primordial clouds at any time. Based upon this result, we find the cosmological seed magnetic field generated by most
of the mechanisms proposed so far is not sufficient to form jets as well as to activate MRI in the star-forming cloud.
Only one mechanism proposed by Langer, Puget & Aghanim (2003) is able to create sufficient field strength.
We thank Kazu Omukai for stimulating discussions. Noriaki Shibazaki and Ken Ohsuga are acknowledged for
continuous encouragement. The analysis has been made with computational facilities at Rikkyo University. This work
was supported in part by Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Grants-in-Aid,
Specially Promoted Research 16002003 and Young Scientists (B) 17740110.
Appendix 1. Density distribution of free-falling matter
In this Appendix, we derive the explicit formula for
which gives the density distribution of free-falling
matter in equation (8).
Equation of motion (5) is integrated as
ṙ2 =
−E0, E0 =
u20, (E0 > 0), (A1)
where u0 is the velocity of a fluid element at (t0,r0) where r0 denotes the position of the element at some initial time
t0. A solution of this equation is given as
r(α,α0) =
(1− cosα), (A2)
t(α,α0) = t0 −
(2E0)3/2
[(α− sinα)− (α0 − sinα0)] , (A3)
where α is the so-called development angle, and α0 is its value at t= t0. α0 also satisfies following relation
(1− cosα0). (A4)
Using above relations, (∂r/∂r0)t is derived as
M0,E0
, (A5)
where
12 Maki & Susa [Vol. ,
(1− cosα) = r
, (A6)
=−GM0
(1− cosα) =− r
, (A7)
M0,E0
sinα, (A8)
sinα0
, (A9)
= 4πr20ρ0, (A10)
, (A11)
(t− t0)
− 3(t− t0)
. (A12)
References
Abel, T., Bryan, G. L., & Norman, M. L. 2000, ApJ, 540, 39
Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295,
Bodenheimer, P., Burkert, A., Klein, R., & Boss, A. 2000,
Protostars and Planets IV 2000, eds. V. Mannings, A. P.
Boss, and S. S. Russel, U. of Arizona Press, 675
Bromm, V., Coppi, P. S., & Larson, R. B. 1999, ApJ, 527, L5
Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23
Christlieb, N. et al. 2004, ApJ, 603, 708
Colella, P., & Woodward, P.R. 1984, Jour. Comp. Phys., 54,
Flower, D. R. 2002, MNRAS, 333, 763
Frebel, A. et al. 2005, Nature, 434, 871
Galli, D., & Palla, F. 1998, A&A, 335, 403
Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595
Hollenbachm D., & McKee,C.F. 1979,ApJS, 41, 555
Iwamoto, N. et al. 2005, Science, 309, 451
Johnson, J.L. & Bromm, V. 2006, MNRAS, 366, 247
Kamaya, H., & Nishi, R. 2000, ApJ, 543, 257
Langer, M., Puget, J., & Aghanim, N. 2003, Phys. Rev. D, 67,
43505
Larson, R. 1969, MNRAS, 145, 271
Lenzuni, P., Chernoff, D. F., & Salpeter, E. E. 1991, ApJS, 76,
Lepp, S., Stancil, P. C., & Dalgarno, A. 2002, J. Phys. B, 35,
Levermore,C.D & Pomraning,G.C. 1981, ApJ, 248, L321
Machida, M. N., Omukai, K., Matsumoto, T., & Inutsuka, S.
2006, ApJL, 647, 1
Maki, H., & Susa, H. 2004, ApJ, 609, 467
Nakamura, F., & Umemura, M. 1999, ApJ, 515, 239
Nakamura, F., & Umemura, M. 2001, ApJ, 548, 19
Nakano, T., & Umebayashi, T. 1986, MNRAS, 218, 663
Omukai, K. 2000, ApJ, 534, 809
Omukai, K., & Nishi, R. 1998, ApJ, 508, 141
O’Shea, B. W., Abel, T., Whalen, D., Norman, M. L. 2005,
ApJ, 628, 5L
O’Shea, B. W. & Norman, M. L. 2006, ApJ, 648, 31
O’Shea, B. W. & Norman, M. L. 2007, ApJ, 654, 66
Penston, M. V. 1969, MNRAS, 144, 425
Saigo, K, Matsumoto, T & Umemura, M. 2005, ApJ, 615, 65
Sano, T., Inutsuka, S., & Miyama, S. M. 1998, ApJ, 506, L57
Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000,
ApJ, 543, 486
Sano, T., & Inutsuka, S. 2001, ApJ, 561, L179
Spitzer, L. 1978, Physical Processes in the Interstellar Medium,
(New York: Wiley Interscience) p.107
Stancil, P. C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1
Stone, J., Gammie, C., Balbus, S., & Hawley, J. 2000,
Protostars and Planets IV 2000, eds. V. Mannings, A. P.
Boss, and S. S. Russel, U. of Arizona Press, 589
Susa, H. 2007, ArXiv Astrophysics e-prints, arXiv:astro-
ph/0701172
Suto, Y., & Silk, J. 1988, ApJ, 326, 527
Tan, J. C., & Blackman, E. G. 2004, ApJ, 603, 401
Tan, J. C., & McKee, C. F. 2004, ApJ, 603, 383
Widrow, L. M. 2002, Rev. Mod. Phys., 74, 775
Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003,
ApJ, 592, 645
Yoshida, N., Omukai, K., Hernquist, L. & Abel, T. 2006, ApJ,
652, 6
