# GPU-Net: Lightweight U-Net with More Diverse Features

###### Abstract

Image segmentation is an important task in the medical image field and many
convolutional neural networks (CNNs) based methods have been proposed, among
which U-Net and its variants show promising performance. In this paper, we
propose GP-module and GPU-Net based on U-Net, which can learn more diverse
features by introducing Ghost module and atrous spatial pyramid pooling
(ASPP). Our method achieves better performance with more than $4\times$ fewer
parameters and $2\times$ fewer FLOPs, which provides a new potential direction
for future research. Our plug-and-play module can also be applied to existing
segmentation methods to further improve their performance.

Index Terms—  segmentation, U-Net, ghost module, atrous spatial pyramid
pooling

##  1 Introduction

Recently, many impressive CNN-based image segmentation models have been
proposed. However, it is hard to apply these natural image segmentation models
into medical image segmentation directly since there exist domain gaps. Aiming
at this, researchers proposed U-Net [1] and U-Net based variants, which
achieve remarkable performance. These variants mainly focus on optimization of
network structure, which introduce more parameters and generate useful feature
maps inefficiently. In this paper, we explore the enhancement that well-
learned features can bring and boost U-Net with more higher quality diverse
features and fewer parameters. Our main contribution is proposing a
lightweight version of U-Net which can achieve competitive and even better
performance while significantly reduces the number of parameters and FLOPS. We
name our method GPU-Net and experiments show GPU-Net can achieve state-of-the-
art performance. To the best of our knowledge, it is the first paper that
explore the possibility of applying ghost-module and its variants into U-Net.
We believe our method can be an interest topic to discuss for the medical
image segmentation community.

##  2 Related Work

Image segmentation requires classifying each image pixel on an individual
basis. While the deep learning research in natural image segmentation is in
full swing, U-Net [1] is proposed specifically for medical image segmentation.
U-Net can be trained with a relatively small number of medical data from
scratch and achieve competitive performance. Based on U-Net, a series of
methods have been proposed. R2U-Net [2] and Attention U-Net [3] apply
recurrent module and attention mechanism into U-Net, respectively. These
variants all develop the potential of U-Net from different angles but ignore
the problem of feature redundancy. Redundancy in feature maps can be important
but it is better to achieve this kind of redundancy using more efficient ways
like getting a set of intrinsic feature maps first and then generating many
ghost feature maps based on them as proposed in GhostNet [4]. The intrinsic
feature maps can have no redundancy so as to ensure efficiency. The increase
in parameters also impair the innovation of architectures in their methods. In
this paper, our target is to explore how U-Net can benefit from sufficient
diverse features. To get more high quality features with simple operations in
network, we borrow and improve the ghost module in GhostNet [4] and propose
GP-module, a lightweight module that can boost U-Net with more valuable
features and fewer parameters.

(a) Image

(b) First Convolution Feature

(c) Last Convolution Feature

Fig. 1: Visualization of Feature Redundancy.

##  3 Methodology

###  3.1 GP-module

U-Net and its existing variants have feature redundancy problem as shown in
Fig 1, where we visualize the first and last convolution layer of U-Net
segmentation results. This redundancy includes two aspects. On the one hand,
some feature maps contain little or even no useful (green parts) information.
On the other hand, there exist similar feature maps (as red, pink and cyan
parts). Redundancy in feature maps can benefit model performance to some
extent. However, the existing methods obtain sufficient feature maps at the
expense of efficiency. To solve the problems, we propose GP-module, an
enhanced version of ghost module that is first introduced in GhostNet [4].
Ghost module uses cheap linear operations base on a handful of intrinsic
feature maps to generate comparable features with traditional convolution
operation as shown in Fig 2. To be specific, ordinary convolution operation
can be formulated as $Y=X\ast f$, where $X\in\mathbf{R}^{c\times w\times h}$
is the input feature map with $c$ channels and $h$ and $w$ are the height and
width, $\ast$ is the convolution operation, $f\in\mathbf{R}^{c\times k\times
k\times n}$ is the convolution filters to produce $n$ feature maps and
$k\times k$ is the kernel size, $Y\in\mathbf{R}^{w^{\prime}\times
h^{\prime}\times n}$ is the output feature map with $n$ channels and
$h^{\prime}$ and $w^{\prime}$ are the height and width. In a convolution
operation, the number of parameters and FLOPs required can be calculated as
Eq.1 and Eq.2, respectively. Both of them can be very large when $c$ and $n$
are very large, which is the usual case.

$N^{Para}_{conv}=c\cdot k\cdot k\cdot n$ (1)

$N^{FLOPs}_{conv}=c\cdot k\cdot k\cdot w^{\prime}\cdot h^{\prime}\cdot n$ (2)

We can find out from Fig 1 that there are many similar features and they can
be generate using fewer parameters and FLOPs. Ghost module [4] deal with this
problem by dividing feature maps into two parts. One part is a small number of
intrinsic feature maps and the other part of feature maps is ”ghosts” of the
intrinsic feature maps. They are called ghost feature maps and are produced by
using some cheap transformations based on the intrinsic feature maps. The idea
is generating $m$ intrinsic feature maps and for each intrinsic feature map,
applying several cheap linear operations to get $s$ ghost feature maps. By
this way we obtain the $n=m\cdot s$ feature maps desired. The way to generate
intrinsic feature maps is the same as ordinary convolution operation and the
hyper-parameters (kernel size, padding, stride, etc.) are consistent with
ordinary convolution to keep the same output spatial size ($Y^{\prime}=X\ast
f^{\prime}$). The only difference is that $f^{\prime}\in\mathbf{R}^{c\times
k\times k\times m}$ and $Y^{\prime}\in\mathbf{R}^{w^{\prime}\times
h^{\prime}\times m}$, where $m\leq n$ so the number of parameters can be
greatly reduced. After getting intrinsic feature maps, ghost features can be
produced by applying a series of cheap operations on each intrinsic feature as
Eq.3:

$y_{ij}=G_{ij}(y^{\prime}_{i}),\quad\forall i=1,...,m,j=1,...,s$ (3)

where $y^{\prime}_{i}$ is the $i$-th intrinsic feature map of $Y^{\prime}$ and
$G_{ij}$ is the $j$-th linear operation (e.g. $3\times 3$ and $5\times 5$
linear kernels) applied on $y^{\prime}_{i}$ to generate the $j$-th ghost
feature map $y_{ij}$. Each $y^{\prime}_{i}$ can get $s$ ghost feature maps
except that the last one is the identity mapping for preserving the intrinsic
feature map for each $y^{\prime}_{i}$. So we can obtain $n=m\cdot s$ output
feature maps $Y=[y_{11},y_{12},...,y_{1s},...y_{ms}]$. Note that the input and
output of $G_{ij}$ are all single-channel feature maps, Eq.3 can be easily
implemented by using depth-wise convolution. Ghost module can produce
sufficient ghost features by selecting $s$ and kernel size in $G_{ij}$. But it
uses the same single size kernel in all the $s$ linear operations for each
intrinsic feature map $y^{\prime}_{i}$, which is not conducive to generating
diverse and informative feature maps. In order to better handle this problem,
we introduce atrous spatial pyramid pooling (ASPP) mechanism proposed in [5]
and get our GP-module as shown in Fig 2. Our GP-module can be expressed as
Eq.4:

$y_{ij}=G^{\prime}_{i}(y^{\prime}_{i},j),\quad\forall i=1,...,m,j=1,...,s$ (4)

$G^{\prime}_{i}$ is explicitly related to $j$, which means for each intrinsic
feature map $y^{\prime}_{i}$, when generating more than $s>2$ ghost feature
maps, GP-module applies the same kernel size but different dilation rates
(except the identity mapping). By introducing ASPP, GP-module can expand the
receptive field and capture multi-scale contextual information so as to
generate heterogeneous representative ghost feature maps without additional
parameters and FLOPs given that GP-module uses the same kernel size as ghost
module. For each ghost module/GP-module, there is one ordinary convolution
operation to generate $m$ intrinsic feature maps and
$m\cdot(s-1)=\frac{n}{s}\cdot(s-1)$ linear operations (one of the $s$ is the
identity mapping operation) to get ghost feature maps. Note the kernel size in
linear operations is $d\times d$ and the number of parameters and FLOPs
required in ghost module/GP-module can be calculated as Eq.5 and Eq.6,
respectively. In practice, $d\times d$ and $k\times k$ have the similar
magnitude and $s$ can be much smaller than $c$. So the parameters compression
ratio and FLOPs acceleration ratio can be calculated as Eq.7 and Eq.8,
respectively.

$N^{Para}_{gp}=c\cdot k\cdot k\cdot\frac{n}{s}+\frac{n}{s}\cdot(s-1)\cdot
d\cdot d$ (5)

$N^{FLOPs}_{gp}=c\cdot k\cdot k\cdot\frac{n}{s}\cdot w^{\prime}\cdot
h^{\prime}+\frac{n}{s}\cdot(s-1)\cdot d\cdot d\cdot w^{\prime}\cdot
h^{\prime}$ (6)

$r_{P}=\frac{N^{Para}_{conv}}{N^{Para}_{gp}}\approx\frac{c\cdot s\cdot k\cdot
k}{c\cdot k\cdot k+(s-1)\cdot d\cdot d}\approx\frac{c\cdot s}{c+s-1}\approx s$
(7)

$r_{F}=\frac{N^{FLOPs}_{conv}}{N^{FLOPs}_{gp}}\approx\frac{c\cdot
s}{c+s-1}\approx s$ (8)

The parameters compression ratio is equal to FLOPs acceleration ratio which
means our GP-module does not cost extra computing resources based on ghost
module and greatly reduced parameters and FLOPs compared to ordinary
convolution operation.

(a) Ordinary Convolution

(b) Ghost Module

(c) GP-module

Fig. 2: Visualization of Feature Redundancy.

###  3.2 Network Architecture

Based on GP-module, we introduce GP-bottleneck (GP-bneck) as shown in Fig 3,
which is similar as G-bneck in [4]. GP-bottleneck has a residual structure
proposed in ResNet [6] and two stacked GP-modules, which corresponds to the
two continuous convolution operations of each level in U-Net. Batch
normalization (BN) [7] is used after each GP-module and ReLU is used only
after the first GP-module in GP-bottleneck. After getting GP-bottleneck, we
replace the convolution operations in U-Net with our GP-bottleneck and name it
GPU-Net, which is a lightweight and powerful (like its name) network for
medical image segmentation as shown in Fig 3, where the red module indicates
the location of the replacement. We use binary cross entropy as loss function
to train the whole network.

(a) GP-bneck

(b) U-Net & Ghost U-Net & GPU-Net

Fig. 3: Visualization of Feature Redundancy.

##  4 Experiments and Results

###  4.1 Dataset

Skin Lesion Segmentation Skin cancer is the most common cancer and accurate
predictions of lesion segmentation boundaries in dermoscopic images can
benefit clinical treatment. The dataset is from MICCAI 2018 Workshop [8] (ISIC
for short). ISIC dataset contains 2594 images in total and is split into
training set (1814 samples, $\sim 70\%$ ), validation set (260 samples, $\sim
10\%$), and test set (520 samples, $\sim 20\%$). All the samples are resized
into $192\times 256$ given that the original samples are slightly different.

Lung Segmentation Lung segmentation is important for analyzing lung related
diseases, and can be further applied into lung lesion segmentation and other
problems. This dataset is from Kaggle Challenge (LUNA for short). LUNA dataset
contains 267 2D CT images and is split into training set (186 samples, $\sim
70\%$ ), validation set (27 samples, $\sim 10\%$), and test set (54 samples,
$\sim 20\%$). The original image size is $512\times 512$ and we resize all the
images into $256\times 256$ in the experiment.

Nuclei Segmentation Identifying the cells’ nuclei is the starting point of
many research and can help researchers better understand the underlying
biological processes. This dataset comes from Kaggle 2018 Data Science Bowl
(DSB for short). DSB dataset has 670 nucleus images and corresponding masks.
The whole dataset is split into training set (469 samples, $70\%$ ),
validation set (67 samples, $10\%$), and test set (134 samples, $20\%$). The
original image size is $96\times 96$ and remain the same in the experiment.

$AC=\frac{TP+TN}{TP+TN+FP+FN}$ (9)

$F1=2\frac{\left|GT\cap{SR}\right|}{\left|GT\right|+\left|SR\right|},JS=\frac{\left|GT\cap{SR}\right|}{\left|GT\cup{SR}\right|}$
(10)

Table 1: Experimental Results on Three Datasets Dateset | Methods | AC | F1 | JS | Params(M) | FLOPs(G)
---|---|---|---|---|---|---
ISIC | U-Net | 0.9554 | 0.8791 | 0.8071 | 34.53 | 49.10
| Ghost U-Net | 0.9605 | 0.8852 | 0.8140 | 9.31 | 18.78
| GPU-Net | 0.9613 | 0.8926 | 0.8232 | 8.27 | 17.56
LUNA | U-Net | 0.9813 | 0.9785 | 0.9644 | 34.53 | 65.47
| Ghost U-Net | 0.9879 | 0.9796 | 0.9665 | 9.31 | 25.05
| GPU-Net | 0.9892 | 0.9811 | 0.9675 | 8.27 | 23.42
DSB | U-Net | 0.9719 | 0.8698 | 0.7944 | 34.53 | 9.21
| Ghost U-Net | 0.9716 | 0.8731 | 0.7954 | 9.31 | 3.52
| GPU-Net | 0.9727 | 0.8819 | 0.8049 | 8.27 | 3.29

###  4.2 Experimental Setup

We implement all the experiments on a NVIDIA TITAN V GPU. We use a batch size
of 4 for ISIC dataset, 2 for LUNA dataset and 16 for DSB dataset since they
have different input sizes. We train all the models for 100 epochs and the
initial learning rate is $0.001$. We compare our GPU-Net with original U-Net
and Ghost U-Net (replace convolution operations in U-Net with G-bneck in [4]).
We set $s=2$ and $d=3$ in ghost module of G-bneck, which can achieve best
performance according to [4]. For our GP-module, we set $s=6$, which means
$s-1=5$ cheap operations. We set the same $d=3$ for 4 of the 5 cheap
operations and their dilation rates are 1, 6, 12, 18, respectively. The last
cheap operation is $1\times 1$ depth-wise convolution. The other hyper-
parameters in GPU-Net/Ghost U-Net are consistent with U-Net. Our code is
available at https://github.com/Heng14/GPU-Net.

(a) ISIC

(b) LUNA

(c) DSB

Fig. 4: Segmentation result samples. From left to right: input images, ground
truth, U-Net, Ghost U-Net and GPU-Net outputs.

###  4.3 Results and Discussion

To make a detailed comparison of the model performance, we take Accuracy (AC,
Eq.9), F1-score and Jaccard similarity (F1 and JS, Eq.10) as quantitative
analysis metrics. Variables involved in these formulas are: True Positive
(TP), True Negative (TN), False Positive (FP), False Negative (FN), Ground
Truth(GT) and Segmentation Result (SR).

Experimental results are listed as Table 1 and our method can achieve better
segmentation performance with significantly reduced parameters and FLOPs. We
also show some segmentation results in Fig 4, where we can find that our GPU-
Net can better capture details and have better segmentation results at the
edges (red circles). To better demonstrate the important role of our GP-
module, we visualize the feature maps in the first level of the three networks
as in Fig 5. We can see that U-Net has many similar and blank feature maps
while Ghost U-Net and our GPU-Net can make full use of all feature maps. Our
GPU-Net can learn more diverse feature maps since ASPP module guarantees the
change of the receptive field. The feature maps in the front part tend to
learn textures, and the feature maps in the back tend to learn edges. In this
way our method can achieve better performance.

##  5 Conclusion

In this paper, we propose a lightweight model named GPU-Net, which can learn
diverse feature maps and get better segmentation performance. We test our
method on three dataset and visualize several feature maps to proof the
effectiveness and efficiency of our method. Our plug-and-play GP-module can
also be applied to existing segmentation methods to further improve their
performance with fewer parameters and fewer FLOPs, which shed the light on the
further research.

(a) U-Net

(b) Ghost U-Net

(c) GPU-Net

Fig. 5: Visualization of Feature Maps.

## References

  * [1] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in International Conference on Medical image computing and computer-assisted intervention. Springer, 2015, pp. 234–241.
  * [2] Md Zahangir Alom, Mahmudul Hasan, Chris Yakopcic, Tarek M Taha, and Vijayan K Asari,  “Recurrent residual convolutional neural network based on u-net (r2u-net) for medical image segmentation,”  arXiv preprint arXiv:1802.06955, 2018.
  * [3] Ozan Oktay, Jo Schlemper, Loic Le Folgoc, Matthew Lee, Mattias Heinrich, Kazunari Misawa, Kensaku Mori, Steven McDonagh, Nils Y Hammerla, Bernhard Kainz, et al.,  “Attention u-net: Learning where to look for the pancreas,”  arXiv preprint arXiv:1804.03999, 2018.
  * [4] Kai Han, Yunhe Wang, Qi Tian, Jianyuan Guo, Chunjing Xu, and Chang Xu,  “Ghostnet: More features from cheap operations,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 1580–1589.
  * [5] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille,  “Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs,”  IEEE transactions on pattern analysis and machine intelligence, vol. 40, no. 4, pp. 834–848, 2017.
  * [6] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [7] Sergey Ioffe and Christian Szegedy,  “Batch normalization: Accelerating deep network training by reducing internal covariate shift,”  in International conference on machine learning. PMLR, 2015, pp. 448–456.
  * [8] Noel Codella, Veronica Rotemberg, Philipp Tschandl, M Emre Celebi, Stephen Dusza, David Gutman, Brian Helba, Aadi Kalloo, Konstantinos Liopyris, Michael Marchetti, et al.,  “Skin lesion analysis toward melanoma detection 2018: A challenge hosted by the international skin imaging collaboration (isic),”  arXiv preprint arXiv:1902.03368, 2019.


# Amplification and saturation of turbulent magnetic field in collapsing
primordial gas clouds

Sho Higashi Department of Physics, Konan University, Okamoto, Kobe, Japan
Hajime Susa Department of Physics, Konan University, Okamoto, Kobe, Japan
Christoph Federrath Research School of Astronomy and Astrophysics, Australian
National University, Canberra, ACT 2611, Australia  Australian Research
Council Centre of Excellence in All Sky Astrophysics (ASTRO3D), Canberra, ACT
2611, Australia  Gen Chiaki National Astronomical Observatory of Japan,
Mitaka, Tokyo, Japan  National Institute of Technology, Kochi College, Monobe,
Kochi, Japan

###### Abstract

Recent numerical studies suggest that magnetic fields play an important role
in primordial star formation in the early universe. However, the detailed
evolution of the magnetic field in the collapse phase still has uncertainties
because of the complicated physics associated with turbulence in a collapsing
magnetized system. Here, we perform a suite of numerical MHD simulations that
follow the collapse of magnetized, turbulent primordial gas clouds to
investigate the evolution of the magnetic field associated with the
turbulence, assuming a polytropic equation of state with exponent $\gamma_{\rm
eff}$ and with various numerical resolutions. In addition, we generalize the
analytic theory of magnetic field growth/saturation so that it can deal with
various exponents $\gamma_{\rm eff}$ and turbulence energy spectra. We find
that the numerical results are well reproduced by the theory for various
$\gamma_{\rm eff}$ through the collapse phase during the formation of the
first stars. The magnetic field is eventually amplified by a factor of
$10^{12}$ – $10^{15}$ due to kinematic and non-linear turbulent dynamo effects
and reaches 3% – 100% of the equipartition level, depending on $\gamma_{\rm
eff}$. We also find that the transition between the kinematic and non-linear
stages can be analytically estimated. These results indicate that the strong
magnetic field accompanied by supersonic turbulence is a general property and
suggest that it can play a crucial role in the formation of the first stars.

††software: athena++ (Stone et al., 2020; Tomida & Stone, 2023), yt (Turk et
al., 2011)

##  1 Introduction

Theoretical studies in the $\Lambda$CDM paradigm predict that the first
generation of stars, known as population III stars, are formed in host dark
matter minihalos at redshifts $z\gtrsim 10$ (Haiman et al., 1996; Tegmark et
al., 1997; Nishi & Susa, 1999; Fuller & Couchman, 2000; Abel et al., 2002;
Bromm et al., 2002; Yoshida et al., 2003). In such a primordial environment,
the star-forming gas only contains light elements. The cooling rates per
volume due to the chemical species of light elements (e.g., ${\rm H_{2}}$) are
lower than those of heavier elements (e.g., CO). As a result, stars are born
in an environment 2–3 orders of magnitude higher in mass and 2 orders of
magnitude higher in temperature than the present-day star-forming environment,
such as the Milky Way. Therefore, the stars are expected to be typically high
mass $(\sim 10-1000M_{\odot})$ due to the large fragmentation scale and high
mass accretion rate if they form as a single stars in the cloud (e.g., Susa et
al., 1996; Omukai & Nishi, 1998; Yoshida et al., 2008).

On the other hand, a number of recent numerical studies have indicated that
the first stars are likely formed as multiple systems, such as a binary,
triplet, etc. The multiple systems form through the fragmentation of the
accretion disk, which emerges after the formation of the central protostar
(Greif et al., 2011, 2012; Susa, 2013; Susa et al., 2014; Susa, 2019; Hirano &
Bromm, 2017; Inoue & Yoshida, 2020; Sharda et al., 2020, 2021; Chiaki &
Yoshida, 2022; Sugimura et al., 2020, 2023). Studying this fragmentation
behavior is considered to be crucial for determining the initial mass function
(IMF) of the first stars.

Turbulence and magnetic fields are fundamental physical phenomena in the study
of the disk fragmentation. For instance, disk fragmentation is promoted by
turbulence (Clark et al., 2011; Riaz et al., 2018; Wollenberg et al., 2020;
Sugimura et al., 2023) or suppressed by magnetic fields (Machida & Doi, 2013;
Sharda et al., 2020; Sadanari et al., 2021). The efficiency of these effects
depends on the properties of the magnetic field (its strength and structure);
therefore, investigating and clarifying their details are essential for the
study of the first-star formation.

As for the turbulence, previous studies indicated that it is amplified by
gravitational contraction (e.g., Vázquez-Semadeni et al., 1998; Sur et al.,
2010, 2012; Federrath et al., 2011b; Robertson & Goldreich, 2012; Birnboim et
al., 2018; Mandal et al., 2020; Hennebelle, 2021). In addition to these
studies, we have recently confirmed that the turbulence becomes supersonic,
even if the initial velocity fluctuations are relatively small at the onset of
a collapse phase (Higashi et al., 2021). We also found that the amplified
turbulent velocities are saturated with a constant Mach number by a balance
between the energy injection via gravitational contraction and dissipation
(Higashi et al., 2022). Since the subsonic turbulence with a Mach number
$\mathcal{M}\sim 0.5$ is often generated due to the accretion flow of gas
along the filaments into halos at the onset of the cloud collapse (Wise &
Abel, 2007; Greif et al., 2008; Klessen & Hennebelle, 2010), the gas cloud
cores in minihalos may usually have saturated supersonic turbulent velocities
during the collapse phase.

As for the magnetic field, one of the representative generation mechanisms of
a very weak seed magnetic field in the early universe is introduced by
Biermann (1950, so-called Biermann battery mechanism). In the primordial
environment, a non-parallel component of the gradient of electron density and
pressure due to the virial shock at the halo formation or the spatial
variations of ionization in the presence of temperature gradients can generate
the seed field ($\sim 10^{-18}$ G) (e.g., Kulsrud & Anderson, 1992; Lazarian,
1992; Gnedin et al., 2000), or strong radiation flux from the neighboring
star-forming clouds may also generate the seed field (Ando et al., 2010; Doi &
Susa, 2011). Although this magnetic field was considered to be not dynamically
important for the first-star formation, recent studies with an
analytical/numerical approach show that the turbulent motions of gas convert
the kinetic energy into magnetic energy by randomly stretching, twisting, and
folding the magnetic field lines on the viscous scale. This is so-called the
small-scale (turbulent) dynamo effect(e.g., Brandenburg & Subramanian, 2005;
Schekochihin et al., 2007; Federrath, 2016). Due to the short time scales
associated with the small spatial scales, the dynamo operates highly
efficiently and increases the magnetic energy with an exponential growth rate
(e.g., Kazantsev, 1968; Kulsrud & Anderson, 1992; Federrath et al., 2011a;
Schober et al., 2015; Xu & Lazarian, 2016). In collapsing primordial gas
clouds, earlier numerical studies show that the magnetic field is amplified
due to both a small-scale dynamo and due to compression during the collapse
(e.g., Sur et al., 2010, 2012; Federrath et al., 2011b; Turk et al., 2012).
Additionally, recent studies further investigate the detailed properties of
the evolution of the magnetic field in collapsing gas clouds (Xu & Lazarian,
2020; Sharda et al., 2021; Hirano & Machida, 2022; McKee et al., 2020; Stacy
et al., 2022; Sadanari et al., 2023).

In this paper we comprehensively investigate the growth/saturation of magnetic
fields due to turbulence for various effective polytropic exponent
$\gamma_{\rm eff}$ of the barotropic EoS, $P\propto\rho^{\gamma_{\rm eff}}$,
by means of high-resolution numerical simulations and analytic calculations.
We show that the magnetic field always grows to an equipartition level in the
collapsing primordial star-forming cloud.

This paper is structured as follows. We describe the setup of the numerical
simulations in Section 2 and present our results in Section 3. We introduce
earlier analytical/theoretical estimates and generalize them in Section 4. We
compare the analytic calculations to our simulation results in Section 5, and
present the results of the numerical resolution study in Section 6. Section 7
is devoted to discussion. We finally summarize the main points of this work in
Section 8.

##  2 Numerical Method

We perform a suite of high-resolution three-dimensional collapse simulations
by using the adaptive mesh refinement (AMR) magneto-hydrodynamics (MHD) code
athena++ (Stone et al., 2020)111https://github.com/PrincetonUniversity/athena.
In order to solve the MHD equations with self-gravity in Cartesian
coordinates, we choose a 2nd-order accurate Piecewise Linear Method (PLM) for
the spatial reconstruction and a 2nd-order accurate Runge-Kutta method for the
time integration, using an HLLD Riemann solver (Miyoshi & Kusano, 2005). The
basic equations are

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mbox{\boldmath$v$})$ $\displaystyle=$ $\displaystyle 0,$
(1) $\displaystyle\frac{\partial\rho\mbox{\boldmath$v$}}{\partial
t}+\rho(\mbox{\boldmath$v$}\cdot\nabla)\mbox{\boldmath$v$}$ $\displaystyle+$
$\displaystyle\nabla P=$ (2)
$\displaystyle(\mbox{\boldmath$B$}\cdot\nabla)\mbox{\boldmath$B$}-\rho\nabla\phi,$
$\displaystyle\frac{\partial E}{\partial
t}+\nabla\cdot[\left(E+P\right)\mbox{\boldmath$v$}$ $\displaystyle-$
$\displaystyle\mbox{\boldmath$B$}(\mbox{\boldmath$B$}\cdot\mbox{\boldmath$v$})]=$
(3) $\displaystyle-$
$\displaystyle\rho\mbox{\boldmath$v$}\cdot\nabla\phi-\Lambda+\Gamma,$
$\displaystyle\frac{\partial\mbox{\boldmath$B$}}{\partial
t}-\nabla\times(\mbox{\boldmath$v$}\times\mbox{\boldmath$B$})$
$\displaystyle=$ $\displaystyle 0.$ (4)

In these equations, $E,~{}\rho,~{}\mbox{\boldmath$v$},~{}\mbox{\boldmath$B$}$,
and $\phi$ are the total energy density, mass density, velocity, magnetic
field, and gravitational potential, respectively. The total energy density is
given by $E=e+\rho v^{2}/2+B^{2}/2$, where $e$ is the thermal energy density.
The total pressure $P$ is composed of $P=p+B^{2}/2$, where $p$ is the gas
thermal pressure. The gravitational potential $\phi$ is obtained from the
Poisson’s equation, using the full multigrid (FMG) method, implemented by
Tomida & Stone (2023) with the multipole boundary condition,

$\nabla^{2}\phi=4\pi G\rho.$ (5)

To satisfy the divergence-free constraint of the magnetic field
$\left(\nabla\cdot\mbox{\boldmath$B$}=0\right)$, athena++ adopts the
constrained transport (CT) scheme. It enables us to preserve the divergence of
the magnetic field to zero, to the machine precision.

Our simulations in this paper are performed with the assumption of ideal MHD.
This means that we explicitly neglect the non-ideal MHD effects of the
magnetic field, such as Ohmic dissipation and/or ambipolar and Hall diffusion,
because these dissipation effects can be safely neglected, at least at the
core scale during the collapse phase of the primordial gas cloud (Maki & Susa,
2004, 2007; Susa et al., 2015; Nakauchi et al., 2019; Sadanari et al., 2023).

We also assume a barotropic Equation of State (EoS) for simplicity. In this
model, the temperature evolution of the collapsing gas is calculated by using
a polytropic model as in Higashi et al. (2021, 2022). The relation between the
thermal pressure $p$ and the density $\rho$ of the gas is expressed as

$p\propto\rho^{\gamma_{\rm eff}},$ (6)

with the effective polytropic exponent $\gamma_{\rm eff}$. We perform
simulations with 3 difference polytropic indices of

$\displaystyle\gamma_{\rm eff}=\begin{cases}1.2,&\\\
1.1&(\text{primordial}),\\\ 1.0&(\text{isothermal}).\end{cases}$ (7)

Here, the model for $\gamma_{\rm eff}=1.1$ mimics the primordial gas in
minihalos (Omukai & Nishi, 1998). The other two models are used for
comparison.

As an initial condition, we generate a gravitationally unstable, magnetized
Bonner-Ebert sphere with a central density $\rho_{\rm peak,0}=4.65\times
10^{-20}~{}\mathrm{g~{}cm^{-3}}$, and a uniform temperature of
$T_{0}=200~{}\mathrm{K}$. The initial cloud mass is $1700M_{\odot}$. The
computational domain has a size of 5 times the initial cloud radius with
periodic boundary conditions.

By following earlier numerical results (e.g., Federrath et al., 2010;
Federrath, 2013; Federrath et al., 2016), we add an initial velocity field
with a turbulence spectrum typical of supersonic turbulence, $E(k)\propto
k^{-2}$ (Burgers, 1948). Here, $E(k)$ is the 1D turbulent kinetic energy power
spectrum with wavenumber $k$ in Fourier space. This turbulent velocity field
corresponds to a root-mean-square Mach number of 0.5, and is composed of a
mixture of solenoidal (divergence-free) and compressive (solenoidal-free)
modes, with a 2:1 ratio, called the natural mixture (Federrath et al., 2008,
2010).

We also add a uniform magnetic field along the $z$-axis with the initial
strength $\mbox{\boldmath$B$}_{0}=10^{-9}$ G. Note that the strength and
statistical properties of the dynamo-generated magnetic fields are independent
of the seed field structure (Seta & Federrath, 2020). Additionally,
theoretical studies (e.g., McKee et al., 2020) have shown that the very weak
seed field ($\lesssim 10^{-19}$ G) generated by the Biermann battery
(Biermann, 1950) can be amplified quickly to moderate values ($\sim 10^{-8}$
G) due to the turbulent dynamo. Hence our choice of initial magnetic field
strength and structure is reasonable, but not critical to the final magnetic
field strength and structure.

All of our numerical simulations start with a base grid of $512^{3}$ grid
cells. Using the AMR technique, we enforce to resolve the Jeans length by
progressively refining cells as the collapse proceeds. In the highest-
resolution simulations, the Jeans length is resolved by at least 128 grid
cells. We call this number the “Jeans parameter”, $N_{\rm J}$, as in Higashi
et al. (2021, 2022). We note that previous work has shown that this is
sufficient to resolve dynamo amplification during the collapse (Sur et al.,
2010; Federrath et al., 2011b). When the peak densities reach
$10^{-4}~{}\mathrm{g~{}cm^{-3}}$, the minimum spatial resolutions are
$1.91^{-5}$ AU for $\gamma_{\rm eff}=1.0$, $7.65\times 10^{-5}$ AU for
$\gamma_{\rm eff}=1.1$, and $1.22\times 10^{-3}$ AU for $\gamma_{\rm
eff}=1.2$, respectively. We also perform higher/lower-resolution runs, as a
resolution study for the $\gamma_{\rm eff}=1.1$ case, with $N_{\rm J}=$ 32,
64, and 256.

We terminate the simulations when the peak density, $\rho_{\rm peak}$, reaches
$10^{-4}~{}\mathrm{g~{}cm^{-3}}$. We continue only the simulation with
$\gamma_{\rm eff}=1.1$ and $N_{\rm J}=128$ for another 2 orders of magnitude
in higher density, because this is our main simulation model, which allows us
to study the saturation phase of the magnetic field. We use the yt toolkit
(Turk et al., 2011)222https://yt-project.org/ to analyze the simulation data.

##  3 Results

When we analyze our simulation results, we use physical quantities only in a
sphere with a radius of half the Jeans length, $L_{\rm J}$, as Higashi et al.
(2021, 2022) did. Hereafter, we call this sphere “Jeans volume”. Also, all
averaged quantities in this paper are cell volume-weighted averages.

Here, we briefly describe the methodology to evaluate $L_{\rm J}$ in every
data snapshot (see also Higashi et al., 2022).

  1. 1.

Cut out an overdense region, with densities of $\rho\geq\rho_{\rm peak}/16$
centered on the position of the density maximum.

  2. 2.

Calculate the center of mass $\mbox{\boldmath$r$}_{\rm c}$ and ‘tentative’
Jeans length $\left(=\sqrt{\pi c_{\rm s,t}^{2}/G\rho_{\rm mean,t}}\right)$
from the average density $\rho_{\rm mean,t}$, and sound speed $c_{\rm s,t}$ in
the overdense region, and then assume half of this Jeans length to be the
tentative core radius $r_{0}$.

  3. 3.

Recalculate the tentative Jeans length and radius $r_{1}$ as done in step 2,
but now within the sphere of radius $\mbox{\boldmath$r$}_{\rm c}$ centered on
$r_{0}$.

  4. 4.

With the new cloud center and radius, go back to step 3, and repeat it until
$|r_{1}-r_{0}|$ becomes smaller than the smallest cell width.

  5. 5.

Obtain $r_{1}$ as the core radius $r_{1}(=L_{\rm J}/2)$ after its value
converges.

Figure 1: Evolution of the turbulent velocity $v_{\rm turb}$ as a function of
mean density $\rho_{\rm mean}$ in the Jeans volume. The solid curves with
different colors denote the simulation results for different $\gamma_{\rm
eff}$. The dashed lines are the analytic estimates of the saturation
velocities in Higashi et al. (2022, eq. 17).

###  3.1 Overall Evolution of the Fields

Figure 1 shows the evolution of the turbulent velocity as a function of mean
density $\rho_{\rm mean}$ within the Jeans volume in each snapshot. The solid
curves denote the simulation results and the dashed lines indicate our
analytic estimates presented by Higashi et al. (2022, eq. 17).

The turbulent velocity within the Jeans volume $V_{\rm J}$ is defined as
(Higashi et al., 2021, 2022)

$v_{\rm turb}^{2}\equiv\sum_{\leq L_{\rm J}/2}\frac{V_{i}}{V_{\rm
J}}(\mbox{\boldmath$v$}_{i}-\mbox{\boldmath$v$}_{{\rm rad},i})^{2},$ (8)

where $V_{i},~{}\mbox{\boldmath$v$}_{i},~{}\mbox{\boldmath$v$}_{\rm{rad},i}$
are the cell volume, the velocity, and smoothed radial velocity of the $i$th
cell, respectively. $\mbox{\boldmath$v$}_{{\rm rad},i}$ is estimated as
follows:

  1. 1.

Calculate the radial velocity profile within the Jeans volume with $N_{\rm
rad}$ radial bins. To remove the fluctuations of the radial velocity around
the scale of the cell size, $N_{\rm rad}$ should be smaller than the Jeans
parameter so that we set $N_{\rm rad}=16$.

  2. 2.

Interpolate the profile at the position of each cell to estimate the smoothed
radial velocity.

As we found in our previous (purely hydrodynamical) simulations (Higashi et
al., 2021, 2022), the turbulent velocities increase and saturate at the level
of a few times the sound speed, also in the present MHD simulations. Although
the difference in growth rates among the three models is not clear because the
dynamic range of the rapid increase regime is narrow, the growth rates are
almost consistent with our analytic estimate in Higashi et al. (2021, eq. 14).
The reason for this is that the analytically estimated growth rates depend on
$\gamma_{\rm eff}$, but their values are very similar. In addition, the
saturation levels of the turbulent velocities are also consistent with the
analytic estimate in Higashi et al. (2022, eq. 17).

This means the effects of the back-reaction of the magnetic field onto the
flow, the energy transportation between the magnetic field and the turbulent
motion, have little impact on the growth and saturation of the turbulent
velocity.

Figure 2: Upper panel: evolution of the RMS magnetic field strength $B_{\rm
RMS}$ as a function of the mean density $\rho_{\rm mean}$ in the Jeans volume.
Lower panel: magnetic field strength divided by $\rho_{\rm mean}^{2/3}$ to
compensate the effect of the field growth due to flux freezing and contraction
during the collapse. The solid curves with different colors denote the
simulation results for different $\gamma_{\rm eff}$. The dashed lines show the
approximated growth rate for each model.

Figure 2 shows the evolution of the root-mean-square (RMS) magnetic field
strength $B_{\rm RMS}$ within the Jeans volume (upper panel) and its division
by $\rho_{\rm mean}^{2/3}$ to compensate the effect of the amplification due
to flux-freezing (bottom panel). We see that the initially very weak
($=10^{-9}$ G) magnetic field increases during the collapse and is eventually
amplified by a factor of $10^{12}-10^{15}$. The field strength at the end of
the collapse, especially for $\gamma_{\rm eff}=1.1$ (red curve) is consistent
with the results of a previous study (Sadanari et al., 2023).

In the lower panel, focusing on the primordial case, the growth rate is
$B_{\rm RMS}\propto\rho^{0.9}$. This is almost identical to the simulation
result of cosmological simulations by Turk et al. 2012
($B\propto\rho^{0.89}$).

On the other hand, comparing the 3 models, while the growth rates of the
magnetic fields in the very early phase $(\rho_{\rm mean}\lesssim
10^{-17}~{}\mathrm{g~{}cm^{-3}})$ are almost identical, they are different in
the late stage, depending on $\gamma_{\rm eff}$. We will discuss the
differences in the growth rates in this stage (kinematic stage) in detail, in
later Sections 4 and 5.

###  3.2 Turbulence/Magnetic Field Spectra

To take a closer look at the detailed evolution of the velocity and magnetic
fields, we compute both fields in Fourier space (hereafter $k$-space) using a
fast Fourier transform code. To make the analysis easier, we remap the
simulation results onto a uniform grid with the same cell size as the highest
refinement level. We here note that, unlike our earlier work Higashi et al.
(2021), we cut out a cube with $2L_{\rm J}$ on each side to capture the
emergence of the Kazantsev spectrum ($\propto k^{3/2}$, Kazantsev, 1968;
Kulsrud & Anderson, 1992) at low-$k$ end. This procedure is similar to the
Fourier analysis method applied in Federrath et al. (2011b).

We calculate the specific turbulent kinetic energy $\left(\varepsilon_{\rm
turb}\right)$ and the specific magnetic energy $\left(\varepsilon_{\rm
mag}\right)$ spectra as

$\displaystyle\varepsilon_{\rm turb}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\left|\hat{\mbox{\boldmath$v$}}_{\rm
turb}(\mbox{\boldmath$k$})\right|^{2}4\pi k^{2}dk,$
$\displaystyle\varepsilon_{\rm mag}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\frac{\left|\hat{\mbox{\boldmath$B$}}(\mbox{\boldmath$k$})\right|^{2}}{4\pi\rho_{\rm
mean}}4\pi k^{2}dk,$ (9)

where $\mbox{\boldmath$\hat{v}_{\rm turb}(\mbox{\boldmath$k$}$})$ and
$\mbox{\boldmath$\hat{B}(\mbox{\boldmath$k$}$})$ are the turbulent velocity
and the magnetic field in $k$-space, and $k$ is the wave vector.

We plot the obtained spectra for the $\gamma_{\rm eff}=1.1$ run in Figure 3.
The upper panel depicts the specific turbulent kinetic energy spectra, while
the lower panel shows the specific magnetic energy spectra. Each color denotes
the snapshot corresponding to $\rho_{\rm mean}$. We are interested in the
scale around the Jeans length, so we normalize the wavenumber by $k_{\rm J}$,
where $k_{\rm J}$ is the wavenumber at the Jeans length.

In the upper panel, we see that the specific turbulent kinetic energy spectra
roughly keep the initial spectra $\propto k^{-2}$ (black dashed line) in the
early phase $(\rho_{\rm mean}\lesssim 10^{-17}~{}\mathrm{g~{}cm^{-3}})$, as
seen in Higashi et al. (2021). In the lower panel, the magnetic spectra have a
mild peak at a smaller k (larger spatial scale) as collapse proceeds. This is
good evidence that a small-scale dynamo is in action in the collapsing gas
cloud. In the kinematic regime, corresponding to the green and blue curves
$(10^{-17}<\rho_{\rm mean}<10^{-8}~{}\mathrm{g~{}cm^{-3}})$, there is some
evidence of a Kazantsev spectrum ($\propto k^{3/2}$; black dot-dashed line)
emerging in the low-$k$ limit $(\log[k/k_{\rm J}]<0.3)$. At the same scale
range, the turbulent kinetic energy is consistent with the Kolmogorov
($\propto k^{-5/3}$; black dotted line) spectrum. These results are in
agreement with Brandenburg & Subramanian (2005) and Federrath et al. (2011b).
The simulation results also suggest that the spectrum is likely flatter than
the Kazantsev one at the low-$k$ end, after saturation (purple curve) in this
collapse simulation, as also shown in the uniform-grid simulation of Seta &
Federrath (2020). While one can suggest that the growth still continues at a
lower $k$ than $k_{\rm J}$ after the saturation occurs at $k_{\rm J}$, we do
not fully understand the reason for the flattening of the spectrum after the
saturation. We plan to investigate this further in a future study.

Figure 3: Evolution of the specific turbulent kinetic energy spectra (upper
panel) and the specific magnetic energy spectra (lower panel), as a function
of wavenumber $k$, scaled by the Jeans wavenubmer $k_{\rm J}$. The different
colors denote different $\rho_{\rm mean}$ to indicate the time evolution
during the collapse. The black dashed, dotted, and dash-dotted lines denote
$\propto k^{-5/3}$, $\propto k^{-2}$, $\propto k^{3/2}$, respectively.

###  3.3 Saturation Level of the Magnetic Field

When the magnetic field strength amplified by the turbulence reaches a
saturation level, the Lorentz force has become strong enough to counter the
stretching, twisting, and folding actions induced by the turbulence during the
exponential growth phase of the turbulent dynamo (Federrath, 2016). The
saturation level can maximally become as high as the equi-partition between
the specific magnetic energy and the specific turbulent kinetic energy.
However, the actual saturation level depends on the properties of the
turbulence, in particular the sonic Mach number and the driving mode of the
turbulence (Federrath et al., 2011a; Achikanath Chirakkara et al., 2021). To
quantify the saturation level of the magnetic fields, we introduce the ratio
$B_{\rm RMS}^{2}/B_{\rm eq}^{2}$ (see Figure 4). $B_{\rm RMS}$ denotes the
root mean square of the field strength in the Jeans volume. $B_{\rm eq}$ is
the equi-partition level of the magnetic field defined as

$B_{\rm eq}\equiv(4\pi\rho_{\rm mean})^{1/2}v_{\rm turb}.$ (10)

The solid curves are drawn by utilizing $v_{\rm turb}$ from the simulation
results (solid curves in Figure 1), while the dashed curves are estimated by
using the analytic saturated turbulent velocities (Higashi et al., 2022)
(dashed lines in Figure 1). As shown in Figure 1, $v_{\rm turb}$ and $v_{\rm
sat}$ are almost identical in the late phase, so that the two ratio are also
almost identical. The dashed curves are smoother than the solid ones by
construction, which enables us to observe the clear saturation behavior.

For comparison, we plot the shaded areas, which denote the saturation levels
obtained from turbulent-in-a-box simulations with an isothermal EoS for
driving Mach numbers $\mathcal{M}=1.0$ (cyan) and 2.0 (magenta) in Federrath
et al. (2011a). The range of the shaded area is given by the difference in
saturation level between turbulence driving modes, that is, compressive or
solenoidal (Federrath et al., 2010); the former gives an estimate for the
lower limit of the expected saturation level, while the latter provides an
estimate of the upper limit. The mixture of the two extreme driving modes
corresponds to an intermediate value.

In Figure 4, we see that the saturation levels depend on $\gamma_{\rm eff}$.
For $\gamma_{\rm eff}=1.1$ and $1.2$, it can be clearly seen that the specific
energy ratio becomes constant at the end of both simulations, whereas it is
unclear for $\gamma_{\rm eff}=1.0$. While a turnover to a likely saturation is
underway in this model, the saturation process has not been completed. These
results indicate that the lower the $\gamma_{\rm eff}$, the higher the
saturation level. As shown in the upper panel of Figure 2, all of the magnetic
field strengths are almost identical, whereas the turbulent velocities depend
on $\gamma_{\rm eff}$. This indicates that the differences in the saturation
levels originate from the differences in the turbulent velocity among the
three runs. In summary, we find that the magnetic energies reach some
fractions of the equi-partition level of about 3% – 100% depending on
$\gamma_{\rm eff}$, until the protostellar density $\left(\simeq
10^{-5}~{}\mathrm{g~{}cm^{-3}}\right)$ is reached. This shows that the
magnetic field strength can amplify to a significant level in the collapsing
phase, which will likely affect the dynamics of the core, and ultimately also
the parental cloud.

Figure 4: Evolution of the ratio between the specific magnetic energy and the
specific turbulent kinetic energy. The solid curves are drawn by utilizing
$v_{\rm turb}$ from simulation results shown in Figure 1. The dashed curves
are estimated using $v_{\rm sat}$ from Higashi et al. (2022). The shaded areas
denote the saturation levels obtained from the results of turbulence-in-a-box
simulations with an isothermal EoS for driving Mach numbers $\mathcal{M}=1.0$
(cyan) and $2.0$ (magenta) in Federrath et al. (2011a).

##  4 Analytic Estimate of the Growth rate

In this section, we introduce analytic estimates for the growth rate of the
magnetic field through gravitational contraction, based on the arguments of Xu
& Lazarian (2020), McKee et al. (2020), and Stacy et al. (2022).

###  4.1 Kinematic Dynamo Stage

In the stage when the magnetic energy is much smaller than the kinetic energy,
the magnetic field lines are stretched, twisted, and folded by turbulence on
the viscous scale $\ell_{\nu}$ (i.e., smallest eddy scale), where the
subscript $\nu$ labels the quantities at the viscous scale. This stage is
called the kinematic dynamo stage, and then the magnetic field is expected to
be amplified with the exponential growth rate, $\Gamma$, as

$B^{2}\propto\exp\left(\int\Gamma dt\right).$ (11)

We assume that the turbulent velocity on the viscous scale $\ell_{\nu}$ is
given by the extrapolation from a certain scale $L$ with a characteristic
scaling law, that is,

$v_{\nu}=v_{\rm L}\left(\frac{\ell_{\nu}}{L}\right)^{\vartheta},$ (12)

where $L$ is the characteristic scale of the turbulence, which is typically of
the order of the Jeans length $L_{\rm J}$ in collapsing gas clouds, and
$\vartheta$ is in the range from 1/3 (Kolmogorov turbulence, Kolmogorov 1941)
to 1/2 (compressible turbulence, Burgers 1948); see also the discussion and
simulations in Federrath et al. (2021). The growth rate of the magnetic
fields, $\Gamma_{\nu}$, at the viscous scale (Schober et al., 2012; Bovino et
al., 2013) is given as

$\Gamma_{\nu}=\frac{v_{\nu}}{\ell_{\nu}}=\left(\frac{v_{\rm turb}}{L_{\rm
J}}\right)\mathrm{Re}^{(1-\vartheta)/(1+\vartheta)},$ (13)

where $\mathrm{Re}$ is the Reynolds number (Haugen et al., 2004; Federrath et
al., 2011b) defined as

$\mathrm{Re}=L_{\rm J}v_{\rm turb}/\nu\simeq N_{\rm J}^{4/3}.$ (14)

Additionally, in collapsing primordial gas clouds, a high ionization degree
enables the condition in which the gas and magnetic flux are fully coupled at
the Jeans scale (Maki & Susa, 2004, 2007; Susa et al., 2015; Nakauchi et al.,
2019). For isotropic, uniform collapse with this condition, i.e., the magnetic
field is randomly oriented in the core, the magnetic field will increase as
$B\propto\rho^{2/3}$ even without the dynamo growth, i.e., purely due to flux
conservation during collapse. Thus, the sum of the growth effects due to flux-
freezing and kinematic dynamo gives

$B_{\rm
kin}^{2}=B_{0}^{2}\xi^{4/3}\exp\left(C_{\Gamma}\int\Gamma_{\nu}dt\right),$
(15)

where $B_{0}$ denotes the initial field strength, and $\xi$ is the ratio
between the density $\rho$ and the initial density $\rho_{0}$. In the above
equation, $C_{\Gamma}$ is a constant, introduced by Stacy et al. (2022),
possibly determined by the properties of turbulence, such as the driving mode
mixture and the magnetic Plandtl number (Federrath et al., 2011a, 2014). In
Stacy et al. (2022), they adopt $C_{\Gamma}=3/4$, whereas we assume
$C_{\Gamma}=0.075$. The reason of this choice is discussed in Section 5. From
Equation (13), assuming the outer scale of turbulence as $L_{\rm J}$, we have

$\int\Gamma_{\nu}dt=\left(\frac{3}{32}\right)^{1/2}\mathrm{Re}^{(1-\vartheta)/(1+\vartheta)}\langle\mathcal{M}_{\rm
turb}\rangle\tau,$ (16)

where $\langle\mathcal{M}_{\rm turb}\rangle$ is the time-averaged turbulent
Mach number in the region where $\rho_{\rm mean}\leq
10^{-12}~{}\mathrm{g~{}cm^{-3}}$. In the above equation, we use the relation

$L_{\rm J}=\left(\frac{\pi c_{\rm
s}^{2}}{G\rho}\right)^{1/2}=\left(\frac{32}{3}\right)^{1/2}c_{\rm s}t_{\rm
ff}$ (17)

We also define

$\tau\equiv\int(1/t_{\rm ff})dt,$ (18)

as in Sur et al. (2010) and Federrath et al. (2011b).

The magnetic field keeps growing exponentially due to the kinematic dynamo
until the magnetic energy becomes comparable to the turbulent kinetic energy
at the viscous scale.

###  4.2 Non-linear Dynamo Stage

When the magnetic energy becomes comparable to the turbulent kinetic energy at
the viscous scale, the suppression of the turbulent dynamo action kicks in.
After the suppression of the growth at the smaller scale, the larger scale
eddies take over the dynamo growth, until the two energies are comparable at
that scale. It is theoretically expected that this cycle continues until the
equipartition is achieved up to the driving scale of the
turbulence.333However, the equipartition level at the driving scale is usually
lower than unity, depending on the magneto-hydrodynamic properties. For
example, Federrath et al. (2011a, 2014) indicated that the equipartition level
depends on the driving mode, the Mach number, and the Plandtl number of the
turbulence.

This stage is called the non-linear dynamo stage. In this stage, the treatment
of the kinematic stage, where only the microscopic diffusion exists and the
magnetic flux is frozen-in on larger scales than the dissipation scale of the
magnetic energy, is no longer valid. The developed MHD turbulence operates the
fast turbulent reconnection and magnetic flux diffusion becomes a non-
negligible effect (see Lazarian et al., 2015, for a review in detail). As a
result, the flux-freezing condition is violated, and the growth rate becomes
smaller (Xu & Lazarian, 2016, 2020).

Xu & Lazarian (2020) derived the equation for a non-linear dynamo in a time-
dependent medium with a self-consistent treatment, and Stacy et al. (2022)
extended their treatment to describe the evolution of the dynamo up to the
late stage of the collapse phase. However, both studies assumed an isothermal
EoS, even though the gas is not necessarily isothermal in the primordial
environment. For example, in Stacy et al. (2022), the thermal evolution for a
number density $\gtrsim 10^{10}{\rm cm^{-3}}$ is almost isothermal (see fig.
11 in their paper), but it is not in Yoshida et al. (2008), at the same
density range. In fact, the detailed thermal evolution shown in different
works in the literature slightly differs from each other, depending on the
implemented chemistry network and radiation transfer prescription. Thus, we
further generalize their arguments by considering the dependence on
$\gamma_{\rm eff}$.

In the non-linear stage, the dynamo growth can occur at $k<k_{\rm p}$, where
$k_{\rm p}$ denotes the wavenumber at the kinetic-magnetic energy
equipartition is achieved. We obtain the time-dependent specific magnetic
energy $\varepsilon_{\rm B}$ by the integral form of the Kazantsev spectrum
(Kazantsev, 1968; Kulsrud & Anderson, 1992) as,

$\varepsilon_{\rm B}=\varepsilon_{\rm B,ref}\left(\frac{k_{\rm p}}{k_{\rm
ref}}\right)^{5/2}\exp\left(\frac{3}{4}\int_{t_{1}}^{t}\Gamma dt\right),$ (19)

where the quantities with the subscript ‘ref’ denote the reference values in
the contracting system, for example, the time-varying core radius of the
collapsing cloud. The quantities with the subscript ‘1’ denote the values at
the onset of the non-linear stage. The flux-freezing condition gives the
growth of the reference specific magnetic energy as $\varepsilon_{\rm
B,ref}=\varepsilon_{\rm B1}\left(\xi/\xi_{1}\right)^{1/3}$. The reference
wavenumber $k_{\rm ref}$ scales inversely with the length as the contraction
proceeds. Considering that the reference scale is the Jeans length, it is
given as $k_{\rm ref}\propto L_{\rm J}^{-1}\propto\xi^{1-\gamma_{\rm eff}/2}$.
So, we obtain

$k_{\rm ref}=\left(\xi/\xi_{1}\right)^{1-\gamma_{\rm eff}/2}k_{\rm ref,1}.$
(20)

On the other hand, using the scaling law in Equation (12) for the turbulent
velocity, we obtain $v_{\rm p}=v_{\rm turb}\left(k_{\rm ref}/k_{\rm
p}\right)^{\vartheta}$. The magnetic energy at $k_{\rm p}$ is given by the
turbulent kinetic energy at the same scale because the energy equipartition is
established for $k>k_{\rm p}$:

$\varepsilon_{\rm B}=\frac{1}{2}v_{\rm p}^{2}=\frac{1}{2}v_{\rm
turb}^{2}\left(\frac{k_{\rm ref}}{k_{\rm p}}\right)^{2\vartheta}.$ (21)

The eddy turnover rate at $k_{\rm p}$ is

$\Gamma=\Gamma_{\rm p}=v_{\rm p}k_{\rm p}.$ (22)

We note that we consider only $\vartheta=1/3$ in the non-linear stage, unlike
the kinematic stage. The reason of this is discussed later (Section 5).

By combining the above two equations and assuming the turbulent velocity is
saturated (i.e., $v_{\rm turb}=v_{\rm sat}$) at the onset of the non-linear
stage, we obtain

$\Gamma\varepsilon_{\rm B}=\frac{1}{2}v_{\rm sat}^{3}k_{\rm ref}$ (23)

Considering the density dependence of $k_{\rm ref}$ and $v_{\rm sat}\propto
c_{\rm s}$, the above equation can be rewritten as

$\Gamma\varepsilon_{\rm
B}=\frac{1}{2}\varepsilon_{1}\left(\frac{\xi}{\xi_{1}}\right)^{\gamma_{\rm
eff}-1/2},$ (24)

where $\varepsilon_{1}=v_{\rm sat,1}^{3}k_{\rm ref,1}$.

Applying $d\ln/dt$ to both sides of Equations (19) and (21), vanishing $k_{\rm
p}$ by combining them, and using Equation (24), we find

$\frac{d\varepsilon_{\rm B}}{dt}=\varepsilon_{\rm
B}a_{\gamma}\frac{d}{dt}\left[\ln\left(\frac{\xi}{\xi_{1}}\right)\right]+\chi\varepsilon_{1}\left(\frac{\xi}{\xi_{1}}\right)^{\beta_{\gamma}},$
(25)

where

$\displaystyle a_{\gamma}$ $\displaystyle=$
$\displaystyle\frac{15}{19}\gamma_{\rm eff}-\frac{41}{57},$
$\displaystyle\beta_{\gamma}$ $\displaystyle=$ $\displaystyle\gamma_{\rm
eff}-1/2,$ $\displaystyle\chi$ $\displaystyle=$ $\displaystyle\frac{3}{38}.$
(26)

We finally obtain the following equation by integrating Equation (25) in the
same manner as in Stacy et al. (2022):

$\varepsilon_{\rm
B}=\left(\frac{\xi}{\xi_{1}}\right)^{a_{\gamma}}\varepsilon_{\rm
B,1}+\frac{\chi\varepsilon_{1}\xi^{a_{\gamma}}}{\xi_{1}^{\beta_{\gamma}}}\int_{t_{1}}^{t}\xi(t^{\prime})^{\beta_{\gamma}-a_{\gamma}}dt^{\prime}.$
(27)

The first term in Equation (27) expresses the growth due to cloud collapse
including the effects of turbulent reconnection diffusion, and the second is
for the dynamo growth. Assuming $\gamma_{\rm eff}=1$, Equation (27) recovers
the identical expression in Xu & Lazarian (2020) and Stacy et al. (2022).

In order to estimate the field strength, we use the same approximation as in
Stacy et al. (2022). Using eq. (B17) in McKee et al. (2020), we have

$\displaystyle\int_{t_{1}}^{t}\xi(t^{\prime})^{\beta_{\gamma}-a_{\gamma}}$
$\displaystyle dt^{\prime}$ $\displaystyle\simeq\frac{2\phi_{\rm ff}t_{\rm
ff,0}}{3\pi(\frac{1}{2}-\beta_{\gamma}+a_{\gamma})}$ (28)
$\displaystyle\times$
$\displaystyle\left(\frac{1}{\xi_{1}^{\frac{1}{2}-\beta_{\gamma}+a_{\gamma}}}-\frac{1}{\xi^{\frac{1}{2}-\beta_{\gamma}+a_{\gamma}}}\right),$

where $t_{\rm ff,0}=\xi^{1/2}t_{\rm ff}$ is the free-fall time at the initial
density and $\phi_{\rm ff}$ is the ratio between the times that the gas
collapses to a star $t_{\rm coll}(\rho)$ and the free-fall time $t_{\rm
ff}(\rho)$. This integration is correct if $\beta_{\gamma}-a_{\gamma}<1/2$ is
satisfied, which is true in the present case (c.f., Equations 26).
Substituting the above expression into Equation (27) and arranging the second
term with Equation (17), we obtain

$\displaystyle\varepsilon_{\rm B}$ $\displaystyle=$
$\displaystyle\left(\frac{\xi}{\xi_{1}}\right)^{a_{\gamma}}\varepsilon_{\rm
B,1}+\frac{\chi\phi_{\rm ff}\mathcal{M}_{\rm sat}v_{\rm
sat,1}^{2}}{2\sqrt{6}\pi\left(\frac{1}{2}-\beta_{\gamma}+a_{\gamma}\right)}$
(29) $\displaystyle\times$ $\displaystyle\left(k_{\rm ref,1}L_{\rm
J,1}\right)\left[\left(\frac{\xi}{\xi_{1}}\right)^{a_{\gamma}}-\left(\frac{\xi}{\xi_{1}}\right)^{\beta_{\gamma}-1/2}\right],$

where $\mathcal{M}_{\rm sat}(\equiv v_{\rm sat,1}/c_{\rm s,1})$ is the
saturation Mach number in Higashi et al. (2022, eq. 18). In the above
equation, all density-dependent quantities are converted to the form of
$Y(\xi)=Y_{1}(\xi/\xi_{1})^{Z}$. We also arrange the equation as

$\varepsilon_{\rm
B}=\left(\frac{\xi}{\xi_{1}}\right)^{a_{\gamma}}\varepsilon_{\rm
B,1}\mathcal{A}_{\gamma,1}^{2},$ (30)

where

$\displaystyle\mathcal{A}_{\gamma,1}^{2}$ $\displaystyle=$ $\displaystyle
1+\frac{\chi\phi_{\rm ff}\mathcal{M}_{\rm
sat}}{\sqrt{6}\pi\left(\frac{1}{2}-\beta_{\gamma}+a_{\gamma}\right)}\left(\frac{\frac{1}{2}v_{\rm
sat,1}^{2}}{\varepsilon_{\rm B,1}}\right)$ (31) $\displaystyle\times$
$\displaystyle\left[1-\left(\frac{\xi}{\xi_{1}}\right)^{\beta_{\gamma}-1/2-a_{\gamma}}\right].$

Here, we utilize the relation $k_{\rm ref,1}L_{\rm J,1}=1$, obtained by
Equation (20). Unless the specific magnetic energy ($\varepsilon_{\rm B,1}$)
is much smaller than the specific turbulent kinetic energy ($v_{{\rm
sat},1}^{2}/2$), the contribution due to dynamo amplification is not
significant (i.e., $\mathcal{A}_{\gamma,1}\sim 1$).

Rewriting Equation (29), we have the analytic estimate of the magnetic field
strength in the non-linear stage,

$B_{\rm
nl}=B_{1}\mathcal{A}_{\gamma,1}\left(\frac{\xi}{\xi_{1}}\right)^{(1+a_{\gamma})/2}.$
(32)

Finally, we need to assess $\xi_{1}$ to distinguish the two stages. We can
complete this task with the argument of the energy ratio. With the definition
for the beginning of the non-linear stage, the specific energies
$\varepsilon_{\nu}$ of the magnetic field and the turbulence are identical at
the viscous scale. Assuming the Kazantsev spectrum ($\propto k^{3/2}$) for the
specific magnetic energy and the Kolmogorov spectrum ($\propto k^{-5/3}$) for
the specific turbulent kinetic energy, we can define the energy ratio at the
onset of the non-linear stage $R_{\rm nl,1}$ as

$R_{\rm nl,1}=\frac{\int_{k_{\rm
J}}^{k_{\nu}}\varepsilon_{\nu}(k/k_{\nu})^{3/2}dk}{\int_{k_{\rm
J}}^{k_{\nu}}\varepsilon_{\nu}(k/k_{\nu})^{-5/3}dk}.$ (33)

Normalizing the wavenumbers by $k_{\rm J}$, the above equation can be written
as

$\displaystyle R_{\rm nl,1}$ $\displaystyle=$ $\displaystyle N_{\rm
J}^{-19/6}\frac{\int_{1}^{N_{\rm J}}N^{3/2}dN}{\int_{1}^{N_{\rm
J}}N^{-5/3}dN}$ (34) $\displaystyle=$ $\displaystyle\frac{4}{15}N_{\rm
J}^{-19/6}\frac{N_{\rm J}^{5/2}-1}{1-N_{\rm J}^{-2/3}},$

where $k/k_{\rm J}=N$ and $k_{\nu}/k_{\rm J}=N_{\rm J}$. For $N_{\rm J}=128$,
Equation (34) gives us $R_{\rm nl,1}\simeq 0.011$. Hence we can evaluate
$\xi_{1}$ as $\xi$ such that $\varepsilon_{\rm mag}/\varepsilon_{\rm
turb}=R_{\rm nl,1}$ is satisfied, i.e., $\xi_{1}=\xi\left(\varepsilon_{\rm
mag}/\varepsilon_{\rm turb}=R_{\rm nl,1}\right)$.

###  4.3 Equipartition

When the specific magnetic energy reaches a certain fraction of the specific
turbulent kinetic energy at the driving scale of the turbulence, the field
growth due to the non-linear dynamo ends. Afterward, the magnetic field
strength can be estimated solely from the specific turbulent kinetic energy.
In addition, the turbulent velocity at the late phase of collapse is generally
saturated and given by a few times the sound speed (Higashi et al., 2022).
Thus, the dependence of the equipartition field strength is

$B_{\rm eq}\propto\left(4\pi\rho\right)^{1/2}v_{\rm sat}\propto\xi^{1/2}c_{\rm
s}\propto\xi^{\gamma_{\rm eff}/2}.$ (35)

For comparison, we summarize the growth rates of the magnetic field in the
form of $B\propto\xi^{X}$ in the non-linear stage without dynamo growth and
the equipartition stage in Table 1.

Table 1: Growth rates of the magnetic field $\left(B\propto\xi^{X}\right)$ in the non-linear stage without dynamo growth and the equipartition stage. $\gamma_{\rm eff}$ | Non-linear | Equipartition
---|---|---
1.0 | $\xi^{0.535}$ | $\xi^{0.500}$
1.1 | $\xi^{0.575}$ | $\xi^{0.550}$
1.2 | $\xi^{0.614}$ | $\xi^{0.600}$

##  5 Comparison between Simulation Results and Analytic Estimates

First, we compare our simulation results with the analytic values obtained by
Equation (15) and Equation (32) in Figure 5. From top to bottom, the panels
depict the comparison for $\gamma_{\rm eff}=1.2$, 1.1, and 1.0, respectively.
The threshold $\xi_{1}$ shown at the upper left in each panel is calculated
using Equation (34).

Figure 5: Comparison between the simulation results and the analytic
estimates. From top to bottom, the panels depict the comparison for
$\gamma_{\rm eff}=1.2$, 1.1, and 1.0, respectively. The $\xi_{1}$ values
(separating the kinematic from the non-linear stage) are shown in the upper
left in each panel. The solid curves are the simulation results. The dashed
curves are obtained by Equation (15). The dotted curves are obtained by
Equation (32).  Figure 6: Scatter plot of the calculated turbulent velocity
spectral indices for $\gamma_{\rm eff}=$1.2 (top), 1.1 (middle), and 1.0
(bottom). The solid curves are a quadratic fit in the range $\xi<\xi_{1}$ by
using the least-square method. The dashed lines denote $\vartheta=1/3$.

In order to measure the spectral indices $\vartheta$ in Equation (16), we fit
the simulation results. In the kinematic stage, the dynamo growth occurs at
$k_{\nu}$, while it occurs at $k_{\rm p}$ in the non-linear stage. In order to
assess them from simulation results, we adopt the relation
$k_{\nu}=0.025k_{\rm turb}\mathrm{Re}^{3/4}$ from Kriel et al. (2022),
assuming $k_{\rm turb}/k_{\rm J}=1$. $k_{\rm p}$ is the wavenumber where
$\varepsilon_{\rm mag}$ is the maximum. As a result, the estimated spectral
indices $\vartheta$ are shown in Figure 6. We also fit these points in the
range $\xi<\xi_{1}$ by a quadratic function using the least-square method and
plot them as solid curves. 444Note that the fitting by a quadratic function
can be valid only in the kinematic stage for $\xi_{1}\ll\infty$. For
$\xi_{1}\rightarrow\infty$, we might need another fitting function. Overall,
the spectral indices (initially $\vartheta=1/2$) gradually decrease as the
collapse proceeds and roughly decline to $\vartheta\sim 1/3$ (dotted lines) in
all models. This is the reason why we assume $\vartheta=1/3$ in Section 4.2.

Going back to Figure 5 and focusing on the dashed curves, these curves are
obtained by using the fitting function in Figure 6. We can see that the
magnetic fields in the kinematic stage for all models are well reproduced by
Equation (15). This result also means that the growth of the magnetic fields
sensitively depends on the spectral index of the turbulent kinetic energy (see
also Schober et al. (2012)).

As shown in Section 4.1, we set $C_{\Gamma}=0.075$. This is the best-fit
parameter chosen for reproducing the simulation results, but possibly
determined by the property of turbulence. To clarify the physical dependence
of $C_{\Gamma}$ is interesting, but is beyond the scope of this work. We will
address this issue in a future study.

In the non-linear stage, the remaining unknown is the value of $\phi_{\rm ff}$
in Equation (29). In Stacy et al. (2022), they assumed $\phi_{\rm ff}$ as a
constant and then calculated its value by integrating the infall velocity of
the gas. In this work, for simplicity, we assume the cloud collapses in a
self-similar fashion, which can be a good approximation of the ideal
collapsing cloud as confirmed by Higashi et al. (2022). While they studied the
non-MHD case, the initial magnetic field strength is very weak in the present
work, and hence the same approximation holds here, i.e., the collapse dynamics
are not much different from the non-magnetized collapse. Since the gas core
density is given by the collapse time $t_{\rm ss}$ as $\rho_{\rm
mean}=\alpha(0)/4\pi Gt_{\rm ss}^{2}$, where $\alpha(0)$ is a constant
depending on $\gamma_{\rm eff}$, $\phi_{\rm ff}$ can be analytically estimated
as

$\phi_{\rm ff}=t_{\rm ss}/t_{\rm
ff}=\left(\frac{32\alpha(0)}{12\pi^{2}}\right)^{1/2}.$ (36)

From Higashi et al. (2022), $\alpha(0)=\left[1.67,3.46,7.19\right]$ for
$\gamma_{\rm eff}=\left[1.0,1.1,1.2\right]$, so we have $\phi_{\rm
ff}=\left[0.67,0.94,1.39\right]$.

Focusing on the dotted curves in Figure 5, we see that the magnetic fields in
the non-linear stage for all models can be reproduced to within a factor of
2.05 by our extended estimate based on Equation (32).

To summarize, the overall magnetic field evolution can be reproduced by our
extended analytic estimates (Equations (15) and (32)) for various different
values of $\gamma_{\rm eff}$. The transition between the kinematic stage and
the non-linear stage can also be estimated analytically from Equation (34)
using the numerical resolution.

##  6 Resolution Study

In this section, we present the results of the resolution study for
$\gamma_{\rm eff}=1.1$ to reinforce the arguments in the earlier sections.
First of all, we compare the turbulent velocities for different numerical
resolutions (Jeans resolution $N_{\rm J}$) in Figure 7. As shown in Section 2,
$N_{\rm J}$ is the number of grid cells per Jeans length , i.e., this number
doubles when the minimum cell width is halved.

Figure 7: The same as Figure 1, but for different numerical Jeans resolutions.
The different colors denote different numbers of grid cells to resolve the
Jeans length.

We find that the turbulent velocity converges to the saturation velocity
(magenta solid line) in all the resolution cases. As discussed in Higashi et
al. (2021), turbulence with a lower initial Mach number requires a higher
Jeans resolution to capture the growth accurately; for example, $N_{\rm J}\geq
128$ is required for $\mathcal{M}_{0}=0.05$ in our previous results. In this
work, we find that $N_{\rm J}=32$ is enough resolution to accurately capture
the growth of the turbulent velocity for $\mathcal{M}_{0}=0.5$.

Figure 8: Upper panel: evolution of specific energy ratios between the
magnetic energy and turbulent kinetic energy for different Jeans resolutions,
as indicated in the legend. The dot-dashed lines show the estimated energy
ratio at the onset of the non-linear stage. Lower panel: comparison between
the simulation results and the analytic estimates in the kinematic (dashed)
and non-linear (dotted) stages. We multiply factors to separate the different
Jeans resolution cases, to facilitate the comparison between the simulations
and the models. For $N_{\rm J}$ = 256, 128, 64, and 32, the factors are +2,
+1, -1, and -2 in the logarithmic scale, respectively.

Figure 8 shows the evolution of the specific energy ratio (upper panel) and
the comparison between the analytic estimates and our simulation results
(lower panel) for various Jeans resolutions $N_{\rm J}$.

In the upper panel, the solid curves denote our simulation results, and the
dash-dotted lines are the corresponding transition ratios obtained by Equation
(34), which are distributed between 0.0068 and 0.029. We can see that the
models for $N_{\rm J}\geq 64$ reach the non-linear stage, but not for $N_{\rm
J}=32$.

In light of the above results, we present the comparison between our
simulations and the estimated field strength obtained by Equations (15) and
(32) in the lower panel. The transition points in each color between the
dashed and dotted curves correspond to the intersections between the solid
curves and the dash-dotted lines in the upper panel. We add factors to the
results for ease of viewing. We can see that our simulation results can be
well reproduced by our analytic estimates for various resolutions, both in the
kinematic and non-linear stages. The transition between the two stages is also
consistent with the analytic estimate for all resolutions. This fact indicates
the robustness of our models.

##  7 Discussion

The magnetic field plays important roles also in the mass accretion phase
after the collapse phase. In a general scenario, a collapsing primordial cloud
becomes adiabatic when the density reaches $\gtrsim 10^{-5}$ g cm-3, which is
determined by the end of H2 dissociation. This threshold density is almost
independent of the metallicity (Omukai et al., 2005). As a result, a
hydrostatic core forms, known as a protostar, accompanied by an accretion disk
(Larson, 1969; Penston, 1969; Greif et al., 2012). This disk can fragment into
pieces when it becomes gravitationally unstable due to the mass accretion from
the surroundings, resulting in a multiple primordial stellar system (e.g.,
Greif et al., 2012; Stacy et al., 2016; Susa, 2019), similar to present-day
star formation and fragmentation in disks (Kuruwita & Federrath, 2019).

Most of the recent numerical studies (Sharda et al., 2020, 2021; Stacy et al.,
2022) have shown that strong magnetic fields amplified by the dynamo effect
can moderately suppress the fragmentation of the disk, resulting in a
relatively top-heavy mass function555Note that Hirano & Machida (2022) claimed
the magnetic field amplification around the formed protostars may be so
significant that fragmentation may be totally suppressed..

Our results indicate that strong magnetic fields always exist in PopIII
prestellar cores irrespective of $\gamma_{\rm eff}$, i.e., different
cooling/heating prescriptions or finite metallicity cases. Thus, the
theoretical argument that the expected IMF of the first stars would be more
top-heavy than those in non-magnetized calculations becomes more robust than
what is considered so far.

This fact might give some suggestions for the observational traces of the
first stars. If the dynamo-generated magnetic field indeed can suppress the
disk fragmentation, the formed first stars would become more massive than
expected by earlier numerical studies without magnetic effects. This means
that the formation of the low-mass first stars, which can survive until the
present-day universe ($\lesssim 0.8M_{\odot}$), tends not to be favored.

On the other hand, the existence of the strong magnetic field at the onset of
the accretion phase gives advantages to the formation of a massive tight
population III star binary, which is considered to be a progenitor of the
massive black hole binary that emits gravitational waves (GWs). In order to
coalesce to emit GWs within a Hubble time, the separation of the black hole
binary has to be $\lesssim 100R_{\odot}$, which is as large as the expanded
size of the progenitor protostars. In the absence of strong magnetic fields,
transportation of angular momentum from the spin of protostars to the orbital
one results in a widening of the binary system(Kirihara et al., 2023). Thus,
it is hard to form such a close binary. The magnetic breaking in the presence
of a strong magnetic field could potentially counteract this widening
mechanism by the spin-orbit angular momentum transfer. So far, no studies have
examined in detail the conditions for first stellar binary mergers
incorporating the effects of magnetic braking. More in-depth research on this
issue is needed.

##  8 Summary

We have performed high-resolution numerical simulations to comprehensively
study the magnetic field evolution associated with turbulence during the
collapse of prestellar cores in the primordial star formation environment. We
have also generalized analytic estimates for the growth rate of the magnetic
field during the gravitational contraction and compared them with our
numerical simulation results. Finally, we performed additional collapse
simulations for the resolution study to reinforce our arguments.

Our main findings are summarized below.

  1. 1.

Amplification and saturation of the turbulent velocity through the
gravitational contraction are not significantly affected by the back-reaction
from the strong magnetic field (Figures 1 and 7)

  2. 2.

The Kazantsev spectrum ($\propto k^{3/2}$) emerges in the collapsing clouds
until the saturation of the magnetic field (Figure 3).

  3. 3.

The magnetic/kinetic energy ratio at the saturation becomes a constant even if
the turbulent velocity keeps increasing during the collapse, because the
growth rate of the magnetic energy equals that of the turbulent kinetic energy
at late times, and its value depends on $\gamma_{\rm eff}$. The smaller the
$\gamma_{\rm eff}$, the larger the energy ratio, because the saturation
turbulent velocity is smaller for lower $\gamma_{\rm eff}$ (Figure 4).

  4. 4.

The growth rates of the magnetic field due to the turbulent dynamo in both the
kinematic stage and the non-linear stage depend on $\gamma_{\rm eff}$
(Equations 15, 25, and Figure 2).

  5. 5.

The transition specific energy ratio between the magnetic energy and the
turbulent kinetic energy can be obtained solely from the numerical Jeans
resolution (Equation 34), and it does not depend on $\gamma_{\rm eff}$.

  6. 6.

The magnetic field evolution in the collapsing gas clouds is well described by
Equations (15) and (32) (c.f., Figures 5 and 8), for various $\gamma_{\rm
eff}$ and numerical Jeans resolutions.

Combined with the series of our previous studies (Higashi et al., 2021, 2022),
we confirm the growth of tangled magnetic fields as well as the supersonic
turbulence in collapsing primordial clouds. This result is supported not only
by high-resolution numerical simulations, but also by detailed analytic
calculations. We find that the simulation and analytical results agree well
with each other. Moreover, we extended the theory to general $\gamma_{\rm
eff}$ cases, for which the agreement still holds at a very good level,
providing support for the validity of the simulation and analytical results.

We are gratful to the anonymous referee for the careful reading of the
manuscript and the constructive comments. We thank Kengo Tomida, Kazuyuki
Sugimura, Tsuyoshi Inoue, Kazutaka Kimura, Kohei Inayoshi, Sadanari Kenji
Eric, Shingo Hirano, and Raiga Kashiwagi for fruitful discussions and useful
comments. This work was supported by JST SPRING, Grant Number JPMJSP2117. We
are thankful for the support by Ministry of Education, Science, Sports and
Culture, Grants-in-Aid for Scientific Research No.22K03689. C. F. acknowledges
funding provided by the Australian Research Council (Future Fellowship
FT180100495 and Discovery Project DP230102280), and the Australia-Germany
Joint Research Cooperation Scheme (UA-DAAD). C. F. further acknowledges high-
performance computing resources provided by the Leibniz Rechenzentrum and the
Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale
project 10391), the Australian National Computational Infrastructure (grant
ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the
framework of the National Computational Merit Allocation Scheme and the ANU
Merit Allocation Scheme. Numerical calculations in this work were carried out
on Yukawa-21 at the Yukawa Institute Computer Facility and Cray XC50 at Center
for Computational Astrophysics, National Astronomical Observatory of Japan.
Computations and analysis described in this work were performed using the
publicly-available Athena++ and yt codes, which are the products of a
collaborative effort of many independent scientists from numerous institutions
around the world.

## References

  * Abel et al. (2002) Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93, doi: 10.1126/science.295.5552.93
  * Achikanath Chirakkara et al. (2021) Achikanath Chirakkara, R., Federrath, C., Trivedi, P., & Banerjee, R. 2021, Phys. Rev. Lett., 126, 091103, doi: 10.1103/PhysRevLett.126.091103
  * Ando et al. (2010) Ando, M., Doi, K., & Susa, H. 2010, ApJ, 716, 1566, doi: 10.1088/0004-637X/716/2/1566
  * Biermann (1950) Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Birnboim et al. (2018) Birnboim, Y., Federrath, C., & Krumholz, M. 2018, MNRAS, 473, 2144, doi: 10.1093/mnras/stx2426
  * Bovino et al. (2013) Bovino, S., Schleicher, D. R. G., & Schober, J. 2013, New Journal of Physics, 15, 013055, doi: 10.1088/1367-2630/15/1/013055
  * Brandenburg & Subramanian (2005) Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1, doi: 10.1016/j.physrep.2005.06.005
  * Bromm et al. (2002) Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23, doi: 10.1086/323947
  * Burgers (1948) Burgers, J. 1948, in Advances in Applied Mechanics, Vol. 1, Advances in Applied Mechanics, ed. R. Von Mises & T. Von Kármán (Elsevier), 171–199, doi: https://doi.org/10.1016/S0065-2156(08)70100-5
  * Chiaki & Yoshida (2022) Chiaki, G., & Yoshida, N. 2022, MNRAS, 510, 5199, doi: 10.1093/mnras/stab2799
  * Clark et al. (2011) Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V. 2011, ApJ, 727, 110, doi: 10.1088/0004-637X/727/2/110
  * Doi & Susa (2011) Doi, K., & Susa, H. 2011, ApJ, 741, 93, doi: 10.1088/0004-637X/741/2/93
  * Federrath (2013) Federrath, C. 2013, MNRAS, 436, 1245, doi: 10.1093/mnras/stt1644
  * Federrath (2016) —. 2016, Journal of Plasma Physics, 82, 535820601, doi: 10.1017/S0022377816001069
  * Federrath et al. (2011a) Federrath, C., Chabrier, G., Schober, J., et al. 2011a, Phys. Rev. Lett., 107, 114504, doi: 10.1103/PhysRevLett.107.114504
  * Federrath et al. (2021) Federrath, C., Klessen, R. S., Iapichino, L., & Beattie, J. R. 2021, Nature Astronomy, 5, 365, doi: 10.1038/s41550-020-01282-z
  * Federrath et al. (2016) Federrath, C., Klessen, R. S., Iapichino, L., & Hammer, N. J. 2016, arXiv e-prints, arXiv:1607.00630, doi: 10.48550/arXiv.1607.00630
  * Federrath et al. (2008) Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79, doi: 10.1086/595280
  * Federrath et al. (2010) Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. M. 2010, A&A, 512, A81, doi: 10.1051/0004-6361/200912437
  * Federrath et al. (2014) Federrath, C., Schober, J., Bovino, S., & Schleicher, D. R. G. 2014, ApJ, 797, L19, doi: 10.1088/2041-8205/797/2/L19
  * Federrath et al. (2011b) Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011b, ApJ, 731, 62, doi: 10.1088/0004-637X/731/1/62
  * Fuller & Couchman (2000) Fuller, T. M., & Couchman, H. M. P. 2000, ApJ, 544, 6, doi: 10.1086/317187
  * Gnedin et al. (2000) Gnedin, N. Y., Ferrara, A., & Zweibel, E. G. 2000, ApJ, 539, 505, doi: 10.1086/309272
  * Greif et al. (2012) Greif, T. H., Bromm, V., Clark, P. C., et al. 2012, MNRAS, 424, 399, doi: 10.1111/j.1365-2966.2012.21212.x
  * Greif et al. (2008) Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2008, MNRAS, 387, 1021, doi: 10.1111/j.1365-2966.2008.13326.x
  * Greif et al. (2011) Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75, doi: 10.1088/0004-637X/737/2/75
  * Haiman et al. (1996) Haiman, Z., Thoul, A. A., & Loeb, A. 1996, ApJ, 464, 523, doi: 10.1086/177343
  * Haugen et al. (2004) Haugen, N. E., Brandenburg, A., & Dobler, W. 2004, Phys. Rev. E, 70, 016308, doi: 10.1103/PhysRevE.70.016308
  * Hennebelle (2021) Hennebelle, P. 2021, A&A, 655, A3, doi: 10.1051/0004-6361/202141650
  * Higashi et al. (2021) Higashi, S., Susa, H., & Chiaki, G. 2021, ApJ, 915, 107, doi: 10.3847/1538-4357/ac01c7
  * Higashi et al. (2022) —. 2022, ApJ, 940, 38, doi: 10.3847/1538-4357/ac9b0c
  * Hirano & Bromm (2017) Hirano, S., & Bromm, V. 2017, MNRAS, 470, 898, doi: 10.1093/mnras/stx1220
  * Hirano & Machida (2022) Hirano, S., & Machida, M. N. 2022, ApJ, 935, L16, doi: 10.3847/2041-8213/ac85e0
  * Inoue & Yoshida (2020) Inoue, S., & Yoshida, N. 2020, MNRAS, 491, L24, doi: 10.1093/mnrasl/slz160
  * Kazantsev (1968) Kazantsev, A. P. 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Kirihara et al. (2023) Kirihara, T., Susa, H., Hosokawa, T., & Kinugawa, T. 2023, ApJ, 950, 188, doi: 10.3847/1538-4357/acd1e0
  * Klessen & Hennebelle (2010) Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17, doi: 10.1051/0004-6361/200913780
  * Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Kriel et al. (2022) Kriel, N., Beattie, J. R., Seta, A., & Federrath, C. 2022, MNRAS, 513, 2457, doi: 10.1093/mnras/stac969
  * Kulsrud & Anderson (1992) Kulsrud, R. M., & Anderson, S. W. 1992, ApJ, 396, 606, doi: 10.1086/171743
  * Kuruwita & Federrath (2019) Kuruwita, R. L., & Federrath, C. 2019, MNRAS, 486, 3647, doi: 10.1093/mnras/stz1053
  * Larson (1969) Larson, R. B. 1969, MNRAS, 145, 271, doi: 10.1093/mnras/145.3.271
  * Lazarian (1992) Lazarian, A. 1992, A&A, 264, 326
  * Lazarian et al. (2015) Lazarian, A., Eyink, G., Vishniac, E., & Kowal, G. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140144, doi: 10.1098/rsta.2014.0144
  * Machida & Doi (2013) Machida, M. N., & Doi, K. 2013, MNRAS, 435, 3283, doi: 10.1093/mnras/stt1524
  * Maki & Susa (2004) Maki, H., & Susa, H. 2004, ApJ, 609, 467, doi: 10.1086/421103
  * Maki & Susa (2007) —. 2007, PASJ, 59, 787, doi: 10.1093/pasj/59.4.787
  * Mandal et al. (2020) Mandal, A., Federrath, C., & Körtgen, B. 2020, MNRAS, 493, 3098, doi: 10.1093/mnras/staa468
  * McKee et al. (2020) McKee, C. F., Stacy, A., & Li, P. S. 2020, MNRAS, 496, 5528, doi: 10.1093/mnras/staa1903
  * Miyoshi & Kusano (2005) Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315, doi: 10.1016/j.jcp.2005.02.017
  * Nakauchi et al. (2019) Nakauchi, D., Omukai, K., & Susa, H. 2019, MNRAS, 488, 1846, doi: 10.1093/mnras/stz1799
  * Nishi & Susa (1999) Nishi, R., & Susa, H. 1999, ApJ, 523, L103, doi: 10.1086/312277
  * Omukai & Nishi (1998) Omukai, K., & Nishi, R. 1998, ApJ, 508, 141, doi: 10.1086/306395
  * Omukai et al. (2005) Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627, doi: 10.1086/429955
  * Penston (1969) Penston, M. V. 1969, MNRAS, 144, 425, doi: 10.1093/mnras/144.4.425
  * Riaz et al. (2018) Riaz, R., Bovino, S., Vanaverbeke, S., & Schleicher, D. R. G. 2018, MNRAS, 479, 667, doi: 10.1093/mnras/sty1635
  * Robertson & Goldreich (2012) Robertson, B., & Goldreich, P. 2012, ApJ, 750, L31, doi: 10.1088/2041-8205/750/2/L31
  * Sadanari et al. (2021) Sadanari, K. E., Omukai, K., Sugimura, K., Matsumoto, T., & Tomida, K. 2021, MNRAS, 505, 4197, doi: 10.1093/mnras/stab1330
  * Sadanari et al. (2023) —. 2023, MNRAS, 519, 3076, doi: 10.1093/mnras/stac3724
  * Schekochihin et al. (2007) Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007, New Journal of Physics, 9, 300, doi: 10.1088/1367-2630/9/8/300
  * Schober et al. (2012) Schober, J., Schleicher, D., Federrath, C., et al. 2012, ApJ, 754, 99, doi: 10.1088/0004-637X/754/2/99
  * Schober et al. (2015) Schober, J., Schleicher, D. R. G., Federrath, C., Bovino, S., & Klessen, R. S. 2015, Phys. Rev. E, 92, 023010, doi: 10.1103/PhysRevE.92.023010
  * Seta & Federrath (2020) Seta, A., & Federrath, C. 2020, MNRAS, 499, 2076, doi: 10.1093/mnras/staa2978
  * Sharda et al. (2020) Sharda, P., Federrath, C., & Krumholz, M. R. 2020, MNRAS, 497, 336, doi: 10.1093/mnras/staa1926
  * Sharda et al. (2021) Sharda, P., Federrath, C., Krumholz, M. R., & Schleicher, D. R. G. 2021, MNRAS, 503, 2014, doi: 10.1093/mnras/stab531
  * Stacy et al. (2016) Stacy, A., Bromm, V., & Lee, A. T. 2016, MNRAS, 462, 1307, doi: 10.1093/mnras/stw1728
  * Stacy et al. (2022) Stacy, A., McKee, C. F., Lee, A. T., Klein, R. I., & Li, P. S. 2022, MNRAS, 511, 5042, doi: 10.1093/mnras/stac372
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Sugimura et al. (2020) Sugimura, K., Matsumoto, T., Hosokawa, T., Hirano, S., & Omukai, K. 2020, ApJ, 892, L14, doi: 10.3847/2041-8213/ab7d37
  * Sugimura et al. (2023) —. 2023, arXiv e-prints, arXiv:2307.15108, doi: 10.48550/arXiv.2307.15108
  * Sur et al. (2012) Sur, S., Federrath, C., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2012, MNRAS, 423, 3148, doi: 10.1111/j.1365-2966.2012.21100.x
  * Sur et al. (2010) Sur, S., Schleicher, D. R. G., Banerjee, R., Federrath, C., & Klessen, R. S. 2010, ApJ, 721, L134, doi: 10.1088/2041-8205/721/2/L134
  * Susa (2013) Susa, H. 2013, ApJ, 773, 185, doi: 10.1088/0004-637X/773/2/185
  * Susa (2019) —. 2019, ApJ, 877, 99, doi: 10.3847/1538-4357/ab1b6f
  * Susa et al. (2015) Susa, H., Doi, K., & Omukai, K. 2015, ApJ, 801, 13, doi: 10.1088/0004-637X/801/1/13
  * Susa et al. (2014) Susa, H., Hasegawa, K., & Tominaga, N. 2014, ApJ, 792, 32, doi: 10.1088/0004-637X/792/1/32
  * Susa et al. (1996) Susa, H., Uehara, H., & Nishi, R. 1996, Progress of Theoretical Physics, 96, 1073, doi: 10.1143/PTP.96.1073
  * Tegmark et al. (1997) Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1, doi: 10.1086/303434
  * Tomida & Stone (2023) Tomida, K., & Stone, J. M. 2023, ApJS, 266, 7, doi: 10.3847/1538-4365/acc2c0
  * Turk et al. (2012) Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ, 745, 154, doi: 10.1088/0004-637X/745/2/154
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * Vázquez-Semadeni et al. (1998) Vázquez-Semadeni, E., Cantó, J., & Lizano, S. 1998, ApJ, 492, 596, doi: 10.1086/305064
  * Wise & Abel (2007) Wise, J. H., & Abel, T. 2007, ApJ, 671, 1559, doi: 10.1086/522876
  * Wollenberg et al. (2020) Wollenberg, K. M. J., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2020, MNRAS, doi: 10.1093/mnras/staa289
  * Xu & Lazarian (2016) Xu, S., & Lazarian, A. 2016, ApJ, 833, 215, doi: 10.3847/1538-4357/833/2/215
  * Xu & Lazarian (2020) —. 2020, ApJ, 899, 115, doi: 10.3847/1538-4357/aba7ba
  * Yoshida et al. (2003) Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645, doi: 10.1086/375810
  * Yoshida et al. (2008) Yoshida, N., Omukai, K., & Hernquist, L. 2008, Science, 321, 669, doi: 10.1126/science.1160259