text:
# Fragmentation of a Molecular Cloud Core versus Fragmentation of the Massive Protoplanetary Disk in the Main Accretion Phase

## 1 Introduction

It is widely accepted that binary and multiple stars form as a result of fragmentation in collapsing molecular cloud cores (e.g., Bodenheimer et al., 2000). In the last two decades, fragmentation of the molecular cloud core has been investigated through numerical simulations by many authors. Criteria for the fragmentation of isothermally collapsing clouds has been investigated by Miyama, Hayashi, & Narita (1984); Boss (1993); Boss & Myhill (1995); Tsuribe & Inutsuka (1999). Their criteria have converged as far as the isothermal phase is concerned; a cloud with $`\alpha 0.20.5`$ collapses into fragments depending little on $`\beta `$, even when the initial density and velocity distributions differ. Here, $`\alpha `$ and $`\beta `$ denote the ratio of thermal energy to gravitational energy and that of rotation energy to gravitational energy, respectively. Tsuribe & Inutsuka (1999) pointed out that the criteria for fragmentation corresponds to the formation of a flat disk with flatness greater than $`4\pi `$. On the other hand, a cloud with $`\alpha 0.20.5`$ collapses self-similarly and shows no sign of fragmentation. Hanawa & Matsumoto (1999) and Matsumoto & Hanawa (1999) investigated deformation of the self-similarly collapsing cloud in search of the possibility that deformation of the central cloud to a bar might trigger fragmentation. The growth of the bar mode is slow compared to the timescale of the collapse, i.e., $`\mathrm{\Delta }\rho _{\mathrm{max}}^{0.177}`$, where $`\mathrm{\Delta }`$ and $`\rho _{\mathrm{max}}`$ denote the amplitude of the bar mode and the maximum density of the cloud, respectively. The bar may indeed fragment, but only at a later stage.

As the collapse proceeds, the cloud core becomes optically thick and the efficiency of radiative cooling decreases. The temperature starts increasing when the central density exceeds the critical density of $`10^{13}\mathrm{g}\mathrm{cm}^3`$. This increase in temperature results in the formation of a quasi-static core, i.e., the first core of Larson (1969). The first core grows by accreting gas, and this accretion phase persists for long enough for the core to fragment. Thus, the dynamics of the cloud changes qualitatively at the critical density. Stability against fragmentation is also likely to change at this critical density. In fact, the first core has been shown by recent simulations to be very unstable taking account of the change in temperature (Burkert, Bate, & Bodenheimer, 1997; Nelson, 1998; Sigalotti, 1998; Klapp & Sigalotti, 1998; Boss et al., 2000). Their simulations, however, assume rather small $`\alpha `$. When $`\alpha `$ is small, the cloud is Jeans unstable and fragments easily (Tsuribe & Inutsuka, 1999). The fragmentation of the first core may be due to the small $`\alpha `$ assumed. It is still unanswered whether the first core fragments when the initial cloud has a moderately large $`\alpha `$. Thus, we investigate the fragmentation of the cloud with focus on the cloud with large $`\alpha `$ of 0.765. The initial cloud is only 1.1 times more massive than the critical Bonnor-Ebert sphere (Ebert, 1955; Bonnor, 1956), which is an equilibrium state of the isothermal cloud.

The Bonnor-Ebert sphere provides a good fit to the density distribution of a dark globule. According to recent near-infrared observations (Alves, Lada, & Lada, 2001), the model with $`\xi =6.9\pm 0.2`$ gives the best fit for B68, where $`\xi `$ denotes the non-dimensional radius. Similarly, those with $`\xi =12.5\pm 2.6`$ and $`7.0\pm 0.3`$ give a good fit for B335 (Harvey et al., 2001) and the Coalsak (Racca, Gómez, & Kenyon, 2002), respectively. When $`\xi >6.45`$, the Bonnor-Ebert sphere is unstable against collapse. Thus, our initial model can be applied to these globules.

The model can also be applied to cloud cores embedded with molecular clouds. The masses of such cores evaluated from C18O luminosity are similar to the virial masses (e.g., Onishi et al., 1996). This implies that the cores are gravitationally bound and nearly in equilibrium, and accordingly, the parameter $`\alpha `$ is only slightly less than unity.

Recent numerical studies have lacked any survey of model parameters. In this study, 225 models with different rotation speed, rotation law, and amplitude of bar mode perturbation are considered. The simulations show many types of fragmentation, some of which are new. The main features of each type and the territory of each type in the parameter space are discussed.

In this paper, the collapse and fragmentation of molecular cloud cores is investigated using a nested grid. The nested grid has high spatial resolution near the center of the computation domain and allows fragmentation to be followed without violating the Jeans condition (Truelove et al., 1997). In §2, the models of cloud cores are introduced. In §3, the methods of numerical simulations are presented. In §4, the results are shown and the fragmentation is classified. In §5, the origins of different types of fragmentation are discussed, the simulations are compared with earlier numerical works, the implications on binary formation are related. The paper is concluded in §6.

## 2 Models

As a model for molecular cloud cores, we consider Bonnor-Ebert spheres, which belong to a sequence of equilibrium-state isothermal spherical clouds confined by external pressure (Ebert, 1955; Bonnor, 1956). Given the external pressure ($`P_{\mathrm{ex}}`$) and the sound speed ($`c_s`$), the Bonnor-Ebert sphere is stable against collapse only when the central density is lower than the critical value, $`14.0P_{\mathrm{ex}}/c_s^2`$. The critical Bonnor-Ebert sphere is used as a template for the model clouds examined in this study. In the models, the initial density distribution is given by

$$\rho (r)=\rho _c\varrho _{\mathrm{BE}}(r/a)$$
(1)
and

$$a=c_s\left(\frac{f}{4\pi G\rho _c}\right)^{1/2},$$
(2)
where $`r`$, $`f`$, $`\rho _c`$ and $`G`$ denote the radius, density enhancement factor, the initial central density, and the gravitational constant, respectively. The function $`\varrho _{\mathrm{BE}}`$ denotes the density distribution of the critical Bonnor-Ebert sphere, and can be approximated as

$$\varrho _{\mathrm{BE}}(\xi )=1\frac{\xi ^2}{6}+\frac{\xi ^4}{45}+𝒪(\xi ^6).$$
(3)
The critical Bonnor-Ebert sphere has radius of $`\xi =6.45`$. A density enhancement factor of $`f=1.1`$ is assumed in typical models because observed molecular clouds are nearly in the virial equilibrium. This slight density enhancement collapses rotating clouds when the initial cloud rotates slowly.

The initial central density is set at $`\rho _c=1\times 10^{19}\mathrm{g}\mathrm{cm}^3`$, which corresponds to a number density of $`n_c=2.6\times 10^4\mathrm{cm}^3`$ for the assumed mean molecular weight of 2.3. An initial temperature of $`T=10\mathrm{K}`$ is assumed, and hence $`c_s=0.19\mathrm{km}s^1`$. The radius and mass of the cloud are thus $`R_c=0.144\mathrm{pc}`$ and $`M=3.24M_{}`$ for $`f=1.1`$.

The initial velocity includes only the $`\phi `$-component, and the angular velocity depends on $`R`$ and $`\phi `$ in cylindrical coordinates $`(R,\phi ,z)`$. The angular velocity $`\mathrm{\Omega }`$ is expressed as

$$\mathrm{\Omega }(R,\phi )=\left[\mathrm{\Omega }_0+\mathrm{\Omega }_2\mathrm{cos}\left(2\phi \right)+\mathrm{\Omega }_3\left(\frac{R}{a}\right)\mathrm{cos}\left(3\phi \right)\right]\left[1+2C\left(\frac{R}{a}\right)^2\right]^{1/2},$$
(4)
where $`\mathrm{\Omega }_0`$ denotes the amplitude of the global rotation, and $`\mathrm{\Omega }_2`$ and $`\mathrm{\Omega }_3`$ denote the amplitudes of the velocity perturbation of $`m=2`$ and 3. The parameter $`C`$ specifies the dependence of $`\mathrm{\Omega }`$ on $`R`$. When $`C=0`$, the angular velocity is independent of $`R`$ and rotation is “rigid”. When $`C`$ is larger, the angular velocity decreases more rapidly with increasing in $`R`$. As shown later, fragmentation of the cloud core depends strongly on the parameter $`C`$. A small amplitude for the perturbation of $`m=3`$ is set, such as $`\mathrm{\Omega }_3t_{\mathrm{ff}}=1\times 10^3`$, in all the models, where $`t_{\mathrm{ff}}`$ denotes the initial freefall timescale at the center and is defined as $`(3\pi /32G\rho _c)^{1/2}`$. This $`m`$ = 3 mode breaks the point symmetry with respect to $`R`$ = 0, and accordingly the fragments are slightly asymmetric in this simulation. The model parameters $`\mathrm{\Omega }_0`$, $`\mathrm{\Omega }_2`$, and $`C`$ are varied to investigate the effects of rotation speed, amplitude of the bar mode, and the rotation law on fragmentation of the cloud cores.

To compare our initial models with those of earlier simulations, the ratios $`\alpha =E_{\mathrm{th}}/|E_{\mathrm{grav}}|`$ and $`\beta =E_{\mathrm{rot}}/|E_{\mathrm{grav}}|`$ are computed, where $`E_{\mathrm{th}}`$, $`E_{\mathrm{rot}}`$, and $`E_{\mathrm{grav}}`$ are thermal energy, rotation energy, and gravitational energy, respectively (e.g., Tohline, 1981; Miyama et al., 1984). In our model, the initial cloud has $`\alpha =0.765(1.1/f)`$. The parameter $`\beta `$ is independent of $`f`$, and when $`\mathrm{\Omega }_3=0`$, $`\beta `$ is described by

$$\beta =\beta _Ct_{\mathrm{ff}}^2\left(\mathrm{\Omega }_0^2+\frac{1}{2}\mathrm{\Omega }_2^2\right),$$
(5)
where the coefficient $`\beta _C`$ is a function of $`C`$ as shown in Figure 1. When a cloud rotates rigidly ($`C=0`$), it has $`\beta _C=0.892`$. When $`C1`$, $`\beta _C`$ decreases approximately in proportion to $`C^{1/2}`$. The decrease in $`\beta _C`$ is due to slow rotation in the outer part of the cloud.

When $`\alpha `$ and $`\beta `$ are small, the cloud is unstable against fragmentation. When the cloud has uniform density and rotates rigidly, the criterion for fragmentation is $`\alpha 0.5`$ (e.g., Tsuribe & Inutsuka, 1999). When the cloud is centrally peaked and the axis ratio is 1.5, the criterion is $`\alpha 0.45`$ for low $`\beta `$ (Boss, 1993). When the cloud is more oblate, i.e., the axis ratio is 2.0, the criterion is $`\alpha 0.33`$ for low $`\beta `$. This suggests that a cloud of $`\alpha 0.5`$ is stable against fragmentation. In this study, the possibility of fragmentation of a cloud with $`\alpha >0.7`$ is examined.

The dynamical evolution of a cloud is followed taking account of the self-gravity and gas pressure. The magnetic field is neglected for simplicity. The gas temperature is assumed to be 10 K below the critical density $`\rho _{\mathrm{cr}}=2\times 10^{13}\mathrm{g}\mathrm{cm}^3`$ ($`n_{\mathrm{cr}}=5.24\times 10^{10}\mathrm{cm}^3`$), and to increase adiabatically in proportion to $`\rho ^{2/5}`$ above it. In other words, a barotropic equation of state is assumed, as expressed by

$$P=\{\begin{array}{cc}c_s^2\rho \hfill & \mathrm{for}\rho <\rho _{\mathrm{cr}}\hfill \\ c_s^2\rho _{\mathrm{cr}}(\rho /\rho _{\mathrm{cr}})^{7/5}\hfill & \mathrm{for}\rho \rho _{\mathrm{cr}}\hfill \end{array}.$$
(6)
This change in temperature reproduces the formation of the adiabatic core, which corresponds to the first core of Larson (1969). The value of the critical density $`\rho _{\mathrm{cr}}`$ is taken from the numerical results of Masunaga, Miyama, & Inutsuka (1998), who studied the spherical collapse of molecular cloud cores with radiative hydrodynamics.

## 3 Numerical Methods

In the simulations, the hydrodynamical equation and Poisson equation are solved by a finite difference method with second-order accuracy in time and space. A nested grid is employed to solve the central region with higher spatial resolution. The hydrodynamic code for the nested grid was developed by extending the simulation code of Matsumoto & Hanawa (1999). The nested grid consists of concentric hierarchical rectangular grids (Yorke, Bodenheimer, & Laughlin, 1993), and the cell width of each grid decreases successively by a factor of two. In the following, the coarsest grid is labeled level $`l`$ = 1. The $`l`$-th level grid has $`2^{l1}`$ times higher spatial resolution than the coarsest grid. All the fluxes are conserved at the interface between the coarse and fine grids as in the standard adaptive mesh refinement (AMR; Chiang, van Leer, & Powell, 1992). Thus, the total mass is conserved in our computations. The numerical fluxes are obtained by the method of Roe (1981) with modification to solve the isothermal and polytrope gas. A MUSCL approach and predictor-corrector method are adopted for time integration (e.g., Hirsch, 1990). The Poisson equation is solved by a multigrid iteration on a nested gird (Matsumoto & Hanawa, 2003). This code solves self-gravity consistently over all grids with different levels such that “the gravitational field line” is continuous at interfaces between coarse and fine grids. This consistency ensures that the obtained gravitational potential is accurate at least up to the quadrapole moment of a binary. Thus, the gravitational torque induced by a binary is accurately taken into account in our simulation.

Mirror symmetry with respect to the $`z=0`$ plane is employed to reduce computation cost. A fixed boundary condition is set for the surface of $`r=R_c`$, representing a constant external pressure that confines the cloud during evolution. Gas is considered only in $`rR_c`$ when solving the Poisson equation.

In this paper, each grid has $`256\times 256\times 32`$ cubic cells in high-resolution models, and $`128\times 128\times 16`$ cubic cells in low-resolution models in $`(x,y,z)`$. The model parameters of the high-resolution models are shown in Table 1. The other models shown in this paper are the low-resolution models. The nested grid consists of grids of 5 levels at the initial stage. A new finer grid is introduced to maintain the Jeans condition of $`\lambda _\mathrm{J}/4>h`$ with ample margin (Truelove et al., 1997), where $`\lambda _\mathrm{J}`$ and $`h`$ are the Jeans length and cell width, respectively. Whenever an eighth of the minimum Jeans length ($`\lambda _{\mathrm{J},\mathrm{min}}`$) became smaller than the cell width in the finest grid, a new finer grid was added to the nested grid. This means that a finer grid was added with ample margin of factor 2. Typical models have 14 grid levels at the last stage. The Jeans condition was only violated in these simulations when a high density fragment escaped from the region covered by the finest grid, and the computation was terminated in the stage that this occurred. In the model shown in §4.1.1, evolution was successfully computed up to the stage in which the mass of an adiabatic core (total mass in the region of $`\rho \rho _{\mathrm{cr}}`$) reached 0.07 $`M_{}`$.

## 4 Results

### 4.1 Rigidly Rotating Cloud

In this subsection, a total of 27 models of a rigidly rotating cloud ($`C=0`$) in the region $`0.03\mathrm{\Omega }_0t_{\mathrm{ff}}0.3`$ and $`0\mathrm{\Omega }_2t_{\mathrm{ff}}0.3`$ are presented to study the dependence on $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_2`$.

Figure 2 summarizes the last stages of the 27 models. Each panel denotes the density distribution in the $`z=0`$ plane. The panels are arranged such that $`\mathrm{\Omega }_0`$ increases from left to right and $`\mathrm{\Omega }_2`$ increases from bottom to top. The evolutions of the clouds are classified into five types in the parameter space of the initial rotation ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$) and the initial amplitude of bar mode ($`\mathrm{\Omega }_2t_{\mathrm{ff}}`$). When $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.03`$ (left column), the cloud collapses to form a single disk (disk type collapse). When $`\mathrm{\Omega }_0t_{\mathrm{ff}}=0.3`$ (right column), the cloud never collapses, and instead oscillates. When $`0.05\mathrm{\Omega }_0t_{\mathrm{ff}}0.2`$ (middle three columns), the cloud collapses into several fragments. The last type is further subdivided into three types: disk-bar, ring-bar, and satellite types. In the following, each type is discussed by showing a typical model.

#### 4.1.1 Disk type collapse

In disk type collapse, the cloud collapses almost spherically in the isothermal collapse phase and forms a rotating disk after the central density exceeds the critical density $`\rho _{\mathrm{cr}}`$. The disk grows by accretion and exhibits no sign of fragmentation. The model of ($`\mathrm{\Omega }t_{\mathrm{ff}},\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$) = (0.03, 0.03, 0.0) is shown in Figure 3 as typical model of disk type collapse.

Figure 3a shows the initial stage, showing only the finest three grids ($`3l5`$), i.e., only 1/64 of the full computation volume. The initial cloud has a spherical density distribution, and the cloud collapses almost spherically during the isothermal collapse phase ($`\rho <\rho _{\mathrm{cr}}`$) as a result of the very slow rotation. When $`\rho \rho _{\mathrm{cr}}`$, the central cloud becomes slightly flattened by the rotation and deforms non-axisymmetrically due to perturbation of the bar mode (Figure 3b).

The deformation is evaluated by measuring the moment of inertia,

$$I_{ij}=_{\rho 0.1\rho _{\mathrm{max}}}(r_ir_{g,i})(r_jr_{g,j})\rho (𝒓)𝑑𝒓,$$
(7)
and the total mass,

$$M=_{\rho 0.1\rho _{\mathrm{max}}}\rho (𝒓)𝑑𝒓,$$
(8)
for the gas of $`\rho 0.1\rho _{\mathrm{max}}`$. The subscripts $`i`$ and $`j`$ are coordinate labels, i.e., $`x=r_1`$, $`y=r_2`$, and $`z=r_3`$. The barycenter is defined as

$$r_{g,i}=\frac{1}{M}_{\rho 0.1\rho _{\mathrm{max}}}\rho r_i𝑑𝒓.$$
(9)
The long axis ($`a_l`$), short axis ($`a_s`$), and length along the $`z`$-axis ($`a_z`$) are defined by

$$\left(\begin{array}{c}a_l^2\hfill \\ a_s^2\hfill \\ a_z^2\hfill \end{array}\right)=\frac{1}{2M}\left\{\begin{array}{c}I_{11}+I_{22}+\left[(I_{11}I_{22})^2+4I_{12}^2\right]^{1/2}\hfill \\ I_{11}+I_{22}\left[(I_{11}I_{22})^2+4I_{12}^2\right]^{1/2}\hfill \\ 2I_{33}\hfill \end{array}\right\}.$$
(10)
Figure 4 shows eccentricity ($`a_l/a_s1`$) and flatness \[$`(a_la_s)^{1/2}/a_z1`$\] as functions of the maximum number density $`n_{\mathrm{max}}`$. The eccentricity oscillates with significant amplitude in the range $`10^4\mathrm{cm}^3n_{\mathrm{max}}10^6\mathrm{cm}^3`$, and increases roughly in proportion to $`n_{\mathrm{max}}^{1/6}`$ in the range $`10^6\mathrm{cm}^3n_{\mathrm{max}}10^{10}\mathrm{cm}^3`$ due to bar mode instability (Hanawa & Matsumoto, 1999; Matsumoto & Hanawa, 1999). At the end of the isothermal collapse phase, the long axis is 12% longer than the short axis (eccentricity is 0.122). The flatness, $`a_l/a_z1`$, increases due to spin-up in the central cloud. In the range $`10^8\mathrm{cm}^3n_{\mathrm{max}}10^{11}\mathrm{cm}^3`$, the flatness increases rapidly in proportion to $`n_{\mathrm{max}}^{0.7}`$. At the end of the isothermal collapse phase, the flatness is 0.412 and the long axis is 50% longer than thickness in the $`z`$-direction.

Figure 4 also shows the central angular velocity in unit freefall time $`\mathrm{\Omega }t_{\mathrm{ff}}`$. The angular velocity and the freefall time are measured as $`\mathrm{\Omega }=J_{\mathrm{spin},z}/[M(a_l^2+a_s^2)]`$ and $`t_{\mathrm{ff}}=(3\pi /32G\rho _{\mathrm{max}})^{1/2}`$, where $`J_{\mathrm{spin},z}`$ denotes the $`z`$-component of the total spin angular momentum $`𝑱_{\mathrm{spin}}`$, which is defined as

$$𝑱_{\mathrm{spin}}=_{\rho 0.1\rho _{\mathrm{max}}}(𝒓𝒓_g)\times (𝒗𝒗_g)\rho 𝑑𝒓,$$
(11)
where

$$𝒗_g=\frac{1}{M}_{\rho 0.1\rho _{\mathrm{max}}}𝒗\rho 𝑑𝒓.$$
(12)
This angular velocity represents the average angular velocity in the region of $`\rho 0.1\rho _{\mathrm{max}}`$, and is denoted $`\mathrm{\Omega }_{0.1}`$. Another average angular velocity $`\mathrm{\Omega }_{0.5}`$, defined in the region of $`\rho 0.5\rho _{\mathrm{max}}`$, is also introduced in order to measure the value near the center. At the initial stage, $`\mathrm{\Omega }_{0.1}`$ and $`\mathrm{\Omega }_{0.5}`$ have the same value because the initial cloud rotates rigidly. In the early isothermal collapse phase, the central part rotates faster ($`\mathrm{\Omega }_{0.5}>\mathrm{\Omega }_{0.1}`$). For $`n_{\mathrm{max}}10^8\mathrm{cm}^3`$, the cloud spins up according to $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}n_{\mathrm{max}}^{1/6}`$, in other words, $`\mathrm{\Omega }_{0.5}n_{\mathrm{max}}^{2/3}`$. This spin-up rate coincides with that expected for spherical collapse (Hanawa & Nakayama, 1997). The cloud collapses almost spherically in the isothermal collapse phase, and significant deformation only occurs near the end of the isothermal collapse phase. The maximum value of $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}`$ is 0.0973.

When the central density exceeds $`\rho _{\mathrm{cr}}`$, the infall decelerates near the center and an adiabatic core, the first core, forms. The first core of Larson (1969) is quasi-static and spherical, whereas this adiabatic core is rotating and disk-like. The adiabatic core formation ends the isothermal collapse phase. Infall still continues in the region far from the center, and the adiabatic core accretes gas from the envelope. Thus, the period after adiabatic core formation is called the accretion phase. Figure 3c shows the adiabatic core at the stage of $`tt_{\mathrm{cr}}=216\mathrm{yr}`$, where $`t_{\mathrm{cr}}`$ denotes the time at which $`\rho _{\mathrm{max}}`$ exceeds $`\rho _{\mathrm{cr}}`$. The adiabatic core consists of a flattened central kernel with an envelope of adiabatic gas. The kernel has a radius of 2 AU and a thickness of 1.5 AU. The mass of the kernel is $`M_{13}=7.4\times 10^3M_{}`$, where $`M_N`$ denotes the mass measured for $`n10^N\mathrm{cm}^3`$.

Figure 3d shows the adiabatic core at the stage of $`tt_{\mathrm{cr}}=775\mathrm{yr}`$, after it has begun to accrete gas from the isothermal infalling envelope. The adiabatic core at this stage consists of a flattened central kernel with an extended disk and spiral arms. One of the dense spiral arms evolves into a dense clump, as can be seen to the upper left of the central kernel in Figure 3d. The dense clump falls into the central kernel and merges into it. Figure 3e shows the adiabatic core after the merge of the dense clump ($`tt_{\mathrm{cr}}=838\mathrm{yr}`$). The radius of the adiabatic disk increases to 15 AU by this stage. During the period $`1.0\times 10^3\mathrm{yr}tt_{\mathrm{cr}}1.1\times 10^3\mathrm{yr}`$, this formation and merging of a clump occurs again. During the formation of the clump, the disk shrinks to 10 AU. However after merging, the radius of the adiabatic disk increases to 20 AU. Figure 3f shows the last stage of the simulation ($`tt_{\mathrm{cr}}=1.5\times 10^3\mathrm{yr}`$). The radius of the adiabatic disk has increased to 20 AU, although the radius of the central kernel remains. The masses of the adiabatic core and the kernel are $`6.9\times 10^2M_{}`$ and $`4.2\times 10^2M_{}`$ at this stage.

Figure 5 shows the sizes of the adiabatic core $`(a_l,a_s,a_z)`$ as a function of time. These sizes were evaluated in the same manner as in Figure 4 except that the volume of integration was bounded by $`\rho \rho _{\mathrm{cr}}`$ in equations (7), (8), and (9). Note that the sizes evaluated from the moment of inertia correspond to the scale length of the cloud and are several times smaller than the cloud radius. For example, $`a_l`$, $`a_s`$, and $`a_z`$ are by a factor of $`5^{1/2}=0.447`$ smaller than the radius $`a`$ for a spherical cloud having uniform density. The vertical scale height remains approximately $`a_z=11.5`$ AU, whereas $`a_l`$ and $`a_s`$ increase while oscillating significantly. The average growth rate is roughly $`da_l/dt=da_s/dt3\times 10^3\mathrm{AU}\mathrm{yr}^1`$. These sizes undergo two types of oscillation beyond the steady growth. First, large peaks appear near $`tt_{\mathrm{cr}}=800`$ yr and 1100 yr. This synchronous oscillation corresponds to the formation and merging of the dense clumps in the adiabatic disk. The disk has a small radius during formation, and then expands greatly after merging. Second, $`a_l`$ and $`a_s`$ exhibit an anti-correlation with small amplitude on a short timescale. This anti-correlation is due to intermittent excitation of spiral arms, which transfer angular momentum from the adiabatic disk to the outer infalling envelope. These recurrently excited spiral arms are also seen in Saigo, Hanawa & Matsumoto (2003), in which growth of the first core was investigated.

Figure 5 also shows the mass of the adiabatic core as a function of time. The mass increases monotonically due to accretion. The average accretion rate is approximately $`5\times 10^5M_{}`$ yr-1.

The formation of the dense clump and excitation of the spiral arms are examined here based on the linear stability of Toomre (1964). A thin disk is unstable when the following two conditions are satisfied simultaneously. First, the Toomre $`Q`$-value defined as

$$Q=\frac{c\kappa }{\pi G\mathrm{\Sigma }}$$
(13)
must be smaller than unity in the region of interest, where $`c`$, $`\kappa `$ and $`\mathrm{\Sigma }`$ are the sound speed, the epicyclic frequency, and the surface density, respectively. Second, the unstable region should be larger than the critical Jeans length,

$$\lambda _c=\frac{2c^2}{G\mathrm{\Sigma }}\left[1+\left(1Q^2\right)^{1/2}\right]^1.$$
(14)
These criteria are applied to our simulation by evaluating

$$\mathrm{\Sigma }(x,y)=_{\rho >\rho _{\mathrm{cr}}}\rho (x,y,z)𝑑z,$$
(15)
$$c(x,y)=\frac{1}{\mathrm{\Sigma }(x,y)}_{\rho >\rho _{\mathrm{cr}}}\left[\frac{dP}{d\rho }(x,y,z)\right]^{1/2}\rho (x,y,z)𝑑z,$$
(16)
$$\kappa (x,y)=\left(4\overline{\mathrm{\Omega }}^2+R\frac{d\overline{\mathrm{\Omega }}^2}{dR}\right)^{1/2},$$
(17)
and

$$\overline{\mathrm{\Omega }}(x,y)=\frac{1}{\mathrm{\Sigma }(x,y)}_{\rho >\rho _{\mathrm{cr}}}\frac{v_\phi (x,y,z)}{R}\rho (x,y,z)𝑑z.$$
(18)
Figure 6a shows the distribution of the critical Jeans length $`\lambda _c(x,y)`$ for the stage shown in Figure 3d. The surface density $`\mathrm{\Sigma }`$ takes a positive value within the domain indicated by the closed thick curve, i.e., edge of the adiabatic disk. The Toomre $`Q`$-value is less than unity in the gray regions, and larger than unity in the white regions. It is less than unity only in the central kernel, dense clump and spiral arms. The rest of the disk has $`Q>1`$ and the disk is globally stable against the ring mode. The critical Jeans lengths are less than 1 AU in both the central kernel and the dense clump, which is roughly 1 AU in size. Thus, both are self-gravitationally bounded. Figure 6b shows the same results for the last stage. The $`Q`$-value is less than unity only in the kernel and the spiral arms. The spiral arms have width of 2.5 AU in the central region of $`r5\mathrm{AU}`$, where the critical Jeans length is less than 2.5 AU. In the outer region of $`r10\mathrm{AU}`$, the critical Jeans length is $`57.5`$ AU and the spiral arms have width of 7.5 AU. Thus, the inner and outer spiral arms are self-gravitating and can be supposed to form by gravitational instability.

#### 4.1.2 Disk-bar type fragmentation

In the models of disk-bar type fragmentation, $`(\mathrm{\Omega }_0t_{\mathrm{ff}},\mathrm{\Omega }_2t_{\mathrm{ff}},C)=(0.05,\mathrm{\hspace{0.17em}0.0},\mathrm{\hspace{0.17em}0.0})`$ and $`(0.05,\mathrm{\hspace{0.17em}0.01},\mathrm{\hspace{0.17em}0.0})`$, the cloud collapse is almost axisymmetric in the isothermal collapse phase resulting in the formation of a disk in the beginning of the accretion phase. The disk deforms to a bar-shape in the early accretion phase. Thereafter, the bar fragments into two fragments. This fragmentation is called disk-bar type. The model of ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$) = (0.05, 0.0, 0.0) is shown in Figure 7 as a typical model of disk-bar type fragmentation.

In the initial stage, the cloud undergoes uniform rotation with a small $`m=3`$ perturbation (no $`m=2`$ perturbation). The cloud collapse is almost axisymmetric in the isothermal collapse phase. Figure 7b shows the cloud at the end of the isothermal collapse phase. The central cloud is slightly flattened due to the rotation.

Figure 8 shows the evolution of flatness \[$`(a_la_s)^{1/2}/a_z1`$\], eccentricity ($`a_l/a_s1`$), and angular velocity ($`\mathrm{\Omega }t_{\mathrm{ff}}`$) for the central cloud in the isothermal collapse phase. The flatness increases rapidly in proportion to $`\rho _{\mathrm{max}}^{0.3}`$. At the end of the isothermal collapse phase, a disk forms with flatness of 0.439. The eccentricity remains small because the initial cloud has only a small $`m=3`$ perturbation, with no bar mode perturbation. The angular velocity ($`\mathrm{\Omega }t_{\mathrm{ff}}`$) is relatively large at the beginning, and increases according to $`\mathrm{\Omega }_{0.5}n_{\mathrm{max}}^{1/6}`$ up to $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}=0.452`$ in the isothermal collapse phase. The growth rate of $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}`$ is similar to that shown in Figure 4 for the disk type collapse model.

Figure 7b shows the adiabatic core at the stage of $`tt_{\mathrm{cr}}=301`$ yr. The adiabatic core at this time consists of a flat disk with an envelope. The adiabatic disk deforms to a bar-shape at $`tt_{\mathrm{cr}}=439`$ yr, as shown in Figure 7c, due to self-gravitational instability. The seed of the bar mode is discretization error in this model. Figure 9 shows the distribution of critical Jeans length for the stage shown in Figure 7b. The elliptical disk of 12 AU $`\times `$ 11 AU has a critical Jeans length of $`24`$ AU, except in the region of the central holes, $`r2\mathrm{AU}`$. The disk thus suffers from gravitational instability and deforms to a bar. The bar fragments into two fragments at $`tt_{\mathrm{cr}}500`$ yr, with masses of $`M_{13}=1.0\times 10^2M_{}`$ and $`1.1\times 10^2M_{}`$, and $`M_{14}=4.8\times 10^3M_{}`$ and $`3.0\times 10^3M_{}`$. The fragments rotate around each other and accrete gas from the envelope. Figure 7d shows the fragments at the stage of $`tt_{\mathrm{cr}}=1622`$ yr. Similar to the adiabatic core shown in Figure 3f, each fragment has a central kernel and spiral arms embedded in an extended disk. At this stage, the two fragments have mass of $`M_{13}=3.0\times 10^2M_{}`$ and $`2.4\times 10^2M_{}`$, and $`M_{14}=1.8\times 10^2M_{}`$ and $`1.9\times 10^2M_{}`$.

Figure 10a shows the loci of fragments for the period between the stages of fragmentation and Figure 7d (468 yr $`tt_{\mathrm{cr}}`$ 1622 yr). In this period of $`1.1\times 10^3`$ yr, both fragments rotate approximately three and half times, increasing in separation. The time variation of the separation is shown quantitatively in Figure 11. The separation of the fragments increases from 11.6 AU with significant oscillation due to the eccentricity of the orbits. The fragments attain a maximum separation of 30.6 AU at $`tt_{\mathrm{cr}}=1622`$ yr (at the stage of Figure 7d). As the separation increases, the specific orbital angular momentum also increases by a factor of 4.9.

Figure 10b shows the situation after the maximum separation. After the stage of Figure 7d, the separation begins to decrease, to $`7`$ AU in only $`300`$ yr, as shown in Figure 11. The specific orbital angular momentum also decreases by a factor of 0.28. In this period, the fragments rotate only half a rotation. After the rapid decrease in the separation, the fragments rotate approximately 1.5 revolutions with a nearly constant separation of $`7`$ AU until the last stage. Figures 7e and 7f show the last stage of the simulation at different magnifications. At this stage, the fragments are separated by 6.7 AU, and surrounded by a circumbinary disk with tightly winding spiral arms that transfer the orbital angular momentum of the fragments to the circumbinary disk. The decrease in separation is due to the formation of the circumbinary disk. At the last stage, the two fragments have the same mass of $`M_{14}=2.4\times 10^2M_{}`$.

#### 4.1.3 Ring-bar type fragmentation

In the models showing ring-bar type fragmentation, $`(\mathrm{\Omega }_0t_{\mathrm{ff}},\mathrm{\Omega }_2t_{\mathrm{ff}},C)=(0.1,\mathrm{\hspace{0.17em}0.0},\mathrm{\hspace{0.17em}0.0})`$, $`(0.1,\mathrm{\hspace{0.17em}0.01},\mathrm{\hspace{0.17em}0.0})`$, and $`(0.2,\mathrm{\hspace{0.17em}0.0},\mathrm{\hspace{0.17em}0.0})`$, cloud collapse is almost axisymmetric and a flat disk forms in the isothermal collapse phase. The disk deforms to a ring-shape temporarily, and then to a bar-shape in the early accretion phase. Thereafter, the bar fragments into two or three fragments. This fragmentation is called ring-bar type fragmentation. The model of ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$) = (0.1, 0.0, 0.0) is shown in Figure 12 as a typical example of ring-bar type fragmentation. This type of fragmentation is similar to disk-bar except for the transient formation of a ring and the number of fragments.

The initial stage is the same as for the previous model shown in §4.1.2 except for the initial uniform rotation speed. The cloud collapse is almost axisymmetric in the isothermal collapse phase. Figure 12a shows the cloud at the end of the isothermal collapse phase. The central cloud deforms to a disk-shape due to the rotation.

Figure 13 shows the evolution of flatness \[$`(a_la_s)^{1/2}/a_z1`$\], eccentricity ($`a_l/a_s1`$), and angular velocity ($`\mathrm{\Omega }t_{\mathrm{ff}}`$) for the central cloud in the isothermal collapse phase. These evolutions are similar to those for the disk-bar model. The flatness increases rapidly in proportion to $`\rho _{\mathrm{max}}^{0.4}`$. At the end of the isothermal collapse phase, a thin disk has formed, with flatness of 1.96. The eccentricity remains small because the initial cloud has only a small $`m=3`$ perturbation and no bar mode perturbation. The angular velocity ($`\mathrm{\Omega }t_{\mathrm{ff}}`$) is relatively large from the beginning, and increases according to $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}n_{\mathrm{max}}^{1/6}`$ to a maximum of $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}=0.381`$ at $`n_{\mathrm{max}}=5.46\times 10^9\mathrm{cm}^3`$. When the central cloud is disk-like, $`\mathrm{\Omega }t_{\mathrm{ff}}`$ becomes saturated.

Figure 12b shows the adiabatic core at $`tt_{\mathrm{cr}}=563`$ yr. The central adiabatic core is ring-like and surrounded by a flat isothermal envelope. The ring structure forms due to self-gravitational instability. Figure 14 shows the distribution of critical Jeans length for the stage shown in Figure 12b. The distribution is similar to that for the disk-bar model except that the adiabatic elliptical disk is larger roughly by a factor of two while the critical Jeans length is almost the same. The disk thus suffers from ring instability more strongly than in the disk-bar model.

The ring-shaped adiabatic core deforms to a rotating bar as shown in Figure 12c. The bar is twice as long as that in the disk-bar model, and the bar fragments into three fragments. Figure 12d shows the fragments at $`tt_{\mathrm{cr}}=915`$ yr. The central fragment is the most massive ($`M_{13}=2.0\times 10^2M_{}`$). The other fragments have masses of $`M_{13}=4.8\times 10^3M_{}`$ and $`5.4\times 10^3M_{}`$.

Figure 15 shows the loci of the three fragments. The red and blue fragments rotate at a distance of roughly 15 AU, while the green fragment rotates around the red-blue close binary at distance of roughly 40 AU (see also Figure 12e). The three fragments form a hierarchical triple system.

Figure 12f shows the last stage of the simulation. The three fragments have similar masses of $`M_{13}=1.9\times 10^2M_{}`$ (red), $`M_{13}=1.7\times 10^2M_{}`$ (blue), and $`M_{13}=1.8\times 10^2M_{}`$ (green). The separation between the red and blue fragments is 14 AU, and that between the their barycenters and the green fragment is 41.6 AU at the last stage.

The calculation was terminated at this stage because the green fragment escaped from the region covered by the fine grid ($`l=13`$) to that covered by the coarser grid ($`l=12`$).

#### 4.1.4 Satellite type fragmentation

In the models with $`0.05\mathrm{\Omega }_0t_{\mathrm{ff}}0.2`$ and with large $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, the cloud collapses to form a dense adiabatic core surrounded by an adiabatic disk. The disk suffers from self-gravitational instability and fragments into dense fragments orbiting around the central adiabatic core. The orbiting fragment is called a satellite fragment, and this fragmentation satellite type fragmentation. No appreciable difference between satellite type and disk type is seen up to disk formation. Spiral arms are excited in the adiabatic disk and local dense condensations form. The masses of these condensations exceed the Jeans mass, and evolve into satellite fragments confined by self-gravity. The satellite fragments orbit around the central core and exhibit close encounter, and sometimes merge during the accretion phase. In following, the model of ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$) = (0.1, 0.05, 0.0) is shown as a typical example.

For the model of ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$) = (0.1, 0.05, 0.0), the initial stage is the same as that of the previous model (ring-bar type) except for the amplitude of the bar mode. Figure 16a shows the cloud at the end of the isothermal collapse phase. The dense gas is flattened due to the rotation, and the cross section in the $`y=0`$ plane resembles that of the previous model because of the same initial rotation (see Figure 12a). The central cloud is more elongated than in the previous model due to the large initial amplitude of the bar mode.

Figure 17a shows the distribution of critical Jeans length at the stage of adiabatic bar formation. Only a limited region along the axis is unstable against ring instability in the bar, while the adiabatic core is elongated into a bar.

Figure 16b shows the elongated adiabatic core. The central kernel (dense region in the adiabatic core) rotates differentially and excites spiral arms in the surrounding disk. The elongation of the central kernel precedes excitation of the ring instability as shown by Saigo et al. (2003). In their simulation, a cloud without the bar mode forms a ring-shaped core, while a cloud with the bar mode forms an elongated core with spiral arms. The amplitude of the bar mode separates the satellite type fragmentation from disk-bar and ring-bar types.

The spiral arms sweep up the disk and are wound up. These winding spiral arms then evolve into satellite fragments, which are confined by self-gravity, as shown in Figure 16c. Figure 17b shows the distribution of critical Jeans length at the stage of satellite fragment formation. The adiabatic disk is globally stable against the ring mode, yet is locally unstable in small regions around $`(x,y)(\pm 20\mathrm{AU},7\mathrm{AU})`$. A satellite fragment forms in each of the local unstable regions.

Figure 18a shows the loci of the central fragment and the satellite fragments. The green locus denotes the central fragment, while the blue and red loci indicate the satellite fragments. The blue satellite fragment merges into the central fragment (green) after half a rotation, and the blue and red fragments form a binary system as shown in Figure 16d. The binary excites spiral arms and a new satellite fragment forms 519 yr after the merger. Figure 18b shows the loci after formation of the new satellite fragment. The purple locus denotes the new satellite fragment. The smallest fragment in Figure 16e is the new satellite fragment, which later falls into the green fragment resulting in a binary system again. Figure 16f shows the last stage of the simulation. At this stage, the masses of the fragments are $`M_{13}=2.0\times 10^2M_{}`$ (red) and $`3.0\times 10^2M_{}`$ (green), and the separation between fragments is 26.9 AU.

### 4.2 Dependence on the Rotation Law

In this subsection, the dependence of the collapse and fragmentation on the rotation law specified by the parameter $`C`$ is investigated. The azimuthally averaged angular velocity is independent of $`R`$ (uniform) at $`C=0`$, and decreases with increasing $`R`$ as $`C`$ gets larger.

Figure 19 shows the models of $`\mathrm{\Omega }_0t_{\mathrm{ff}}=0.2`$ for the central 230 AU $`\times `$ 230 AU square at the stage of $`\rho _c\rho _{\mathrm{cr}}`$. Based on the morphology, these models can be classified into three types, i.e., ring, bar, and dumbbell types. Ring-shaped structures are seen in models with large $`C`$, while bar-shaped structures occur in models with large $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$. When both $`C`$ and $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$ are large, the density has two peaks, and forms the dumbbell type. Both the dumbbell and ring forms in the model of $`C0.16`$. Fragmentation, particularly the number of fragments, depends critically on the morphology.

The cloud with $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.03`$ collapses to form an adiabatic disk and exhibits no sign of fragmentation (disk type collapse). On the other hand, a cloud with $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.05`$ fragments by any of the disk-bar, ring-bar, satellite, ring, bar, or dumbbell types as far as it collapses. The parameter $`\mathrm{\Omega }_0t_{\mathrm{ff}}`$ solely specifies whether the cloud fragments, whereas the other parameters specify only the type of fragmentation.

#### 4.2.1 Ring type fragmentation

Ring type fragmentation takes place in models with large $`C`$ and small $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$. A ring forms during the collapse and fragments into more than three fragments.

Figure 20 show the model with $`(\mathrm{\Omega }_0t_{\mathrm{ff}},\mathrm{\Omega }_2t_{\mathrm{ff}},C)=(0.2,\mathrm{\hspace{0.17em}0.0},\mathrm{\hspace{0.17em}1.0})`$ as a typical model of ring type fragmentation. Figure 20a shows the central cloud at the end of the isothermal collapse phase. The evolution is similar for ring type and ring-bar type as far as the isothermal collapse phase is concerned, with a flat disk forming in the cloud center. Figure 20b shows the stage of $`tt_{\mathrm{cr}}=634`$ yr. The central disk suffers from ring instability, and the ring is more prominent than for the ring-bar type. Figure 20c shows the stage of $`tt_{\mathrm{cr}}=879`$ yr, just at the moment of fragmentation. The ring fragments directly into four fragments, whereas the ring deforms into a bar before fragmentation in ring-bar type fragmentation. Figure 20d shows the stage of $`tt_{\mathrm{cr}}=1.26\times 10^3`$ yr (the last stage), in which four fragments can be seen. The loci of the four fragments are shown in Figure 21. The two central fragments (blue and red) exhibit close encounter while the outer fragments (green and purple) rotate with wide orbits. At the last stage, the masses of the fragments are $`M_{13}=2.0\times 10^2M_{}`$ (red), $`9.4\times 10^3M_{}`$ (blue), $`7.0\times 10^3M_{}`$ (green), and $`7.4\times 10^3M_{}`$ (purple). The calculation was terminated here because the Jeans condition was violated after the escape of the outer fragments from the fine grid at $`l=12`$.

#### 4.2.2 Bar type fragmentation

Bar type fragmentation takes place in models with large $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$. The cloud collapses to form a narrow bar. Although the bar fragments into two fragments in many models, the fragments in this model merge to form a central adiabatic core. The merger is due to the small angular momentum of the fragments. After the merger, the adiabatic core excites spiral arms and eventually satellite type fragmentation occurs as shown in §4.1.4. Figure 22 shows the model with $`(\mathrm{\Omega }_0t_{\mathrm{ff}},\mathrm{\Omega }_2t_{\mathrm{ff}},C)=(0.2,\mathrm{\hspace{0.17em}0.2},\mathrm{\hspace{0.17em}0.15})`$ as a typical model of bar type fragmentation followed by satellite type fragmentation.

Figure 22a shows the cloud at the end of the isothermal collapse phase, where it collapses to form a dense bar. Figure 23 shows the evolution of the eccentricity, flatness and rotation of the central part of the cloud in the isothermal collapse phase. The central angular velocity in unit freefall time, $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}`$ and $`\mathrm{\Omega }_{0.1}t_{\mathrm{ff}}`$, increases in proportion to $`n_{\mathrm{max}}^{1/6}`$ and reaches $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}=0.460`$ at $`n_{\mathrm{max}}=9.04\times 10^7\mathrm{cm}^3`$ and $`\mathrm{\Omega }_{0.1}t_{\mathrm{ff}}=0.291`$ at $`n_{\mathrm{max}}=1.14\times 10^8\mathrm{cm}^3`$. Meanwhile, the cloud collapses almost spherically. During $`n_{\mathrm{max}}>10^8\mathrm{cm}^3`$, $`\mathrm{\Omega }t_{\mathrm{ff}}`$ decreases, and the flatness increases in proportion to $`n_{\mathrm{max}}^{1/2}`$ and exceeds unity at $`n_{\mathrm{max}}10^8\mathrm{cm}^3`$. The eccentricity also increases, although with significant oscillation. At the end of the isothermal collapse phase, $`a_l`$, $`a_s`$, and $`a_z`$ are 56.0 AU, 18.4 AU, and 3.84 AU, respectively. The long axis $`a_l`$ is 3.05 times longer than the short axis $`a_s`$.

Figure 22b shows the central cloud at $`tt_{\mathrm{cr}}=210`$ yr. The bar-shaped, adiabatic core is surrounded by the isothermal disk. At this stage, $`a_l`$, $`a_s`$, and $`a_z`$ are 30.4 AU, 2.80 AU, and 2.99 AU, respectively. The long axis $`a_l`$ is 10.9 times longer than the short axis $`a_s`$.

Figure 22c shows the fragmentation of the narrow adiabatic bar at $`tt_{\mathrm{cr}}=510`$ yr. The bar is wound due to the differential rotation, and develops two density peaks. The separation between these peaks is 12.0 AU. Figure 22d shows the cloud at $`tt_{\mathrm{cr}}=1.01\times 10^3`$ yr. These density peaks merge to form a central core surrounded by an adiabatic disk with spiral arms. The spiral arms are the remnants of the wound bar. The disk is supported by centrifugal force and its radius is approximately 20 AU at this stage. The disk radius increases due to the accretion of gas from the infalling envelope, and by $`tt_{\mathrm{cr}}1.5\times 10^3`$ yr, the disk has a radius of $`40`$ AU and deforms into a ring. The ring is connected to the central kernel via the spiral arms, and the three satellite fragments form at the intersections of the ring and the spiral arms, as shown in Figure 22e.

Figure 22f shows the stage of $`tt_{\mathrm{cr}}=1.25\times 10^3`$ yr (the last stage), in which five fragments can be seen. One of the two central fragments is formed by subsequent satellite type fragmentation. These fragments rotate in a close orbit of 17.1 AU. The masses of the fragments of the tight binary at the center are $`M_{13}=1.6\times 10^2M_{}`$ and $`1.7\times 10^2M_{}`$. The masses of the other fragments are $`M_{13}=9.4\times 10^3M_{}`$, $`5.8\times 10^3M_{}`$, and $`1.8\times 10^3M_{}`$ from inner to outer.

#### 4.2.3 Dumbbell type fragmentation

Dumbbell type fragmentation takes place in models with large $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$ and large $`C`$. The cloud collapses to form a dumbbell-shaped dense cloud having two density peaks at the end of the isothermal collapse phase. The dumbbell shape is a hybrid of the ring and bar forms. Although each of the density peaks evolves into a fragment, the fragments often merge as in bar type fragmentation. Satellite fragments form at a later stage whenever dumbbell type fragmentation occurs.

Figure 24 show the model of $`(\mathrm{\Omega }_0t_{\mathrm{ff}},\mathrm{\Omega }_2t_{\mathrm{ff}},C)=(0.2,\mathrm{\hspace{0.17em}0.2},\mathrm{\hspace{0.17em}0.5})`$, as a typical model of dumbbell type fragmentation followed by satellite type fragmentation. The evolution of dumbbell type fragmentation is similar to that of bar type fragmentation. Figure 24a shows the dense dumbbell-shaped cloud at the beginning of the accretion phase ($`tt_{\mathrm{cr}}=151`$ yr). The two density peaks evolve into self-gravitationally bounded fragments, as shown in Figure 24b ($`tt_{\mathrm{cr}}=673`$ yr), surrounded by the isothermal disk. The fragments then merge to form a central kernel. Figure 25 shows the loci of the fragments during the merger.

After the merger, the evolution of dumbbell type fragmentation is very similar to that of bar type. The central kernel is surrounded by a rotation-supported disk, which grows in radius by the accretion of gas from the infalling envelope. Figure 24c shows the disk at $`tt_{\mathrm{cr}}=1.03\times 10^3`$ yr. The disk radius increases to $`40`$ AU at this stage and deforms into a ring. Figure 24d shows the last stage. Satellite fragments form at the intersections of the ring and the spiral arms as in model shown in Figure 22f.

### 4.3 Classification of Fragmentation Processes

The seven types of collapse and fragmentation described above: disk, disk-bar, ring-bar, satellite, ring, bar, and dumbbell types, are summarized as a means of classification. Figure 26 shows the branching of these fragmentation types schematically. In the isothermal collapse phase, cloud collapse is classified into four main types: disk, bar, dumbbell, and ring types. Disk type collapse is subdivided into disk, satellite, bar, and ring in the accretion phase.

For all types, fragments (self-gravitationally confined clumps) form only in the accretion phase. The isothermal collapse phase is so short that the cloud deforms into a disk, bar, dumbbell, or ring-shape, but does not fragment.

Figure 27 summarizes domain of each type of collapse and fragmentation in three-dimensional phase space ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$). The models with $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.03`$ exhibit disk type collapse ($`\times `$) except for one model, while almost all models with $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.05`$ undergo fragmentation. Some exceptional models exhibit oscillation ($``$). The other parameters, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$ and $`C`$, specify the type of fragmentation.

The red symbols denote the models exhibiting satellite type fragmentation. The red cross ($`\times `$) denotes the model exhibiting disk type collapse followed by satellite type fragmentation, as shown in §4.1.1. Similarly the red triangle ($`\mathrm{}`$) denotes the model of bar type fragmentation followed by satellite type fragmentation. Almost all the models proceeds to satellite type fragmentation when the bar mode $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$ of the the initial cloud is significant.

Fragmentation could not be confirmed for the models indicated by filled symbols. Almost all of these models have either a long bar, long dumbbell, or large ring in the beginning of the accretion phase. For the models indicated by filled triangles ($`\mathrm{}`$) and inverted triangles ($`\mathrm{}`$), both the long bar and dumbbell are likely to fragment but could not be confirmed due to violation of the Jeans condition before fragmentation. From comparison with bar and dumbbell type fragmentation as shown in §4.2.2 and §4.2.3, these fragments appear to merge, and satellite type fragmentation should follow the merger. Similarly, for models indicated by filled circles ($``$), the formation of a ring but could be followed but subsequent fragmentation cloud not.

## 5 Discussion

### 5.1 Collapse, Fragmentation, Survival, and Merger

In this subsection, the fate of the collapsing clouds is discussed in detail.

A cloud fragments whenever the cloud collapses and the initial rotation is faster than $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.05`$, independent of the other parameters $`\mathrm{\Omega }_2`$ and $`C`$. The cloud collapses to form a flat disk in the isothermal collapse phase, and fragmentation of the flat disk is occurs by disk-bar, ring-bar, ring, bar, dumbbell, and satellite types. When $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.05`$, a flat disk forms in the accretion phase. The flat disk formed in the isothermal collapse phase fragments, whereas that formed in the accretion phase does not.

The criterion for fragmentation is given by the initial angular velocity. The critical angular velocity, $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.05`$, is evaluated in terms of spin-up of the collapsing cloud. An isothermal cloud in runaway collapse spins up in proportion to $`\mathrm{\Omega }_ct_{\mathrm{ff}}\rho _c^{1/6}`$, where $`\mathrm{\Omega }_c`$ denotes the angular velocity at the center (Hanawa & Nakayama, 1997). After the cloud changes its shape from a sphere to a disk, the angular velocity becomes saturated at $`\mathrm{\Omega }_ct_{\mathrm{ff}}0.5`$ (Matsumoto, Hanawa, & Nakamura, 1997; Matsumoto & Hanawa, 1999). These quantities represent a good index of disk formation.

As shown in §4, the angular velocities $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}`$ in models of disk, disk-bar, ring-bar, and bar types have maximum values of 0.0973, 0.452, 0.381, 0.460, respectively, in the isothermal collapse phase. When $`\mathrm{\Omega }_{0.5}t_{\mathrm{ff}}`$ is close to 0.5 (disk-bar, ring-bar, and bar types), the model forms a disk in the isothermal collapse phase and fragments in the later stages. The small difference between $`\mathrm{\Omega }_c`$ and $`\mathrm{\Omega }_{0.5}`$ is due to differential rotation in the core. The angular velocity is considerably smaller in the rest of the models (disk type).

Applying these quantities, the condition for formation of a disk in the isothermal collapse phase can be evaluated as

$$\mathrm{\Omega }_0t_{\mathrm{ff}}0.5\left(\frac{\rho _c}{\rho _{\mathrm{cr}}}\right)^{1/6}=0.045.$$
(19)
This condition is consistent with the our simulations.

The formation of a flat disk in the isothermal phase depends on the rotation of the central cloud, but not on the rotation law specified by the parameter $`C`$. In the isothermal collapse phase, the central velocity and density become more important to cloud collapse as the cloud shrinks (Matsumoto et al., 1997). During the isothermal (runaway) collapse phase, the mass of the central cloud decreases if defined as the mass contained in the isodensity sphere of $`\rho _{\mathrm{max}}/2`$. The mass of the central cloud is only 0.01 $`M_{}`$ at the end of the isothermal phase. Since fragmentation takes place in the central 0.01 $`M_{}`$, the density and velocity thereof are important. The density and velocity in the envelope have little effect on the initial fragmentation, but is involved in evolution of the fragments in the accretion phase.

In the literature, fragmentation of clouds is typically discussed in terms of the parameters $`\alpha `$ and $`\beta `$ (e.g., Bodenheimer & Burkert, 2001), which are evaluated by volume integration of energy for the entire cloud. As described in equation (5), $`\beta `$ depends not only on the central rotation $`\mathrm{\Omega }_0`$, but also on the rotation in the envelope. As shown in our models, the epoch of flat disk formation depends solely on $`\mathrm{\Omega }_0`$ and is affected little by rotation in the envelope. Thus, the central angular velocity $`\mathrm{\Omega }_0`$, rather than $`\beta `$, describes the formation of the flat disk.

The cloud fragments when it satisfies the criterion of equation (19). The number of fragments and their orbital angular momentum then determines whether the fragments merge or survive. When the cloud fragments into three or more fragments, many of the fragments survive, as shown in ring-bar, satellite and ring types. On the other hand, when the cloud fragments into only two fragments, the fate of the fragments depends on their orbital angular momenta: fragments with high orbital angular momentum survive (e.g., disk-bar type), while fragments with low orbital angular momentum merge (e.g., bar and dumbbell types). In the latter case, satellite fragments form after the merger.

The orbital angular momentum of the fragment depends on the timing of deformation into the bar-shape. When the cloud deforms into a bar in the isothermal collapse phase, the bar does not have sufficient spin angular momentum to be supported by the rotation. In the isothermal collapse phase, the cloud undergoes runaway collapse and is never supported by the rotation (Saigo & Hanawa, 1998). Thus, the orbital motions of the fragments are also never supported by the rotation. Furthermore, the bar loses angular momentum via gravitational torque in the accretion phase. In the model shown in §4.2.3 (dumbbell type), each fragment has a specific orbital angular momentum of $`6.5\times 10^{18}\mathrm{cm}^2\mathrm{s}^1`$ at the stage shown in Figure 24b. This specific orbital angular momentum is only 20% of that required for support by rotation at $`R=20`$ AU. On the other hand, when the cloud deforms into a bar in the accretion phase, the bar is supported by the rotation (e.g., disk-bar type). The adiabatic disk accretes gas with high specific angular momentum and is already supported by rotation prior to its fragmentation.

### 5.2 Comparison with Earlier Numerical Simulations

Bar and dumbbell type fragmentation were also seen in the numerical simulation of Boss et al. (2000). They followed the evolution of clouds having an initial Gaussian density profile by three types of approximations; isothermal equation of state, barotropic equations of state, and the Eddington approximation of radiative transfer. The barotropic model of Boss et al. (2000) exhibits bar type fragmentation. The density profiles shown in their Figs. 5c and 5d are similar to those shown in Figures 22b and 22d here. Using AMR code, they followed the formation of a bar, fragmentation of the bar, and merger of the fragments. This evolution resembles that of the model shown in §4.2.2, although the initial condition is quite different. They terminated the calculation at the stage of adiabatic disk formation after the merger, and the stage of satellite type fragmentation was not shown. Dumbbell type fragmentation is also seen in their Fig. 2a, which is quite similar to Figure 24a here. Dumbbell type fragmentation in their study was computed under the Eddington approximation.

The initial cloud of Boss et al. (2000) is defined by $`\alpha =0.26`$ whereas $`\alpha =0.765`$ in the present study. Their cloud was thus colder, or in other words, more massive. Despite this difference, the fragmentation is very similar. The same mechanism of deformation and fragmentation therefore appears to be valid over a wide range of $`\alpha `$.

Satellite type fragmentation has also been seen in many other simulations (Bonnell, 1994; Burkert et al., 1997; Bate et al., 2002). Bonnell (1994) followed the fragmentation in the second collapse, in which the first core collapses to form the second core. In their simulations, satellite fragments form through interaction of the spiral arms. The satellite fragments form in the same manner both in their and our simulations, even though different situations are considered. Burkert et al. (1997) followed the collapse and fragmentation of molecular cloud cores similar to this paper. They also computed the formation of satellite fragments and followed their orbits using a nested grid. It was not explicitly mentioned whether the nested grid simulation satisfies the Jeans condition. Burkert et al. (1997) and Bate et al. (2002) confirmed the result by independent simulations using an SPH code.

### 5.3 Application to Formation of Binary and Multiple Stars

Some observations have indicated the rotation of molecular cloud cores. Goodman et al. (1993) found that 29 of 43 molecular clouds had significant velocity gradient, corresponding to rigid rotation of $`2\times 10^3<\beta <1.4`$ with typical values of $`\beta 0.02`$. These quantities correspond to $`0.047<\mathrm{\Omega }_0t_{\mathrm{ff}}<1.25`$ with typical values of $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.15`$ for models of $`(\mathrm{\Omega }_2t_{\mathrm{ff}},C)=(0,\mathrm{\hspace{0.17em}0})`$. Unfortunately, the observations were not of sufficient accuracy to specify the rotation law. Our simulations show that collapsing clouds having an initial rotation of $`\mathrm{\Omega }_0t_{\mathrm{ff}}0.05`$ fragment, which is consistent with observed high binary frequency.

Molecular cloud cores have internal motion often interpreted as turbulence. The internal motion should reflect the superposition of various modes of velocity perturbations. When the bar mode of a cloud is a significant, the cloud undergoes satellite type fragmentation as shown in our simulations. Therefore, satellite type fragmentation may be dominant. The satellite fragments merge and scatter while accreting gas. Consequently, the satellites will have various binary separations. This may explain the wide range of separation for young and main sequence binaries (e.g., Mathieu, 1994).

It has recently been suggested that brown dwarfs may be formed by ejection of the seeds of stars from a parent cloud core (Reipurth & Clarke, 2001; Bate et al., 2002). Satellite type fragmentation might be a corresponding case. In many case of satellite type fragmentation, three or more fragments are formed. In these multiple systems, it is possible that a close encounter will eject the fragment from the cloud center. The satellite fragment has speed of $`1\mathrm{km}\mathrm{s}^1`$ at the last stage of the model of ($`\mathrm{\Omega }_0t_{\mathrm{ff}}`$, $`\mathrm{\Omega }_2t_{\mathrm{ff}}`$, $`C`$) = (0.2, 0.2, 0.15) (for a fragment shown in the right side of Figure 22f). The velocity would be reduced substantially before ejection by the gravity of the molecular cloud core. The gravitational potential is evaluated to be $`\psi 2c_\mathrm{s}^2\mathrm{ln}r`$, and the ejection speed would be $`v_{\mathrm{escape}}1\mathrm{km}\mathrm{s}^1`$. The ejected satellite fragment would have a velocity of the order of the escape speed if it exits.

## 6 Summary

The collapse and fragmentation of molecular cloud cores was investigated for the case that the initial cloud is almost in equilibrium, focusing on the effects of rotation speed, rotation law, and bar mode perturbation. The main results are summarized as follows.

A cloud 1.1 times denser than the critical Bonnor-Ebert sphere fragments when rotation of the initial cloud is slowly enough to allow collapse, but still significant, i.e., $`\mathrm{\Omega }_ct_{\mathrm{ff}}0.05`$. The latter condition gives rise to the formation of a flat disk in the isothermal collapse phase. This condition is independent of both the initial amplitude of the bar mode and the initial rotation law.

Six types of fragmentation were identified: disk-bar, ring-bar, satellite, bar, ring, and dumbbell types. The type of fragmentation depends on the initial amplitude of the bar mode and the initial rotation law. The fragments formed via bar or dumbbell types fragmentation merge due to their low angular momenta, and new fragments form via satellite type fragmentation. In other words, a cloud forms satellite fragments whenever the bar mode of the initial cloud is appreciable amplitude. Merger and close encounter of the satellite fragments may result in the wide range of the binary separation.

The authors thank S. Inutsuka, T. Tsuribe, and K. Saigo for valuable discussion. Numerical computations were carried out on VPP5000 at the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan, which is an inter-university research institute of astronomy operated by the Ministry of Education, Culture, Sports, Science and Technology, Japan (MEXT). This study was financially supported in part by Grants-in-Aid for the Encouragement of Young Scientists (12740123, 14740134), for Scientific Research on Priority Areas (A) (13011204) from MEXT, and for Scientific Research (C) (13640237) from the Japan Society of Promotion of Science (JSPS),


format: text

text:
# Scaling Properties of Random Walks on Small-World Networks

## I Introduction

The topological properties of real-world networks have been studied extensively. But an even more intriguing task, and a natural extension of these studies, is to understand how the network structure affects dynamics on the networks.strogatz01  Most people have had the unfortunate experience of catching the flu (an example of disease spreading) or picking up a burning hot plate (an example of neural signal transmission). These phenomena are all examples of dynamics on the special kind of real-world networks that have been found to display “small-world” properties. Of these phenomena, the greatest attention thus far has been given to the study of disease spreading \[see, e. g. albert01 ; kuperman01 ; newman02 ; satorras01 \]. It has even been suggested that the web of sexual contacts has a small-world structure deszo01 ; liljeros01 . Other dynamical models that have been studied on complex networks include the Hodgkin-Huxley model lago00 , Boolean dynamics reka , and the generic synchronization of oscillators pecora . Extensive reviews can be found in Refs. albert01 ; dorogovtsev01 ; strogatz01 .

In this paper, we will present results for a random walk on a small-world network (SWN, defined below). Such random walks may have several applications to real systems. For example, Scala et alscala01  have argued that the conformation space of a lattice polymer has a small-world topology, and hence, that diffusion and random walks on such small-world networks might give insight into relaxation processes such as protein folding.

Much is known about random walks on both regular and random networks (see, e. g., barber70 ; hughes ). In addition, there have been several recent studies of random walks on SWN’s monasson99 ; jespersen00 ; jasch01 ; lahtinen01 ; pandit01 ; scala01 ; octavio ; lahtinen02 . For example, Jasch and Blumenjasch01  and Lahtinen et. al lahtinen01  have studied the average number of distinct sites visited by a walker, and the probability that the walker is at the origin after $`n`$ steps in the limit that the network size $`L\mathrm{}`$ and, hence, the number of shortcuts $`x\mathrm{}`$. In the present work, we present results for the scaling behavior of such quantities as the mean-square displacement, the mean number of distinct sites covered during a random walk, the first-return time abd clarify the procedure for obtaining a scaling collapse for random walks on SWN’s.

## II The random walk

In previous work kulkarni ; almaas02\_1 , we developed a simple probability function approach for the topological properties of complex networks generated according to the small-world modelwatts98 ; newman99 . In this approach, we start from a one-dimensional regular network with periodic boundary conditions and $`L=2N`$ nodes, each node being connected to its $`2k`$ nearest neighbors. Hence, the “degree” of each node is $`2k`$. Next, we add shortcut ends to each node according to a given degree distribution ($`𝒟_q`$), by following the prescription of Refs. molloy ; newman01 . In the present work we use a modified approach almaas02\_1  by using the following degree distribution in one dimension (with $`k=1`$): $`𝒟_q=(1p)𝒫_{(q)}(p)+p𝒫_{(q1)}(p)`$, where $`𝒫_q(\lambda )=\mathrm{exp}(\lambda )\lambda ^q/q!`$ is the Poisson distribution. We then select pairs of shortcut ends at random and connect them to each other, thus creating a shortcut. This network-generating procedure (with the above $`𝒟_q`$) is equivalent to that outlined by Newman and Watts newman99 . The quantity $`kp`$ is the probability that a given site has a shortcut. On average, there will be $`x=kpL`$ shortcuts in the network.

We carry out the random walk on a lattice (Polya walk) as follows: (i) There is only one walker on the small-world network at a time. (ii) The random walker is injected onto a randomly chosen site on the small-world network, a new site for each walker. We will call this site the “origin” of the walk. (iii) At each discrete time step, $`t`$, the walker will jump to a randomly chosen nearest neighbor of its current site $`m`$ with probability $`1/k(m)`$. Here, $`k(m)`$ is the number of nearest neighbors of site $`m`$, i.e. the degree of node $`m`$. (iv) The random walker is allowed to wander the network for a time longer than the “saturation time” for the quantity studied, i.e. the time when that quantity approaches its limiting behavior. (v) We average over different random walkers and realizations of the small-world lattice until the results converge.

## III Scaling Behavior

### III.1 General Form

In earlier work kulkarni ; almaas02\_1 , we demonstrated that the basic probability distribution $`P(m|n;L,p)`$ (the probability that two sites separated by $`n`$ hops before the introduction of shortcuts has a minimal separation of $`m`$ hops when shortcuts are included) scales with $`x`$ in the limit of $`p1`$, for all choices of $`x`$. As a consequence, topological quantities derived from $`P(m|n;L,p)`$ displays a scaling with $`x=pL`$ in the same limit.

We will first state our main result for the scaling of the random walk in general, and then consider specific examples. Let $`O(p,L,t)`$ be some measurable quantity for the random walk on an SWN, which saturates to a finite value $`O_{sat}`$ as $`t\mathrm{}`$. As specific examples, we will discuss the mean-square displacement and the average number of sites covered by the random walker in time $`t`$. We propose that $`O(p,L,t)`$ satisfies the scaling law:

$$O(p,L,t)=O_{sat}(p^2t,pL).$$
(1)
We have numerically confirmed this scaling Ansatz, as described in the rest of this paper. According to Eq. (1), a scaling collapse is observed if one plots the quantity of interest for various choices of the variables $`(p,L)`$, holding $`x=pL=\text{const.}`$, for any choice of $`x`$. This behavior is in contrast to previous statements that a scaling collapse cannot be expected when $`x100`$ lahtinen01 . Note also that the scaling collapse is seen only for fixed values of $`x=pL`$ and that a prerequisite for scaling is that $`p1`$. Previous workers jasch01 ; lahtinen01  have attempted to show scaling collapse for fixed values of $`L`$ (or $`p`$) while varying $`p`$ (or $`L`$), and did not obtain a perfect scaling collapse. By contrast, our present results display a perfect scaling collapse for fixed values of $`x`$. In the following, we present our numerical results for two quantities: the mean number of distinct sites covered, $`N_{\text{cov}}`$, and the mean square displacement, $`r^2`$.

### III.2 Average Number of Distinct Visited Sites

The average number of distinct nodes visited by a random walker, denoted $`N_{\text{cov}}`$, was first studied by Dvoretzky and Erdös dvoretzky  for a random walker on an infinite $`d`$-dimensional regular lattice. They found that $`N_{\text{cov}}\sqrt{t}`$ in 1-$`d`$, $`N_{\text{cov}}t/\mathrm{ln}t`$ in 2-$`d`$ and $`N_{\text{cov}}t`$ for $`d>2`$, in the limit of $`t\mathrm{}`$. How does the small-world network modify this result? Following the arguments of jasch01 , we expect that for short times, $`t\xi ^2`$ (with $`\xi =1/p`$), the random walker will be probing regions of the small-world network that are essentially without shortcuts (linear regions). Hence, the behavior should be similar to the $`p=0`$ case where the walker covers $`N_{\text{cov}}\sqrt{t}`$ sites. For long times, $`tL\xi `$, we expect all the sites in the network to be covered and $`N_{\text{cov}}=L`$. For intermediate times, $`\xi ^2tL\xi `$, the walker spends on the average $`\xi ^2`$ time steps per linear region, before a new region is accessed – the shortcuts act like a kind of branching process for the random walk. As a consequence, in this time regime there should be $`t/\xi ^2`$ segments covered and $`N_{\text{cov}}t`$ jasch01 . Combining these regimes, we obtain the scaling form

$$N_{\text{cov}}=LS(t/\xi ^2;x),$$
(2)
where

$$S(y;x)\{\begin{array}{cc}\hfill \sqrt{y}/x,& y1,\\ \hfill y/x,& 1yx,\\ \hfill 1,& yx,\end{array}$$
(3)
and $`y=t/\xi ^2`$. This is the expected finite-size scaling form for the average number of distinct visited nodes in the small-world model in the limit $`x>1`$.

In Fig. 1, we show a plot of our calculated $`N_{\text{cov}}`$ for $`x=10`$ (circles), $`10^2`$ (squares) and $`10^3`$ (stars). For each value of $`x`$, we have used different values of $`p`$, ranging from $`p=0.0002`$ to $`p=0.01`$, two for each $`x`$-value. The scaling collapse is excellent. These results suggest that the scaling properties of the small-world network determine the scaling of the random walk; hence, to find a scaling collapse, one must keep the average number of shortcuts, $`x`$, constant. This Figure also shows that, as $`x`$ increases, $`N_{\text{cov}}`$ deviates more and more from the $`p=0`$ result.

If, instead, we hold the density of shortcuts constant and vary the system size, we can explore the finite-size effects on the average number of covered sites, $`N_{\text{cov}}`$. In Fig. 2 we plot the succession of curves $`L=10^3`$, $`10^4`$, $`10^5`$ and $`10^6`$, all with $`p=0.01`$. Upon examination, we find that the slope of the resulting scaling curve changes from $`1/2`$ to $`1`$. That is, for small $`t`$ ($`t\xi ^2`$), $`N_{\text{cov}}\sqrt{t}`$, after which, at larger $`t`$, there is a crossover to $`N_{\text{cov}}t`$, and finally, at even larger $`t`$, finite-size effects become apparent. This behavior observed in the simulations agrees very well with the arguments preceding Eq. (3). Note that, when $`\xi L`$, we do not see a crossover to the linear regime in which $`N_{\text{cov}}t`$ before finite size effects start to become dominant.

We now derive an approximate expression for $`N_{\text{cov}}`$ which is consistent with Eq. (3) and also in good agreement with the large-$`x`$ simulation results. We note that, on average, there are $`2x`$ shortcut ends in the network, defining $`2x`$ contiguous “regions” in the 1-$`d`$ network. If a walker takes a shortcut (or “jump”), it transports the walker from one region to another, randomly chosen, region in the network. When the number of jumps $`n_j`$ is small compared to $`2x`$, there is a high probability that the regions visited by the walker are “distinct” (as defined more precisely below). However, when $`n_j2x`$, there is a high probability that some regions will be visited more than once. These multiple visits lead to saturation effects since the already-visited sites do not contribute to $`N_{\text{cov}}`$.

We now derive an expression for the mean number of distinct regions visited by the walker in $`n_j`$ jumps. Because each jump transports the walker from one region to another, randomly chosen region, the number of distinct regions is obtained by solving the following problem: If $`n_j`$ integers are chosen independently and randomly from the set $`\{1,2,\mathrm{},2x\}`$, what is the mean number of distinct integers chosen? This latter problem is easily solved: The probability that a given integer is not chosen in a given trial is $`q=1\frac{1}{2x}`$; hence, the probability that a given integer is not chosen in $`n_j`$ trials is $`q^{n_j}e^{n_j/(2x)}`$ (where we have assumed $`x1`$). Thus, the mean number of distinct integers chosen after $`n_j`$ trials is

$$I_{\text{cov}}=2x\left(1e^{\frac{n_j}{2x}}\right).$$
(4)
For the original problem, the above expression gives us the number of distinct regions visited by the walker after $`n_j`$ jumps, and each distinct region visited by the walker corresponds to covering $`L/(2x)`$ sites. Furthermore, after time $`t`$, the random walker has, on average, taken $`n_j=4p^2t`$ jumps. Hence, the probability that the walker visits a new region after the $`n_j^{\text{th}}`$ jump is $`\mathrm{exp}(2y/x)`$. Also, the time spent by the walker in this new region is $`t^{}=t\frac{1}{4p^2}4p^2t`$, while the number of sites covered in this time interval is given by $`\sqrt{t^{}}`$, i.e. is proportional to the expression for the $`p=0`$ case. Combining these estimates, we find that the mean number of sites visited after time $`t`$ is given by

$$S(y,x)1e^{\frac{2}{x}y}+\frac{2}{x}\sqrt{yy}e^{\frac{2}{x}y}$$
(5)
Expression (5) is clearly only an approximation to the exact function $`S(y,x)`$. However, it does captures the key processes leading to the growth and saturation of $`N_{\text{cov}}`$ and hence is a useful approximation. In Fig. 3, we compare this function to the numerically obtained curves for $`S(y,x)`$ at $`p=0.01`$, taking $`L=10^5`$ and $`L\mathrm{}`$. The agreement is quite good and shows that Eq. (5) is a reasonable approximation to the scaling function.

### III.3 Mean-Square Displacement

We now compute the mean-square displacement, $`r^2(t)`$, of the random walker as a function of time $`t`$. To calculate this quantity, we first, at each time step, find the minimal distance from the current position of the random walker to the origin (i.e. the smallest number of steps needed for the random walker to reach the origin) using a breadth-first search method. Then we allow the walker to move through the network until $`r^2(t)`$ has saturated. Finally, we average over different initial positions of the walker and realizations of the network.

Now, we know that for a random walk on an infinite, hypercubic, $`d`$-dimensional lattice $`r^2(t)=(1/2)Dt`$, as can be shown using e.g. a generating function formalism barber70 ; montroll64 . However, on a finite lattice, $`r^2(t)`$, must approach a constant for large $`t`$. In this limit, each node of the network has equal probability of being occupied by the random walker. This insight immediately gives $`r^2=\overline{\mathrm{}^2}`$, where $`\overline{\mathrm{}^2}`$ is the squared minimum distance between a pair of nodes, averaged over all possible pairs and network realizations.

On an SWN, the other relevant length scale for the random walker is $`\xi =1/p`$, the average distance the walker travels to reach a shortcut. These two lengths suggest the following scaling Ansatz: $`r^2=\overline{\mathrm{}^2}R(t/\xi ^2;pL)`$. We can also infer the behavior of $`R(y,x)`$ using simple arguments. For times $`t\xi ^2`$, the walker is exploring regions of the network without shortcuts, and we expect diffusive behavior similar to that of a regular network, giving $`R(y,x)y/\overline{\mathrm{}^2}`$. When $`t\sqrt{\overline{\mathrm{}^2}}/p`$, we expect the mean-square displacement to saturate, and $`R(y,x)=1`$. The transition between the two types of behavior is not sharp, since the walker may reach a shortcut before it has travelled a distance $`\xi `$.

We have numerically confirmed this scaling collapse for a wide range of $`x`$ values. In Fig. 4, we plot $`R(t/\xi ^2;x)`$ for a sequence of networks with $`x=pL=10`$, $`10^3`$ and $`10^5`$. For each constant-$`x`$ curve, we show use two distinct values of $`p0.01`$. We see that initially, the walker probes a regular network, and $`R(t/\xi ^2;x)`$ is linear in $`t`$. When $`t/\xi ^21`$, the walker begins to reach some shortcuts in the network, and there is a crossover to superdiffusive behavior. Finally, at still larger $`t/\xi ^2`$, finite-size effects become important and $`R(t/\xi ^2;x)`$ saturates.

To further explore finite-size effects on the mean-square displacement, we have studied $`r^2(t)`$ in a succession of networks, each with the same density of shortcuts, $`p=0.01`$, but with different linear size $`L`$. In Fig. 5, we plot the calculated $`r^2(t)`$ for values of $`L`$ differing by factors of $`10`$ and ranging from $`10^3`$ to $`10^7`$. The resulting curves are all very similar to the $`L\mathrm{}`$ curve, until finite size effects produce saturation with $`r^2=\overline{\mathrm{}^2}`$. Only for the largest $`x`$-values are we able to reach the superdiffusive regime, and even for $`x=10^5`$, this regime is so narrow that we cannot with confidence determine the exponent of the expected power-law time dependence.

The development of an approximate analytical expression for $`r^2(t)`$ is difficult, since we need the minimal distance between two lattice points. However, we know the limiting forms for $`r^2(t)`$. Further, by making the approximation that the random walker only uses a short-cut once, we can extend the arguments used to derive the approximate expression for $`N_{\text{cov}}`$, Eq. (5). Hence, we can write down the following ansatz for $`r^2(t)`$:

$$R(y,x)=1e^{\frac{1}{p^2\overline{\mathrm{}^2}}y^2}\left(1\frac{1}{p^2\overline{\mathrm{}^2}}\left(yy\right)\right).$$
(6)
In Fig. 6, we compare this ansatz with the simulations for $`L=10^7`$; evidently it agrees reasonably well with the numerical results.

We can write down the scaling function in several limiting regimes, without necessarily using the above ansatz for the specific functional form. For example, we must have

$$R(y;x)y/(p^2\overline{\mathrm{}^2}),y1,$$
(7)
and

$$R(y;x)=1,yp\sqrt{\overline{\mathrm{}^2}}.$$
(8)
However, the known limiting behaviors cannot give the scaling function in the intermediate regime. We therefore use the above ansatz to propose

$$R(y;x)y^2/(p^2\overline{\mathrm{}^2}),1yp\sqrt{\overline{\mathrm{}^2}}.$$
(9)
We are not able to verify this last behavior numerically from our present results.

### III.4 First Return Time Distribution

Next, we turn to another property of a walker on an SWN, namely the distribution of first passage (or first return) times. This is the probability that a walker will return to a given site $`m`$ for the first time a time $`t`$ after leaving that site. We denote this distribution $`f(m,t;L,p)`$. Note, that this distribution does not saturate to a finite value when $`t\mathrm{}`$, and hence, it shows that the scaling collapse of random walk properties is not only limited to those which saturate to a finite value. In order to find a smooth distribution, we must make a slight change to the random walker rules: For each time-step, the walker is now allowed to stay on the current node with probability $`1/(k(m)+1)`$.

An extensive review of first passage processes has been given by Redner redner . Unlike the other quantities presented in this article, the first-return time does not saturate to a finite value, instead it approaches zero. However, the scaling collapse is still present. This collapse is demonstrated in Fig. 7, where we plot $`f(0,t)`$ for two different $`(p,L)`$ combinations, holding $`x`$ constant.

The “small-world” effect on the first-return time distribution is most clearly seen at large values of $`x`$. In Fig. 8, we plot $`f(0,t)`$ for a network of $`L=10^6`$ sites and a shortcut density of $`p=0.01`$. For short waiting times, as expected, $`f(0,t)`$ behaves like that of a $`p=0`$ network (no shortcuts). However, for intermediate waiting times such that $`\overline{\mathrm{}^2}tL\xi `$, we find that $`f(0,t)`$ is independent of $`t`$. This independence indicates that the position of the random walker is completely randomized: the memory of the walker’s starting position is no longer retained in the system and the walker is equally likely to occupy any site in the network. This occurs when $`t\overline{\mathrm{}^2}`$, as discussed above for $`r^2(t)`$. Also, $`L\xi `$ is the characteristic time for the asymptotic decay of $`f(m,t;L,p)`$, since it is the saturation timescale for $`N_{\text{cov}}`$. Note that the separation of $`\overline{\mathrm{}^2}(L,p)`$ and $`L\xi `$ increases with increasing number of shortcuts $`x`$ in the system, making the “knee” of the first-return distribution more pronounced.

For a random walk on a finite interval $`[0,L]`$ with an absorbing wall at zero and a reflecting wall at $`L`$ there is a similar effect in the first return timeredner ; redner03 ; at $`tL^2`$ there is an enhancement of the first return probability. In our case, the origin of this effect is the splitting of the timescales $`\overline{\mathrm{}^2}`$ and $`L\xi `$, while for the walk on the finite interval the cause for this effect are the contributions from reflected trajectories.

## IV Summary

In summary, we have studied the behavior of a random walker on a small-world network, using a combination of numerical methods and scaling assumptions. We conjecture that the scaling law of Eq. (1) is obeyed by measurable, saturating properties $`𝒪(p,L,t)`$ of a random walk on an SWN. Among these properties are the average number of distinct visited sites, $`N_{\text{cov}}`$, and the mean-square displacement $`r^2`$ for random walks on SWN’s, both of which we have studied over a wide range of node numbers $`L`$ and shortcut densities $`p`$. In both cases, we find that Eq. (1) is satisfied, and the quantities depend only on the single variable $`x=pL`$. Additionally, we find that non-saturating properties also show a scaling collapse, as exemplified by the first-return time.

Thus, we have shown that the dynamical behavior of a random walker on an SWN has the same scaling behavior as that exhibited by purely geometrical properties of the network (as described in e.g. Ref. kulkarni ). This scaling behavior should be useful in interpreting a variety of other properties on SWN’s, and may be of value in studying real-world phenomena for which an SWN is a good model.

###### Acknowledgements.

This work has been supported by NSF through Grant No. DMR01-04987 (E.A. and D.S.) and the U.S. Department of Energy, Office of Science, Division of Materials Research (R.V.K).