text:
# Binary and Multiple Star Formation in Magnetic Clouds: Bar Growth and Fragmentation

## 1 Introduction

Most main-sequence stars are found in binary or low-order multiple-star systems (Duquennoy & Mayor, 1991; Fischer & Marcy, 1992). Recent observations of star-forming regions have revealed that the binary frequency of pre-main sequence stars is as high as or even higher than that of the main-sequence stars (Leinert et al., 1993; Ghez et al., 1997; Kohler et al., 2000; Mathieu et al., 2000). For example, in the well-studied nearby low-mass star forming Taurus-Auriga region, the binary frequency of pre-main sequence stars is a factor of two higher than that of the main-sequence stars. Other star-forming regions such as Ophiucus, Lupus, Chamaeleon, Scorpius-Centaurus, and Corona Australis show a similar excess of binary fraction. Observations of binary stars in young star clusters such as Orion and Pleiades indicate that the binary frequency of the pre-main sequence stars is comparable to that of the main-sequence stars (Prosser et al., 1994; Padgett et al., 1997; Bouvier et al., 1997). These observations suggest that binary formation is the primary mode of star formation and that binaries should form prior to the pre-main sequence phase, most likely through fragmentation during the gravitational collapse of the parent dense molecular cloud cores (Mathieu, 1994; Bodenheimer et al., 2000; White & Ghez, 2000). Recent detections of binaries and multiple systems at the earliest observable phases of star formation (e.g., Class 0 and Class I) support the fragmentation scenario (e.g., Wootten, 1989; Fuller, Ladd, & Hodapp, 1996; Looney et al., 1997, 2000; Terebey & Padgett, 1997). How the fragmentation occurs is a long-standing, unresolved problem.

Many theoretical studies of star formation have concentrated on the formation of single stars (Shu, Adams, & Lizano, 1987). In the standard picture of single star formation, a strongly magnetized dense molecular cloud core is taken as the initial state of star formation. The dense core gradually contracts as the magnetic support weakens in the high density part due to the ambipolar diffusion. When the magnetic field becomes too weak, it is overwhelmed by self-gravity, and the cloud begins a dynamic contraction to form a star. The process of magnetic leakage during the cloud contraction has been extensively investigated both analytically (e.g., Mestel & Spitzer, 1956; Nakano & Tademaru, 1972) and numerically (e.g., Nakano, 1979; Black & Scott, 1982; Lizano & Shu, 1989; Ciolek & Mouschovias, 1993; Basu & Mouschovias, 1994; Li, 2001). Specific models based on this scenario appear able to explain the observations of well-studied dense cores B1 and L1544 (Crutcher et al., 1994; Ciolek & Basu, 2000; Li et al., 2002). However, these studies have assumed axisymmetry of the model cloud, and are unable to address the issue of fragmentation. On the other hand, numerical hydrodynamic calculations of cloud fragmentation have neglected the effects of magnetic fields (e.g., Tohline, 2002), which are potentially important or even critical.

To quantify the role of magnetic fields on cloud fragmentation, we have made a start in investigating the non-axisymmetric evolution and fragmentation of initially magnetically subcritical clouds, taking into account the effect of magnetic diffusion. Our ultimate goal is to elucidate the fragmentation process in magnetically supported clouds, and extend the current standard scenario of low-mass star formation to binaries and multiple systems.

Magnetically subcritical clouds are stable to both dynamical collapse and fragmentation (e.g., Tomisaka, Ikeuchi, & Nakamura, 1988). Obviously, it is difficult to form binary or multiple stars by direct fragmentation of such clouds. Three-dimensional magnetohydrodynamic (MHD) simulations of magnetic cloud collapse have shown that a frozen-in magnetic field stifles cloud fragmentation (Dorfi, 1982; Benz, 1984; Phillips, 1986a, b). In lightly ionized molecular clouds, however, ambipolar diffusion can initiate dynamical contraction in magnetically supercritical cores by redistributing the magnetic flux in the cloud. The supercritical cores can be unstable to fragmentation if they have masses well above the Jeans limit (e.g., Mestel & Spitzer, 1956; Galli et al., 2001).

Recently, Boss (2000, 2002) carried out three-dimensional collapse calculations of initially magnetically supported clouds, taking into account several properties of magnetic fields approximately. He concluded that the cloud fragmentation is enhanced by magnetic fields because the magnetic tension allows more time for the perturbations to grow off center nonlinearly, which presents the formation of a single, dominant object at the center. Li (2001) studied ambipolar diffusion-driven collapse of magnetized cloud cores with one-dimensional axisymmetric calculations but treating the magnetic forces and ambipolar diffusion more accurately. He also found a similar beneficial effect of magnetic tension. He classified the evolution of magnetically subcritical axisymmetric clouds into two types, depending mainly on the initial cloud mass and initial density distribution. When the initial cloud has a relatively low mass and/or a centrally-condensed density distribution, it collapses to form a single supercritical core (core-forming cloud). On the other hand, when the initial cloud is more massive (containing several Jeans masses or more) and has a relatively flat density distribution, it collapses to form a ring after the central region becomes magnetically supercritical (ring-forming cloud). One-dimensional axisymmetric evolution of the core-forming clouds has also been extensively investigated by Ciolek & Mouschovias (1993) and Basu & Mouschovias (1994). Li & Nakamura (2002) examined the ring-forming clouds and showed that they can fragment into multiple cores in the presence of non-axisymmetric density fluctuations during the isothermal collapse phase after the supercritical core formation. This type of ring fragmentation produces fragments of relatively large initial separations ($`10^4`$ AU). It may be responsible for the formation of small stellar groups such as those in the Serpens cloud core (Williams & Myers, 2000) and aggregates of pre-stellar cores such as those in the $`\rho `$ Oph B2 core (Motte, Andre, & Neri, 1998) and around IRAS 03256+3055 in Perseus (Young et al. 2002). These studies illustrate the crucial role of magnetic diffusion in the fragmentation of magnetic clouds, with strong implications for binary formation.

In this paper, we focus on the non-axisymmetric evolution of the core-forming clouds, with an eye of the formation of binary stars. Initial results are presented in Nakamura & Li (2002), where we demonstrated that the core-forming clouds become unstable to a bar mode after the supercritical core formation \[see also Nakamura & Hanawa (1997) for initially magnetically supercritical clouds\]. Here, we extend the calculations to a larger parameter regime. First, we briefly describe our numerical method and model in §2. In §3, we examine the dependence of bar growth on the initial parameters such as magnetic strength, initial cloud mass, density profile, and rotation. We show that the core-forming clouds become unstable to a bar mode during the dynamic collapse phase over a wide range of initial model parameters. Subsequent fragmentation process is investigated in §4. Based on the numerical results, we classify the fragmentation process into three types, which agree with the empirical classification of young binary and multiple systems by Looney et al. (2000) based on high-resolution millimeter observations. This general agreement and other implications of our results on binary formation are discussed in §5.

## 2 Formulation of The Problem

A detailed description of the model cloud and the problem formulation are presented in §2 of Li & Nakamura (2002). Here, we give a brief summary. We consider an initially axisymmetric, magnetically subcritical model cloud that is in a mechanical equilibrium, with a disk-like mass distribution due to settling along ordered magnetic field lines. In such a subcritical cloud, force balance along the field lines is achieved during the quasi-static phase of evolution and is well maintained even during the subsequent dynamic collapse (e.g., Fiedler & Mouschovias, 1993; Nakamura, Hanawa, & Nakano, 1995). Thus, we adopt the standard thin-disk approximation in which exact force balance along the vertical direction is assumed during the entire evolution (e.g., Ciolek & Mouschovias, 1993; Basu & Mouschovias, 1994; Shu & Li, 1997; Nakamura & Hanawa, 1997; Stehle & Spruit, 2001; Li, 2001). The magnetohydrodynamical equations are integrated in the vertical direction and the cloud evolution is followed in the ($`x`$-$`y`$) disk plane of a Cartesian coordinate system $`(x,y,z)`$ with a two-dimensional MHD code. The distribution of magnetic field is computed in three-dimensional space under the assumption that the magnetic field is current-free outside the disk.

We assume an isothermal equation of state below a dimensionless surface density $`\mathrm{\Sigma }_{\mathrm{cr}}`$, and denote the isothermal sound speed by $`c_s`$. Beyond $`\mathrm{\Sigma }_{\mathrm{cr}}`$, we change the equation of state from isothermal to adiabatic, with an index of $`5/3`$, to mimic the transition to optically thick regime of cloud evolution. In this paper, we set $`\mathrm{\Sigma }_{\mathrm{cr}}=1.9\times 10^4\mathrm{\Sigma }_{0,\mathrm{ref}}`$, corresponding to the volume density $`\rho _{\mathrm{cr}}10^{12}`$ cm-3 (see eq.  below for the definition of $`\mathrm{\Sigma }_{0,\mathrm{ref}}`$).

The initial conditions for star formation are not well determined either observationally or theoretically. Following Basu & Mouschovias (1994) and Li (2001), we prescribe a slowly-rotating, axisymmetric reference state for our model cloud, with the distributions of mass, magnetic field, and velocity given by

$`\mathrm{\Sigma }_{\mathrm{ref}}(x,y)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Sigma }_{0,\mathrm{ref}}}{\left[1+(r/r_0)^n\right]^{4/n}}},`$ (1)
$`B_{z,\mathrm{ref}}(x,y)`$ $`=`$ $`B_{z,\mathrm{}}=\mathrm{const}.,`$ (2)
$`V_\varphi (x,y)`$ $`=`$ $`{\displaystyle \frac{4\omega r}{r_0+\sqrt{r_0^2+r^2}}}c_{ms},`$ (3)
where $`r=(x^2+y^2)^{1/2}`$ is the distance from the origin, $`c_{ms}=c_s(1+\mathrm{\Gamma }_0^2)^{1/2}`$ is essentially the magnetosonic speed, and the parameter $`\omega `$ measures the rotation rate. The exponent $`n`$ controls the surface density profile and (indirectly) the amount of mass in the central plateau region where the distribution is more or less uniform. The background field strength $`B_{z,\mathrm{}}`$ is characterized by the dimensionless flux-to-mass ratio $`\mathrm{\Gamma }_0=B_{z,\mathrm{}}/(2\pi G^{1/2}\mathrm{\Sigma }_{0,\mathrm{ref}})`$. The reference clouds are not in an equilibrium state and are allowed to evolve into one with magnetic field frozen-in, before ambipolar diffusion is turned on at time $`t=0`$. To save computational time, we obtain the equilibrium state using the one-dimensional axisymmetric MHD code of Li (2001), which gives results indistinguishable from those from the two-dimensional code assuming axisymmetry.

At $`t=0`$, we impose on top of the equilibrium surface density $`\mathrm{\Sigma }_0(x,y)`$ an $`m=2`$ perturbation of relative amplitude $`A`$,

$$\mathrm{\Sigma }(x,y)=\mathrm{\Sigma }_0(x,y)[1+A\mathrm{cos}(2\varphi )],$$
(4)
where $`\varphi `$ is the azimuthal angle measured from the $`x`$-axis. The subsequent, ambipolar diffusion-driven evolution of the perturbed cloud is followed numerically with the two dimensional MHD code. As a boundary condition, we fixed $`\rho `$ and $`B_z`$ at the cloud outer radius, taken to be twice the characteristic radius $`r_0`$.

The computations are carried out using non-dimensional quantities. The units we adopted are $`c_s`$ for speed, $`\mathrm{\Sigma }_{0,\mathrm{ref}}`$ for surface density, $`2\pi G\mathrm{\Sigma }_{0,\mathrm{ref}}`$ for gravitational acceleration, and $`B_{z,\mathrm{}}`$ for magnetic field strength. The units for length and time are, respectively, $`L_0c_s^2/(2\pi G\mathrm{\Sigma }_{0,\mathrm{ref}})`$ and $`t_0c_s/(2\pi G\mathrm{\Sigma }_{0,\mathrm{ref}})`$. They have typical values of

$`L_0`$ $`=`$ $`8.43\times 10^2\mathrm{pc}\left({\displaystyle \frac{c_s}{0.33\mathrm{km}\mathrm{s}^1}}\right)^2\left({\displaystyle \frac{\mathrm{\Sigma }_{0,\mathrm{ref}}}{10^2\mathrm{g}\mathrm{cm}^2}}\right)^1,`$ (5)
$`t_0`$ $`=`$ $`2.50\times 10^5\mathrm{yr}\left({\displaystyle \frac{c_s}{0.33\mathrm{km}\mathrm{s}^1}}\right)\left({\displaystyle \frac{\mathrm{\Sigma }_{0,\mathrm{ref}}}{10^2\mathrm{g}\mathrm{cm}^2}}\right)^1,`$ (6)
where $`c_s=0.33`$ km s-1 corresponds to an effective temperature of $`T_{\mathrm{eff}}=30`$ K.

The initial state of our model cloud is completely specified by four parameters: the characteristic radius $`r_0`$, exponent $`n`$ in the reference surface density distribution, dimensionless flux-to-mass ratio $`\mathrm{\Gamma }_0`$, and rotation parameter $`\omega `$. We adopt a relative small rotation parameter of $`\omega <0.20.3`$, consistent with the slow rates of rotation observed in molecular cloud cores (Goodman et al., 1993). The slow rotation is thought to be due to magnetic braking by an ambient medium (Nakano & Tademaru, 1972; Basu & Mouschovias, 1994), which is not treated in our paper. We comment briefly on the potential effects of magnetic braking in § 5.3.

Besides the four parameters characterizing the initial cloud, a parameter $`A`$ is needed to specify the fractional amplitude of the $`m=2`$ bar mode perturbation applied after the reference cloud settles into the equilibrium configuration. We have carried out a systematic parameter survey and selected the models listed in Table 1 to illustrate the salient features of magnetic cloud evolution. The numerical results are described in the next two sections (§3 and §4).

## 3 Bar Growth in Magnetically Subcritical Clouds

As mentioned earlier, Li (2001) classified the evolution of magnetically subcritical axisymmetric clouds into two types, depending mainly on the initial cloud mass and initial density distribution: (a) core-forming cloud, which collapses to form a single supercritical core and (b) ring-forming cloud in which a supercritical ring is formed after the central region becomes supercritical owing to ambipolar diffusion. The ring fragments into multiple supercritical cores in the presence of density fluctuations during the isothermal collapse phase (Li & Nakamura, 2002). In this paper, we focus on the core-forming clouds and show that these clouds become unstable to bar formation after the supercritical core formation over a wide range of initial model parameters. The subsequent evolution and fragmentation of the bar are discussed in §4.

### 3.1 Non-Rotating Cloud

To begin with, we consider a non-rotating ($`\omega =0`$) cloud with the reference state specified by the set of parameters $`r_0=10\pi `$, $`n=2`$, and $`\mathrm{\Gamma }_0=1.5`$. The cloud is therefore magnetically subcritical, with the field strength $`50\%`$ above the critical value at the center. The cloud so specified produces a dense supercritical core in the absence of any non-axisymmetric perturbations. The reference cloud is allowed to evolve into an equilibrium configuration, with the magnetic field frozen-in. After the equilibrium is reached, we reset the time to $`t=0`$, turn on ambipolar diffusion, and add an $`m=2`$ perturbation to the surface density distribution, with a fractional amplitude of $`A=0.05`$. The evolution of such a non-axisymmetrically perturbed cloud is followed numerically to progressively higher densities and smaller scales. The results are shown in Fig. 1 and described as follows.

The evolution of the model is qualitatively similar to that of the model shown in Fig. 1 of Nakamura & Li (2002), where a different reference density distribution (less centrally-condensed than the model shown in this subsection) is adopted. The evolution is followed to a higher density here. As is usual with ambipolar diffusion-driven evolution, the cloud spends most of its time in the subcritical phase when the dimensionless minimum flux-to-mass ratio (at the density peak) is greater than unity, $`\mathrm{\Gamma }_c>1`$ (see panels \[a\] and \[b\] of Fig. 1). During this period, gas motions are very subsonic and the central part of the cloud oscillates with a small amplitude. The oscillation indicates that the cloud is stable to the $`m=2`$ perturbation during the subcritical phase, in agreement with linear analysis by Nakano (1988) and others. Indeed, we have verified numerically that for a frozen-in field the subcritical cloud remains stable to this and other higher (e.g., $`m=3,4,\mathrm{}`$) modes. (We followed the evolution of the frozen-in model for several periods of oscillation and confirmed that the cloud oscillates with almost the same amplitude as the initial.) Once the minimum flux-to-mass ratio has dropped below the critical value, the contraction becomes dynamic and the infall speed eventually exceeds the sound speed at the center (see panels \[c\] and \[d\]). By the time shown in panel (c), the $`m=2`$ mode has grown significantly, resulting in a bar-like core at the center with an aspect ratio of roughly 2. The bar growth remains relatively modest (panel \[d\]) until the very end of the runaway collapse, when the growth rate increases dramatically (panel \[e\]). Beyond a dimensionless critical surface density of $`1.9\times 10^4`$, corresponding to a volume density of $`n_H10^{12}`$ cm-3, we change the equation of state from isothermal to adiabatic, with an index of 5/3, to mimic the transition to optically thick regime of cloud evolution. When the maximum surface density reaches $`10^5`$, the contraction along the minor axis of the bar decelerates significantly and an accretion shock is formed at the surface of the bar. This shock-bound bar is analogous to the “first core” of spherical calculations, and will be referred to as the “first bar” below. Immediately after the shock formation, the aspect ratio of the first bar continues to increase. It reaches a maximum of about 10 by the time shown in panel (e). Thereafter, the infall along the major axis becomes significant, forming a round core at the center. The aspect ratio of the central round core approaches unity as the collapse proceeds. Since the round core resembles the first core of the spherical calculations, we call it the “truly first core”. As will be shown in §4, this round core evolves into a rapidly-rotating disk in the presence of an initial, slow cloud rotation.

To quantify the geometrical change of the supercritical core, we plot in Figure 2 the aspect ratio of the darkest region in Figure 1 (with $`\mathrm{\Sigma }>10^{1/2}\mathrm{\Sigma }_{\mathrm{max}}`$) as a function of the maximum surface density, $`\mathrm{\Sigma }_{\mathrm{max}}`$, together with the flux-to-mass-ratio $`\mathrm{\Gamma }`$ at the center. The aspect ratio is defined as $`R=a_y/a_x`$ where $`a_x`$ and $`a_y`$ are, respectively, the effective scale heights along the $`x`$ and $`y`$ axes, and evaluated as

$$a_x=\left(M^1_{\mathrm{\Sigma }10^{1/2}\mathrm{\Sigma }_{\mathrm{max}}}\mathrm{\Sigma }x^2𝑑x𝑑y\right)^{1/2},$$
(7)
and

$$a_y=\left(M^1_{\mathrm{\Sigma }10^{1/2}\mathrm{\Sigma }_{\mathrm{max}}}\mathrm{\Sigma }y^2𝑑x𝑑y\right)^{1/2},$$
(8)
where $`M`$ is the total mass in the region with $`\mathrm{\Sigma }10^{1/2}\mathrm{\Sigma }_{\mathrm{max}}`$. During the subcritical phase, the aspect ratio of the bar does not grow in time; it oscillates around unity, signifying that the subcritical cloud is stable to the $`m=2`$ mode. As the central region becomes more and more supercritical, the bar mode begins to grow significantly. During this supercritical period, the aspect ratio of the bar stays more or less frozen at $`R1.5`$ before the central density reaches $`\mathrm{\Sigma }10^3`$. This freezing of aspect ratio is related to the dynamical state of the disk. As shown in §3.4, the mass of the central flat region is merely $``$ 1.5 times the effective Jeans mass (as measured by the central surface density, including the magnetic contribution). Such a relatively low-mass disk is only marginally unstable to the non-axisymmetric perturbation. In other words, the time scale of bar growth is comparable to that of dynamic contraction. The flux-to-mass ratio does not change much either during this period, as a result of the rapid collapse which prevents the magnetic flux from leaking out significantly. The bar growth speeds up toward the very end of the supercritical collapse when the mass of the central flat region exceeds about twice the effective Jeans mass. It continues after the first bar formation. When the central surface density reaches $`10^5`$, the infall motion along the major axis of the bar becomes significant and the aspect ratio begins to drop steeply to unity, forming the truly first core.

To study the structure of the bar, we show in Figure 3 the distributions of the surface density, magnetic field, and infall velocity along the axes as a function of the distance from the density maximum, at seven different times. The solid and dashed curves indicate the distributions of the physical quantities along the major ($`y`$) and minor ($`x`$) axes, respectively. During the early phase of supercritical collapse when the aspect ratio of the bar is frozen at $`R1.5`$, the surface density distributions along both axes qualitatively resemble the self-similar distribution of the axisymmetric collapse; they are flat near the density maximum and can be approximated by a power-law of $`\mathrm{\Sigma }d^1`$ at the magnetically supercritical envelope, where $`d`$ is the distance from the density peak. The $`\mathrm{\Sigma }d^1`$ distribution corresponds to a volume density distribution of $`\rho d^2`$ because of the force balance in the vertical direction (which yields $`\rho \mathrm{\Sigma }^2`$). While the central density increases rapidly, the surface density distributions change little in the envelope. The central flat region shrinks in size in accordance with the increase in the central density. The distributions of the vertical field $`B_z`$ are qualitatively similar to those of $`\mathrm{\Sigma }`$ in the magnetically supercritical region and almost proportional to the local surface density ($`B_z\mathrm{\Sigma }`$). The $`x`$\- and $`y`$-components of the magnetic field (in the disk plane) increase approximately linearly with distance $`d`$ in the central flat region, peaking near the edge of the flat region, before falling off as $`B_xd^1`$ and $`B_yd^1`$. In the magnetically supercritical envelope, $`B_x`$ and $`B_y`$ are comparable to $`B_z`$, indicating that the field lines near the disk plane incline with an angle of $``$ 45 degree with respect to the disk plane. The angle approaches $`90`$ degree in the outer, magnetically subcritical envelope. The infall speeds are proportional to the distance $`d`$ near the center, becoming more or less flat in the supercritical part of the envelope, before falling to subsonic values in the outer subcritical envelope. At $`\mathrm{\Sigma }_{\mathrm{max}}=10^4`$, the infall speeds in the supercritical envelope reach $`2c_s`$.

As the collapse proceeds and the bar aspect ratio becomes large, the distributions of physical quantities gradually deviate from the self-similar distributions. For example, the surface density and vertical magnetic field ($`B_z`$) profiles along the major axis become steeper as the bar elongates rapidly. When the central surface density reaches $`2\times 10^4`$, a shock is formed at $`x10^4r_0`$. Thereafter, the infall speed rises in the envelope of the bar. The infall speed along the minor axis follows $`v_xd^{1/2}`$ just outside the first bar, indicating a free-fall collapse towards the bar.

Interestingly, after the bar becomes opaque to dust emission, a local density peak is formed at each end of the first bar. This density peak formation is due to the “edge” effect, which was first calculated by Larson (1972) and then extensively investigated by Bastien (1983) and Bonnell & Bastien (1991). In this model, the density peaks do not grow nonlinearly; they tend to disappear in the subsequent collapse because the central truly first core grows more rapidly in mass owing to the supersonic infall along the major axis, which dominates the gravitational potential in the high density region. In §4, we show that for more massive clouds in which the aspect ratios of the supercritical cores are larger, the density peaks can grow in time because the infall speeds along the major axis remain subsonic.

The model discussed in this subsection serves as a standard against which other models will be compared.

### 3.2 Dependence of Model Parameters: $`\mathrm{\Gamma }_0`$, $`r_0`$, and $`n`$

The properties of the bar growth depend on the initial model parameters. Here, we investigate the effects of the initial flux-to-mass ratio $`\mathrm{\Gamma }_0`$, the characteristic radius $`r_0`$, and the density profile parameter $`n`$ for non-rotating clouds. The effect of rotation is examined in separate subsections §3.3 and §4.1. Exact values of $`\mathrm{\Gamma }_0`$ are difficult to determine observationally; they probably lie within a factor of two of unity (Crutcher, 1999), after correcting for likely projection effects (Shu et al., 1999). Thus, we investigated the dependence on $`\mathrm{\Gamma }_0`$ in the range of $`1<\mathrm{\Gamma }_0<3`$. The characteristic radius $`r_0`$ is to be compared with $`2\pi `$, the (dimensionless) critical wavelength for the (thermal) Jeans instability in a disk of uniform mass distribution (Larson, 1985). The parameter $`n`$ controls the profile of the cloud surface density distribution; those with smaller values of $`n`$ are more centrally peaked.

To examine the effect of the initial flux-to-mass ratio $`\mathrm{\Gamma }_0`$, we first hold the cloud radius and profile parameter fixed at the standard values $`r_0=10\pi `$ and $`n=2`$, and reduce the flux-to-mass ratio $`\mathrm{\Gamma }_0`$. Representative results are shown in panel (a) of Fig. 4, at a time when the maximum surface density $`\mathrm{\Sigma }_{\mathrm{max}}=10^4`$. Panel (a) is the snapshot of a cloud with $`\mathrm{\Gamma }_0=1.25`$ (a weaker initial magnetic field), perturbed by an $`m=2`$ mode of the standard fractional amplitude $`A=0.05`$. Comparison with panel (d) of Fig. 1 where the standard case at the same $`\mathrm{\Sigma }_{\mathrm{max}}`$ is plotted shows that a weaker magnetic field promotes the bar growth. This slower bar growth in the stronger field case is quantitatively shown in panel (a) of Figure 2 (compare the solid curve with the dotted-dashed curve). The reason for this behavior appears to be that, for the cloud with a stronger initial magnetic field, the magnetic support remains more important even in the supercritical phase. In other words, the averaged flux-to-mass ratio in the darkest region ($`\mathrm{\Sigma }10^{1/2}\mathrm{\Sigma }_{\mathrm{max}}`$) is larger for the stronger field model. The stronger field tends to retard the deformation of the supercritical core. We have confirmed this negative effect of the magnetic field on bar growth by following the evolution of the clouds with $`\mathrm{\Gamma }_0=2`$ and 3.111We note that the magnetic field has an opposite, beneficial effect on the fragmentation of ring-forming clouds. More strongly magnetized clouds contract more slowly, allowing more time for non-axisymmetric perturbations to grow off-center. As a result, multiple fragmentation is achieved more readily, at a lower surface density. See Li & Nakamura (2002) for details.

We next consider the effect of the characteristic radius $`r_0`$ which controls the total cloud mass. Increasing the radius $`r_0`$ is expected to have a positive effect on bar growth. This is indeed the case, as illustrated in panel (b) of Fig. 4, where the snapshot of a cloud with $`r_0=15\pi `$, 1.5 times the standard value, is shown at $`\mathrm{\Sigma }_{\mathrm{max}}=10^4`$. The aspect ratio of the bar reaches $`7`$ at the time when the surface density reaches $`10^4`$.

The profile parameter $`n`$ has an effect similar to the parameter $`r_0`$. For instance, a reduction in $`n`$ would decrease the number of Jeans masses contained in the central plateau region. We showed in a previous paper (Li & Nakamura, 2002) that for the ring-forming case, the clouds with more peaked density profiles are less prone to fragmentation. This negative effect of the centrally condensed density distribution is also observed for non-magnetized clouds (Boss, 1993). We find the same trend for the core-forming clouds as well. The aspect ratio of the bar increases more rapidly for larger $`n`$. This is illustrated in panel (c) of Fig. 4, where $`n=4`$, twice the standard value, is chosen. The bar is more elongated than the standard case, even though the characteristic radius $`r_0=5\pi `$ is only half of the standard value.

### 3.3 Slowly-Rotating Cloud

To gauge the effect of slow rotation, we repeat the evolution of a cloud with all parameters identical to those of the standard model in §3.1 except the rotation parameter $`\omega `$, which is now set to $`0.1`$ instead of zero. Snapshots of the surface density distribution and velocity field at three selected times are shown in Fig. 5. The general trend of the cloud evolution is similar to that in the non-rotating case; in the subcritical phase, the cloud is stable to non-axisymmetric perturbations and the $`m=2`$ mode begins to grow in time only after the central region becomes magnetically supercritical. Comparison between Fig. 1d and Fig. 5b indicates a beneficial effect of the slow rotation on the bar growth. More quantitatively, the aspect ratio of the bar is always larger for the rotating cloud at the time when the maximum density reaches the same value, before the bar becomes opaque to dust emission. This effect derives from the fact that rotation retards the radial infall, which allows more time for perturbations to grow.

In §4.1 below, we will show that the slow-down effect due to rotation plays a dramatic role in the cloud fragmentation for more rapidly-rotating, less centrally-condensed clouds, which can fragment into two cores during the isothermal collapse phase soon after the supercritical core formation.

### 3.4 Dynamical State of Supercritical Cores: Axisymmetric Collapse

The bar growth is closely related to the dynamical state of the supercritical core. In this subsection, to clarify the physics of the bar growth, we follow the evolution of unperturbed axisymmetric clouds with the one-dimensional code and measure the dynamical forces in the supercritical cores. For simplicity, we do not change the equation of state for the one-dimensional calculations even when the surface density exceeds the critical value of $`1.9\times 10^4`$. Previous axisymmetric calculations of ambipolar diffusion-driven cloud evolution have shown that the collapse of supercritical cores proceeds self-similarly: the distributions of the surface density and magnetic field are almost constant in the central flat region and decrease with radius as $`\mathrm{\Sigma }r^1`$ and $`B_zr^1`$ (Ciolek & Mouschovias, 1993; Basu & Mouschovias, 1994). In the high density region of the supercritical core, the decrease of the flux-to-mass ratio levels off (at $`\mathrm{\Gamma }0.4`$) due to the rapid dynamic collapse which prevents significant magnetic flux leakage. This is shown in panels (a) and (c) of Figure 6, where the flux-to-mass ratio normalized by the critical value is depicted as a function of radius for two representative models. Shu & Li (1997) and Li & Shu (1997) called such a cloud with a constant flux-to-mass ratio “isopedic”. In the isopedic cloud, the gravity ($`F_{\mathrm{grav}}=\mathrm{\Sigma }g_r`$) and gas pressure ($`F_{\mathrm{gas}}=c_s^2d\mathrm{\Sigma }/dr`$) are proportional to the magnetic tension ($`F_{\mathrm{ten}}=B_rB_z/2\pi `$) and pressure ($`F_\mathrm{B}=Hd(B_z^2/4\pi )/dr`$), respectively (Shu & Li, 1997; Li & Shu, 1997; Nakamura & Hanawa, 1997), where $`H`$ is the half thickness of the disk (see eq.  of Li & Nakamura (2002) for the definition of $`H`$). In fact, the supercritical cores have such a characteristic.

Figure 6 shows the evolution of the ratio of the effective gravity to the effective pressure at six different times for two models. Here, the effective pressure and gravity are defined as $`F_{\mathrm{gas}}+F_\mathrm{B}`$ and $`F_{\mathrm{grav}}+F_{\mathrm{ten}}`$, respectively. At the initial state ($`t=0`$), this ratio is equal to unity, indicating a force balance in the initial cloud. As the collapse proceeds, this ratio increases due to ambipolar diffusion, which reduces the magnetic support in the central high density region. When the central surface density is in the range of $`10^310^4`$, the effective gravity is larger than the effective pressure, but only by a factor of about 1.5 $``$ 2. This corresponds to a mass of the central flat region only about 1.5 $``$ 2 times the effective Jeans mass. Such a relatively low-mass core is only marginally unstable to non-axisymmetric perturbations, which is probably the reason why the aspect ratio of the bar remains almost frozen at a modest value of $`1.52`$ in the presence of an $`m=2`$ perturbation. When the force ratio exceeds $`2`$, the bar growth is accelerated remarkably by the Lin-Mestel-Shu type gravitational instability (Lin, Mestel, & Shu, 1965), as seen in Fig. 2. For both clouds, the force ratios approach a higher value of $`34`$ as the collapse proceeds, although a convergence hasn’t been reached by the time the central surface density reaches the critical value of $`1.9\times 10^4`$. (The convergence toward a self-similar solution is slower for more centrally-condensed clouds.) This behavior suggests that the isothermal evolution tends toward a self-similar solution.

A self-similar solution reproducing the pre-protostellar collapse of an isopedic magnetic disk is derived approximately by Nakamura & Hanawa (1997), who considered a frozen-in magnetic field \[see also Saigo & Hanawa (1998) for the nonmagnetized case\]. Basu (1997) constructed a self-similar collapse model including the effect of ambipolar diffusion \[see Li (1998) for the self-similar collapse after the protostellar core formation\]. The behavior of his solution is in good agreement with that of Nakamura & Hanawa (1997) in the late, self-similar phase of dynamic collapse, when the effective gravity is about $`34`$ times the effective pressure force and the asymptotic velocity is about twice the effective sound speed. Nakamura & Hanawa (1997) investigated the dynamical stability of the self-similar solution to non-axisymmetric perturbations and found that the self-similar solution is unstable to the $`m=2`$ mode and stable to the higher modes. This characteristic is consistent with the numerical results presented in the previous subsections. Thus, we conclude that the tendency for the supercritical collapse to approach the gravity-dominated self-similar solution is responsible for the bar formation during the dynamic collapse.

## 4 Fragmentation of Supercritical Cores

In §3, we demonstrated that initially magnetically subcritical clouds tend to be unstable to an $`m=2`$ mode, forming bar-like cores. In this section, we investigate the subsequent evolution of the bar mode and show that the elongated supercritical cores can break up into discrete fragments. We classify the fragmentation process into three types, depending on the initial model parameters: (1) separate core formation, in which the $`m=2`$ mode breaks up into two separate cores in the isothermal collapse phase, (2) bar fragmentation, in which the first bar breaks up into multiple fragments that move mainly along the major axis of the bar, with frequent merging among the fragments, and (3) disk fragmentation, in which the central region evolves into a rapidly-rotating disk where another bar mode is excited, fragmenting into blobs. In the following, we select three models to illustrate the salient features of these cases.

### 4.1 Separate Core Formation: Effect of Faster Rotation

In §3.3, we showed that rotation slows down the dynamic collapse, allowing more time for perturbations to grow. Here, we demonstrate that a relatively fast rotation can affect the outcome of cloud evolution dramatically, particularly for $`n>4`$ clouds whose initial parameter range is intermediate between those for core-forming and ring-forming clouds in the axisymmetric calculations (Li, 2001). Panels (b) and (c) of Fig. 7 show the snapshots of the surface density distribution and velocity field at two selected times for the model with $`n=4`$, $`\omega =0.2`$, $`r_0=5\pi `$, and $`\mathrm{\Gamma }_0=1.5`$. This cloud is categorized as a core-forming cloud from the axisymmetric one-dimensional calculation. A cloud larger in radius by only a few tens of percent would collapse to form a ring instead of a single core. For comparison, the snapshot of a more slowly rotating cloud with $`\omega =0.1`$ is shown in panel (a) at the time when $`\mathrm{\Sigma }_{\mathrm{max}}=10^2`$. The model with $`\omega =0.1`$ shows the evolution qualitatively similar to the model in §3.3, i.e., after the supercritical core formation, a bar mode grows in time without fragmentation. On the other hand, when the rotation rate is doubled ($`\omega =0.2`$), the bar grows more rapidly even at the early time shown in panel (b), when the maximum surface density reaches 10. By this time, the entire bar (the darkest part in panel \[b\]) becomes magnetically supercritical. By the time shown in panel (c) ($`\mathrm{\Sigma }_{\mathrm{max}}=10^2`$), the $`m=2`$ mode breaks up into two separate cores, each of which continues to collapse dynamically. (We checked that this model is stable to higher $`m`$ modes.) The subsequent evolution of each core is qualitatively similar to that of the models shown in the previous section, i.e., the aspect ratio of each core is more or less frozen at $``$ 2 during the early phase of supercritical collapse, and increases rapidly at the very end of the supercritical collapse.

For the separate core formation case, the magnetic field has a positive effect on the fragmentation. We calculated the models of a weaker and a stronger magnetic field for the same set of parameters $`n=4`$, $`r_0=5\pi `$, $`\omega =0.2`$, and $`A=0.05`$, and found that the cloud with a weaker magnetic field does not break up into two cores (it collapses into a single bar instead), whereas the cloud with a stronger magnetic field fragments into two cores at a lower density and with a wider separation. More centrally-condensed clouds with $`n=2`$ and $`r_0<20\pi `$ do not fragment into two cores even in the presence of a rapid rotation ($`\omega 0.3`$). These dependences on $`\mathrm{\Gamma }_0`$ and $`n`$ are the same as those of the multiple fragmentation case investigated by Li & Nakamura (2002). In this sense, the separate core case is similar to the multiple fragmentation case, although the bar fragmentation into two cores appears to depend more strongly on the rotation rate.

At the time shown in panel (c), the separation of two cores is about a tenth of the initial cloud radius, which corresponds to $`5\times 10^4`$ AU with the fiducial choice of $`\mathrm{\Sigma }_0=0.01`$ g cm-2 and $`T_{\mathrm{eff}}=30`$ K. This separate core formation case may thus be responsible for the formation of wide binaries with separations of the order of $`10^4`$ AU, although substantial orbital evolution is likely after the initial fragmentation, since the cores are not kept apart by rotation.

To summarize, when the initial rotation is relatively fast and the density distribution in the central region is relatively flat, the $`m=2`$ mode can break up into two separate cores in the isothermal collapse phase after the central region becomes supercritical (e.g., for $`n=4`$, $`\omega `$ 0.2, $`\mathrm{\Gamma }_01.5`$, and $`r_0>5\pi `$). We note that this separate core formation differs from that found previously for non-magnetic clouds, which requires a rapid rotation or highly unstable initial condition (e.g., Tsuribe & Inutsuka, 1999). In our case, the cloud is initially stable to gravitational collapse, and the fragmentation is induced by ambipolar diffusion. Indeed, we find that clouds with relatively weak magnetic fields do not fragment into separate cores during the isothermal collapse phase.

### 4.2 Bar Fragmentation and Merging of Fragments

According to linear perturbation theory, long filamentary clouds tend to break up gravitationally into fragments if their lengths exceed several times the diameter. This is indeed the case for the first bars formed in relatively massive, centrally-condensed clouds (with a relatively large $`r_0`$ and $`n<4`$, e.g., for $`n=2`$, $`\mathrm{\Gamma }_01.5`$, $`\omega 0.2`$, and $`r_0>15\pi `$), as illustrated in Fig. 8, where the snapshots of a cloud with $`n=2`$, $`\mathrm{\Gamma }_0=1.5`$, $`r_0=15\pi `$, and $`\omega =0.2`$ are shown at three different times. In this model, fragmentation starts when the bar becomes opaque to dust emission, which retards the radial contraction (at $`\mathrm{\Sigma }_c>10^5`$). At first, the bar produces two dense fragments near the center at $`\mathrm{\Sigma }10^5`$, and then the fragmentation apparently propagates towards the ends of the bar. This apparent propagation of the fragmentation is due to the fact that the growth rate of the gravitational fragmentation is larger for the higher density region closer to the center. By the time shown in panel (a), there is a fragment formed at each of the two ends of the bar. This density peak formation at the end was first noticed by Larson (1972) and subsequently investigated extensively by Bastien (1983) and Bonnell & Bastien (1991). It is caused by a local minimum in the gravitational potential. The masses of the fragments formed by bar fragmentation are on the order of $`10^2M_{}`$. These masses agree with those expected from the linear theory of fragmentation of cylindrical clouds (e.g., Larson, 1985; Nakamura, Hanawa, & Nakano, 1993). By the time shown in panel (b), the two pairs of fragments near the center have merged to produce one, more massive, pair. This process of successive bar fragmentation and subsequent merging is repeated, as shown in panel (c), where the two most massive fragments are about to merge at the center. We will discuss further the evolution of this multiple fragment system in §5.

We note that bar fragmentation was also observed in the numerical simulations of non-magnetic cloud collapse (e.g., Boss, 1993; Bate & Burkert, 1997). The non-magnetic bars are typically much shorter than the one presented in this subsection. To form such a long bar in the non-magnetic case, the initial cloud must contain many thermal Jeans masses and thus be highly unstable. In our case, the initial cloud is stablized by a strong magnetic field, whose support weakens with time as a result of ambipolar diffusion. Nakamura & Hanawa (1997) studied bar formation in clouds with frozen-in magnetic fields. They found bars that are prone to gravitational breakup according to the linear theory of fragmentation. These bars are however shorter than the one discussed here. They are expected to produce a smaller number of fragments.

### 4.3 Disk Formation and Fragmentation

When the initial clouds are centrally condensed and/or less massive (e.g., for $`n=2`$, $`\omega 0.1`$, $`\mathrm{\Gamma }_01.5`$, and $`r_0<15\pi `$), the length of the bar does not exceed the critical wavelength for gravitational fragmentation by the time of first bar formation. It would be difficult for such a short bar to break up into blobs gravitationally. In this case, a rapidly-rotating disk is formed in the presence of slow rotation. Here, we show an example of the rotating disk formation. Figure 9 shows the snapshots at six different times of a model with $`n=2`$, $`\mathrm{\Gamma }_0=1.5`$, $`r_0=5\pi `$, and $`\omega =0.1`$, which is identical to the slowly rotating model discussed in §3.3, except for a smaller characteristic radius $`r_0`$ ($`5\pi `$ instead of $`10\pi `$). For comparison, the linear density contours of the central disk is depicted in Fig. 10 at three selected times. In this model, the bar mode does not grow much by the end of isothermal collapse because of the relatively small cloud radius $`r_0`$222 We estimate that the central plateau region of this cloud has a size roughly $`23`$ times the effective Jeans length, judging from the fact that, for $`n=2`$, reference clouds with $`r_0<2\pi `$ expand rather than contract.. Because the length of the bar is smaller than the critical wavelength for gravitational fragmentation, it does not fragment. Instead, a round dense region, which we called the “truly first core”, is formed at the center due to a supersonic infall motion along the bar.

Initially, the truly first core is supported primarily by the thermal pressure. As ambient gas with a larger angular momentum accretes onto it, the rotational support becomes dominant. As a result, the central core evolves into a rapidly-rotating, opaque disk, which we will sometimes refer to as the “first disk”. The interaction between the disk and the large-scale bar mode outside the disk excites another bar-mode growing inside the disk. We believe that the secondary bar growth is due to rotational instability because the ratio of the rotational energy and the gravitational energy of the first disk reaches the critical value of $`\beta _{\mathrm{cr}}0.3`$ (e.g., Durisen et al., 1986; Bate, 1998; Tohline, 2002). By the time shown in panel (c) of Fig. 10, the bar breaks up into two blobs, which appear to be self-gravitating. We stopped our calculation at this epoch of fragmentation because the volume densities of the blobs reach the critical density beyond which H2 dissociation starts and the equation of state of the gas is expected to change again.

The fate of the two blobs formed in the disk is uncertain. They may merge together during the subsequent evolution. Since the rotation period of the disk is faster than that of the large-scale bar outside the disk, the interaction between the two bar modes generates several blobs in the disk, some of which appear to be self-gravitating. These blobs are expected to interact gravitationally and whether any of them can survive the interaction remains to be determined. This type of fragmentation and merging in a rotating disk was also observed in the numerical simulations of non-magnetic cloud collapse (e.g., Bonnell, 1994; Whitworth et al., 1995; Burkert et al., 1997; Matsumoto & Hanawa, 2002).

## 5 Discussion

### 5.1 Fragmentation of Magnetic Clouds

In the previous sections, we have shown that core-forming clouds tend to become unstable to the bar mode after the supercritical core formation. The numerical results presented in this paper and in Li & Nakamura (2002) indicate that fragmentation of magnetic clouds happens readily over a wide range of initial cloud parameters. On the surface, this result appears to contradict earlier numerical studies of magnetic cloud fragmentation, which suggest that strong magnetic fields tend to stifle fragmentation of dynamically collapsing clouds (Dorfi, 1982; Benz, 1984; Phillips, 1986a, b). For example, Phillips (1986a) performed three-dimensional MHD simulations on initially magnetically supercritical clouds and found that none of his clouds show any evidence for fragmentation even in the presence of a large density perturbation with a fractional amplitude of 50%. We emphasize that there are two important differences between our model and the previous ones. First, previous studies considered a frozen-in field, whereas in our study ambipolar diffusion is a crucial ingredient. Our initial cloud, which is supported by a strong magnetic field, is stable to gravitational fragmentation as well as dynamic collapse. The weakening of the magnetic support driven by ambipolar diffusion is what initiates the dynamic collapse and fragmentation. Second, the contribution of the magnetic tension to the cloud support is much weaker in the earlier studies; indeed, most of them took a uniform magnetic field as the initial state (which of course has no magnetic tension). The beneficial effect of magnetic tension on cloud fragmentation was discussed by Boss (2000, 2002), who showed that a strong magnetic tension can slow down cloud contraction appreciably, which tends to enhance cloud fragmentation.

Recently, Tohline (2002) reviewed the theoretical studies on binary star formation, concentrating on the non-magnetic case. He found broad agreement among numerical simulations by various groups that non-magnetic clouds are susceptible to prompt fragmentation if they collapse from an initial configuration that is relatively uniform in density and contains more than a few Jeans masses of material. If supported only by thermal pressure, such a configuration is highly unstable. The question is then (Tohline, 2002): how is nature able to produce such an unstable configuration in the first place? We believe that the key to resolve this apparent difficulty is the presence of a strong magnetic field, which stabilizes the entire cloud even when the cloud mass is well above one thermal Jeans mass 333The cloud can also be supported by turbulence, although the details of cloud evolution driven by turbulence decay remain to be worked out (Myers, 1999). (e.g., Lizano & Shu, 1989). Specific models of core formation in magnetically subcritical clouds appear able to explain the mass distribution and extended subsonic infall motions inferred in several starless cores (e.g., Lee, Myers & Taffala, 2001; Ciolek & Basu, 2000; Indebetouw & Zweibel, 2000) and are consistent with pattern of chemical differentiation observed in L1544 (Li et al., 2002), provided that the initial flux-to-mass ratio lies within a factor of two or so of the critical value, as adopted in this paper.

### 5.2 Elongation of Dense Cloud Cores

Dense cores of molecular clouds are observed to have significant elongation, with a typical aspect ratio of 2 in the plane of sky (Myers et al., 1991). Their true three-dimensional shapes are probably triaxial, as inferred statistically by Jones, Basu & Dubinski (2001) and Goodwin, Ward-Thompson, & Whitworth (2002). We have previously noted that the triaxial core shape can naturally result from nonlinear growth of the bar mode in the direction perpendicular to the field lines, coupled with the flattening of the mass distribution along the field lines (Nakamura & Li, 2002). The models presented in §3 strengthen this conclusion. We find that the aspect ratio of the bar is more or less independent of the initial model parameters. It stays at about 1.5 $``$ 2 during the period when $`\mathrm{\Gamma }`$ decreases from $`1`$ to $`0.5`$, after the core has become supercritical (which makes bar growth possible) but before the ratio of effective gravity to effective pressure force exceeds $`2`$ (when the bar elongation is greatly amplified through the Lin-Mestel-Shu instability). This initial period of supercritical phase is also the phase most likely observed as a starless core, when the core is dense enough to be identified as a core using high density tracers (such as NH3 or CS), but before rapid collapse sets in which quickly produces a stellar object at the center, ending the starless phase.

In this supercritical collapse phase, the decrease of the central flux-to-mass ratio saturates at about half the critical value, almost independent of the initial parameters (see panel \[b\] of Fig. 2). This saturation of flux-to-mass ratio is also observed in axisymmetric models (Basu & Mouschovias, 1995). It means that the supercritical cores remain significantly magnetized, and magnetic fields may not be completely ignored even during the dynamic collapse phase of star formation.

### 5.3 Three Types of Binary and Multiple Star Formation

Observations of nearby star forming regions have established that binary fraction among pre-main sequence stars is at least as high as among the main-sequence stars. In some nearby low-mass star forming regions like Taurus-Auriga, the binary frequency is even higher than that of main-sequence stars. In such low-mass star forming regions, stellar densities ($`n_{}10`$ pc-3) are much too low to allow the formation of binaries by capture or their destruction by gravitational encounters within their ages of $``$ 1 Myr. This indicates that the observed high binary frequencies more or less reflect the pristine population of proto-binaries at the epoch of formation, suggesting that binary and small multiple stars are formed prior to pre-main sequence phase by gravitational fragmentation of parent dense molecular cloud cores. In fact, recent high resolution millimeter observations of dust continuum have revealed evidence of gravitational fragmentation of protostellar cores in the earliest, Class 0 and I phases. Such observations were summarized in Looney et al. (2000), who classified young multiple star systems into three types: (1) separate envelope (separation $`6000`$ AU), (2) common envelope (separation 150 $``$ 3000 AU), and (3) common disk systems (separation $`100`$ AU). Separate envelope systems exhibit clearly distinct centers of gravitationally bound cores with separations greater than 6000 AU and the components are surrounded by low density material (e.g., NGC1333 IRS2, SVS13). Common envelope systems have one primary core of gravitational concentration which apparently breaks up into multiple fragments with separations of 100 $``$ 3000 AU (e.g., IRAS 16293-2422). Common disk systems have circumbinary disk-like envelopes surrounding young stellar objects with separations smaller than 100 AU (e.g., L1551 IRS5).

We have also uncovered three distinct types of fragmentation of magnetic clouds in our numerical experiments. They correspond roughly with the above empirical classification of young binary and multiple star systems. Based upon the results presented in §4, we classify theoretically the process of fragmentation into three types, depending on the initial model parameters: (1) separate core formation, in which the cloud breaks up into two separate cores during the isothermal collapse phase, (2) bar fragmentation, in which the aspect ratio of the bar becomes large enough to break up into multiple fragments along the major axis and the fragmentation takes place after the bar has become opaque to dust emission, and (3) disk fragmentation, in which the central region, after becoming opaque to dust emission, evolves into a rapidly rotating disk, which fragments into multiple blobs. These cases are discussed in turn.

When the initial cloud has a relatively flat central part and is rotating relatively rapidly (e.g., $`\omega 0.2`$ for $`n4`$), it fragments into two separate cores during the isothermal collapse phase in the presence of an $`m=2`$ perturbation after the central region becomes magnetically supercritical (separate core formation). Morphologically, this case resembles the separate envelope systems discussed in Looney et al. (2000). The early phase of fragmentation also looks like the double-core in the extinction map of Fest 1-457 (Alves et al., 2002) (see Laundardt  for another example of double-core, in Bok globule CB 230). When first produced, the dense cores have insufficient orbital rotation to prevent them from spiraling towards each other. We expect their separation to shrink appreciably with time, producing perhaps a highly eccentric system. Whether they will eventually merge together is unclear at the present. We suspect that the cores will survive to become individual star/stellar system, since the density is much higher inside the cores than outside by the end of our calculations. This means that it takes the cores much less time to collapse onto themselves to form stars than to merge together. Moreover, there exists a strong magnetic field outside the cores which can prevent direct core collision. In our separate core formation case, the density of the central inter-core region decreases with time as matter falls towards the cores, whereas the field strength increases due to magnetic flux diffusion out of the cores into the inter-core region. Indeed, the central inter-core region may become magnetically subcritical again at the later stages of evolution. This region serves as a “magnetic cushion” to direct core collision. We propose that the separate core formation is a mechanism for producing wide binaries with separations of $`10^4`$ AU. Each core will have a second chance at fragmentation at higher densities, when a first bar or rotating disk is formed. The secondary fragmentation within each core could lead to the formation of a small (hierarchical) multiple star system. As mentioned earlier, separate core formation in the presence of an $`m=2`$ perturbation has much in common with the multiple fragmentation induced by higher order perturbations that we studied previously in Li & Nakamura (2002). Multiple core formation is an related channel for wide binary and multiple star formation.

When the initial cloud contains a mass well above the thermal Jeans mass and has a moderately centrally condensed density distribution (e.g, $`r_0>15\pi `$ for $`n2`$ and $`\omega 0.10.2`$), a long bar is formed, which can break up into several (pairs of) fragments after the bar becomes opaque to the dust emission (bar fragmentation). The multiple fragments grow in mass due to both merging and accretion. This bar fragmentation case would correspond to the common envelope systems of Looney et al. (2000), since the whole bar is embedded within a common core out of which the bar is formed. The bar fragmentation occurs at a density of $`10^{12}`$ cm-3. According to linear perturbation theory, a highly elongated bar breaks up gravitationally into fragments, with separations comparable to the wavelength of the most unstable linear perturbation, which corresponds to $`10^2`$ AU at $`10^{12}`$ cm-3. The typical masses of the fragments are in the substellar regime. Some of the fragments could be ejected out of the system in the course of dynamical interaction (see also Bate et al., 2002). Those ejected may become isolated brown dwarfs (Reipurth & Clarke, 2001) and perhaps even free-floating massive planets. Interestingly, the cloud parameters derived by Ciolek & Basu (2000) for the well-studied starless core L1544 put the cloud in the bar fragmentation category. Therefore, L1544 appears to be a progenitor of a binary or small multiple star system of the common envelope type.

When the initial cloud mass in the central plateau region is comparable to the thermal Jeans mass (e.g., $`r_0<15\pi `$ for $`n2`$ and $`\omega 0.1`$), the bar mode does not grow much during the isothermal collapse, and it is difficult for the (relatively short) bar to break up into fragments (disk fragmentation). In this case, the central region evolves into a rapidly rotating disk whose size depends on the rate of initial rotation. The interaction of the central disk with the bar generates complex structures in the disk where a secondary bar mode is excited by rotational instability if $`\beta >0.3`$. Typically, the size of centrifugal radius is of the order of 10 AU. When the disk fragments into blobs, the separations of the blobs should be of the same order of the disk radius. This case may be responsible for the formation of binary stars with a separation smaller than $`10^2`$ AU. This case would be related to the common disk system of Looney et al. (2000).

As mentioned in § 2, we did not consider angular momentum transfer in the vertical direction, which can brake the cloud rotation efficiently during the relatively long magnetically subcritical phase of evolution (Nakano & Tademaru, 1972; Basu & Mouschovias, 1994). Magnetic braking is thought to be responsible for the low rates of rotation observed in dense cores (Goodman et al., 1993) and adopted in this paper. Once the central region of the cloud becomes supercritical, the peak density increases in a runaway fashion, leaving little time for the magnetic braking to operate efficiently (Basu & Mouschovias, 1994). We therefore expect the formation of separate cores, which occurs in the runaway collapse phase, to be little affected by magnetic braking. After a bar or a disk becomes dense enough to be opaque to dust emission, it can be supported by thermal pressure for a time long compared to the local dynamic time, and magnetic braking can in principle become efficient again (Krasnopolsky & Konigl, 2002). However, at such high densities, magnetic fields start to decouple from matter (Nishi, Nakano, & Umebayashi, 1991), which could render the braking inefficient. We will postpone a detailed treatment of the magnetic braking including the effects of decoupling to a future work.

## 6 Conclusion

In this paper, we have followed the early phases of cloud fragmentation in magnetic clouds, using a two-dimensional MHD code based on the thin-disk approximation. At the end of our calculations, we obtained dense, discrete fragments out of initially axisymmetric clouds perturbed by an $`m=2`$ mode of modest fractional amplitude $`A=0.05`$. The masses of the fragments are still small, on the order of $`10^210^1`$ M. We expect the fragments to grow to stellar masses, by merging and further accretion, except perhaps those ejected dynamically out of the system soon after the fragmentation. The masses of those ejected may remain frozen in the substellar regime.

We found three types of fragmentation in magnetic clouds. They are (1) a cloud breaking up during the isothermal phase of cloud evolution into separate cores, each of which collapses more or less independently into a star/stellar system (panels \[b\] and \[c\] of Fig. 7), (2) a core collapsing into a needle-like, opaque “first bar”, which breaks up into several fragments (Fig. 8), and (3) a rapidly rotating, opaque “first disk” producing self-gravitating blobs, which may survive to become seeds of relatively close binary and multiple systems (Figs. 9 and 10). The fragments in the above cases share, respectively, a common cloud, common core, and common disk. They correspond loosely to the empirical classification of independent envelope, common envelope and common disk systems for the youngest, embedded binary and multiple systems. Much work is needed to firm up this connection.

In the near future, we plan to extend our numerical simulations beyond the initial fragmentation to the protostellar accretion phase, to determine the physical properties of the binary and multiple systems formed in magnetic clouds, such as their eccentricity, mass ratio, and orbital period. In a longer run, we hope to generalize the calculations to 3D, which will allow for a detailed study of magnetic braking, which may have a profound impact on the final orbits of the formed stellar system.

Numerical computations in this work were carried out at the Yukawa Institute Computer Facilities, Kyoto University. F.N. gratefully acknowledges the support of the JSPS Postdoctoral Fellowships for Research Abroad, and Z.Y.L. is supported in part by NASA.


format: text

text:
# The complete submillimetre spectrum of NGC 891

## 1 Introduction

NGC 891 is one of the most perfect edge-on \[Sofue & Nakai 1993\] spiral galaxies in our neighbourhood. This peculiarity, and its similarity with the Milky Way galaxy in Hubble type (Sb), optical luminosity and rotational velocity \[van der Kruit 1984\], has allowed it to be intensively studied in order to understand the distribution and the characteristics of the galactic interstellar medium. It has been widely studied in radio continuum (e.g. Allen et al. 1978, Dahlem et al. 1994), carbon monoxide (García-Burillo et al. 1992, Scoville et al. 1993, Sofue & Nakai 1993, García-Burillo & Guélin 1995, Sakamoto et al. 1997), atomic hydrogen (for instance Swaters et al. 1997), molecular hydrogen \[Valentijn & van der Werf 1999\] and optical extinction \[Howk & Savage 1997\]. The dust far-infrared continuum emission has been observed by the IRAS satellite (Wainscoat et al. 1997, Rice et al. 1988), the 1.3 mm emission by Guélin et al. (1993), and the submillimetre emission has recently been mapped with SCUBA at 450 and 850 $`\mu `$m by Alton et al. (1998) and Israel et al. (1999). For instance, Alton et al. (1998) showed that the dust submillimetre emission could be fitted by a two-dust component model, assuming a spectral index equal to 2. However, this was a fit with zero degree of freedom, thus it needs to be confirmed by including other submillimetre measurements. The accurate knowledge of dust properties in this galaxy, and therefore in other spiral galaxies such as the Milky Way, needs a better submillimetre spectral coverage, in order to properly derive the dust temperature(s) and spectral index. This is crucial to derive accurate masses of the interstellar medium in galaxies.

In this article, we present new results obtained in the submillimetre range with the ISOPHOT \[Lemke et al. 1996\] instrument on board of ESA’s Infrared Space Observatory \[Kessler et al. 1996\], and with the PRONAOS (PROgramme NAtional d’Observations Submillimétriques) balloon-borne telescope. Section 2 presents the observations and data processing, Section 3 presents the results obtained, and Section 4 presents an analysis of the data set.

## 2 Observations and data processing

### 2.1 PRONAOS data

PRONAOS is a French balloon-borne submillimetre experiment. Four bolometers cooled at 0.3 K measure the submillimetre flux with sensitivity to low brightness gradients of about 4 MJy/sr in band 1 (200 $`\mu `$m) and 0.8 MJy/sr in band 4 (580 $`\mu `$m). The effective wavelengths are 200, 260, 360 and 580 $`\mu `$m, and the angular resolutions are 2 in bands 1 and 2, 2.5 in band 3 and 3.5 in band 4. Details about the instrument can be found in Ristorcelli et al. (1998) and Lamarre et al. (1994). The data which we analyze here were obtained during the second flight of PRONAOS in September 1996, at Fort Sumner, New Mexico. The data processing and the map-making method, including deconvolution from chopped data, are described in Dupac et al. (2001). However, for these NGC 891 data, we use a non prior deconvolution method, implemented the same way as the Wiener filter in Dupac et al. (2001). Indeed, these data are quite noisy, so the optimal Wiener filter may somehow modify the signal, which is not the case with the non prior map-making method. Then we carefully subtract the background around the galaxy to get the corrected surface brightness, and we smooth the maps of the first three bands with adequate profiles in order to obtain the same angular resolution as the fourth band (3.5). After having done this, we check the correlation of the pixels in each band with respect to all the other bands in order to detect possible residual offsets between bands. This happened in band 2 that we rescaled to match the other bands. The noise fluctuations in the final maps define the spatial relative uncertainty which is 5 MJy/sr at 200 $`\mu `$m, 4 MJy/sr at 260 $`\mu `$m and 2 MJy/sr at 360 and 580 $`\mu `$m. The intercalibration error between bands is 5 % (1 $`\sigma `$), and the absolute calibration uncertainty is 8 %.

### 2.2 ISOPHOT data

NGC 891 has been mapped with ISOPHOT using the Astronomical Observation Template PHT32 in chopped mapping mode with the C200 array detector (2$`\times `$2 pixels, 92′′ per pixel) at the reference wavelengths 170 and 200 $`\mu `$m (see details in Laureijs et al. 2001). The data have been processed using the ISOPHOT Interactive Analysis software PIA v9.1 \[Gabriel et al. 1997\]. The data reduction procedure includes ramp linearisation, ramp deglitching, reset interval correction, dark current subtraction, signal linearisation and signal deglitching. Each observation has been bracketted with two measurements of the fine calibration source (FCS). For the flux calibration, we have used the second FCS. The first quartile normalisation flat-fielding method has been used as implemented in PIA in order to correct for the remaining responsivity differences of the individual detector pixels. We smooth the derived 1.5-resolution maps to 3.5 in order to obtain a consistent data set. The absolute photometric error of ISOPHOT is around 20 %.

### 2.3 IRAS and SCUBA data

We use the SCUBA 450 and 850 $`\mu `$m flux profiles, as well as the high resolution (HiRes) 60 and 100 $`\mu `$m flux profiles from the InfraRed Astronomical Satellite survey, as published in Alton et al. (1998). We smooth the original profiles to the angular resolution of the fourth band of PRONAOS (3.5).

## 3 Results

The information obtained is presented as intensity profiles along the major axis in Fig. 1. We have checked the photometry of all 3.5-resolution maps with respect to all the others, by plotting pixel-pixel diagrams. No noise offset between bands could be detected. In particular, both 200 $`\mu `$m bands (ISOPHOT and PRONAOS) are in good agreement with each other, which gives a good global consistency to the data set.

All bands of PRONAOS and ISOPHOT show a good symmetry between both sides of the galaxy, which extends approximatively 4 away from the centre. This extension is similar to the one of the SCUBA maps of Alton et al. (1998) and Israel et al. (1999). Fig. 1 shows a large enhancement of the intensities at an angular distance of about 1.5 away from the galactic centre along the north-eastern side of the galaxy. This is clearly visible in the 200 and 260 $`\mu `$m PRONAOS bands which have the adequate angular resolution (2), as well as in both ISOPHOT bands. This enhancement can also be seen (at the same place) in the SCUBA maps presented by Alton et al. (1998) and Israel et al. (1999). Another enhancement, though less large, is visible at the same angular distance on the other side of the galaxy on SCUBA maps. The explanation for these enhancements could be spiral arms or a molecular ring encircling the galactic bulge.

## 4 Analysis

### 4.1 The dust submillimetre emission

We present in Fig. 2 the spectra obtained towards the galactic centre and 3 away from the centre on either side. The plotted error bars are dominated by the global calibration uncertainty of each instrument. As can be seen in Fig. 2, the submillimetre spectral energy density is very well constrained with the present data.

Several fitting procedures have been applied to the three spectra displayed in Fig. 2. The usual simple modelization of the large-grain emission is the modified blackbody, which obeys the following equation:

$$I_\nu =ϵ_0B_\nu (\lambda ,T)(\frac{\lambda }{\lambda _0})^\beta $$
(1)
where $`\lambda `$ is the wavelength, $`ϵ_0`$ the emissivity of the observed dust column density at $`\lambda _0`$, T the temperature of the grains, $`\beta `$ the spectral index and $`B_\nu `$ the Planck function. We fit the data using a one-component modified blackbody model, with and without the 60 $`\mu `$m point (which may be contaminated by very small grain emission, see Désert et al. 1990) and assuming $`\beta `$=2 or not. We also fit the data using a two-component model, assuming $`\beta `$=2 for both components, which have different temperatures. A summary of the results is shown in Table 1.

The derived parameters (temperatures and spectral indices) are consistent for the three fitted positions (centre, north-east and south-west). However, the central peak exhibits a slightly higher temperature than the wings. This is also true for the two-component fits and the fits with the 60 $`\mu `$m data. The one-component fits with $`\beta `$ free and without the 60 $`\mu `$m data show that the wings may have a spectral index slightly lower than 2, while the centre exhibits a low spectral index (1.4). The results show an anticorrelation between the temperature and the spectral index, which is consistent with previous results obtained by Dupac et al. (2001), Dupac et al. (2002) and Dupac et al. (2003). However, if one considers $`\beta `$=2 to be the rule, then the central peak still exhibits a higher temperature (19 K) than the wings (18 K). This is less than what derived by Israel et al. (1999) under the same assumption of $`\beta `$=2, with IRAS 100 $`\mu `$m and SCUBA data only (21 K). Fitting two dust components to the emission longwards of 100 $`\mu `$m does not give very significant results: the warm temperature is not much changed with respect to the one-component fit ($`\beta `$=2) temperature, and the cold component is indeed extremely cold (5-9 K). Therefore, and given the goodnesses of fit presented in Table 1, we consider that the dust emission longwards of 100 $`\mu `$m is very adequatly fitted by a single dust population.

The fits including the 60 $`\mu `$m emission show that the one-component modeling is still possible, especially toward the centre where the $`\chi ^2`$/d.o.f. value is low for the fit with $`\beta `$ free. In this case, the derived spectral indices are low: 1.15 toward the central peak and around 1.4 toward the wings. However, the $`\chi ^2`$/d.o.f. values are clearly larger than those of the fits longwards of 100 $`\mu `$m. This may indicate that a two-component model is useful for fitting with the 60 $`\mu `$m emission. This fit gives temperatures of 29 and 16 K, as well for the wings as for the centre. Given this quite high temperature of 29 K, it is likely that the 60 $`\mu `$m emission is dominated by big grains, which gives an a posteriori consistency to this two-component model.

### 4.2 Total fluxes and the mass of the interstellar medium

For calculating the total fluxes (presented in Table 2), we integrate the measurements over the three 3 beams along the major axis, which corresponds to 25 kpc. For estimating the mass of the galactic interstellar medium, we follow the simple model described in Dupac et al. (2001), which uses the dust 100 $`\mu `$m opacity in the diffuse interstellar medium from Désert et al. (1990). In this model, the total column density of interstellar medium is simply proportional to $`ϵ_{100}`$, the coefficient being 1.67$`\times `$1024 H.cm-2. Two mass estimations are made: one assumes a single dust component with the emissivity $`ϵ_{100}`$ varying along the major axis as derived by the free-$`\beta `$ fits without the 60 $`\mu `$m data (see Table 1); the other assumes two dust components with varying emissivities as derived by the two-component fits with the 60 $`\mu `$m data. In either case, we assume a distance to NGC 891 of 9.5 Mpc \[van der Kruit & Searle 1981\]. We find a mass of the interstellar medium of 4.6$`\times `$109 M when applying the one-component free-$`\beta `$ model, whereas we find 12$`\times `$109 M when applying the two-component model (including the 60 $`\mu `$m points). This discrepancy is due to the existence of quite large amounts of cold (15.7 K) dust in the two-component results. The value of 4.6$`\times `$109 M is consistent with previous estimates from Guélin et al. (1993): 4$`\times `$109 M, and Israel et al. (1999): 3.9$`\times `$109 M. We can assume from Guélin et al. (1993) that the atomic hydrogen mass is 2.5$`\times `$109 M, which gives a molecular hydrogen amount of 2.1$`\times `$109 M from our measurement (one-component fit). This is in good agreement with the standard assumption of similar quantities of atomic and molecular gas in spiral galaxies. However, if we trust the two-component estimation, the molecular to atomic gas ratio becomes 3.8, which is much larger than standard assumptions. This amount of molecular gas (9.5$`\times `$109 M) is not in good agreement with the standard estimation from CO measurements either, using the CO-to-H2 conversion factor of Strong et al. (1988): about 4.5$`\times `$109 M. From this two-component modeling, we derive cold over warm component mass ratios for the three studied positions: 19 for the centre, 42 for the north-eastern data point and 34 for the south-western one. Therefore, there is a strong trend to consider that the cold (16 K) dust amount is much larger than the warm dust amount, and that the cold to warm dust ratio is twice larger in the wings of this galaxy than in the centre. These ratios are smaller than those derived by Alton et al. (1998) with IRAS and SCUBA data only, but larger than those derived by Israel et al. (1999) with slightly different temperatures (cold to warm dust ratio of about 10).

As our ISOPHOT and PRONAOS measurements constrain very well the peak of the spectra, we are confident in the reliability of the derived dust properties (emissivity, temperature and spectral index). However, the different possibilities to fit (well) the spectra lead to different estimates of the mass. If one accepts the possibility that the spectral index of the dust can be somewhat different from 2, which seems wise to us regarding recent results (Walker et al. 1990, Ristorcelli et al. 1998, Dupac et al. 2001, Dupac et al. 2002, Dupac et al. 2003, Bennett et al. 2003), it does not seem possible to better discriminate between the different fitting procedures given the present data, although the spectral coverage is relatively accurate. Also, the $`ϵ`$/NH ratio used \[Désert et al. 1990\] is debatable, and the mass estimates are of course very much dependent on this value.

We compare these ISM mass estimates to the dust mass (19$`\times `$106 M) of NGC 891 obtained by Alton et al. (2000), assuming given values for the grain radius and density. This is substantially less than the estimation of 50$`\times `$106 Mfrom Alton et al. (1998), who used a two-component model of the submillimetre emission. Comparing the Alton et al. (2000) dust mass to our ISM mass estimates, we obtain a gas to dust mass ratio of 240 for the one-component mass estimate and 640 for the two-component one. The first estimate of this ratio is slightly less than that derived by Alton et al. (2000): 260, a bit closer to the Milky Way value and well in agreement with the canonical dependence on the metallicity \[Issa et al. 1990\]. If we compare the ISM mass estimate from our two-component model (12$`\times `$109 M) to the two-component model estimate from Alton et al. (1998): 50$`\times `$106 M, then the gas to dust ratio is (also) 240. We therefore consider that this value is somewhat robust.

## 5 Conclusion

We have presented PRONAOS and ISOPHOT data of the NGC 891 spiral galaxy, together with IRAS and SCUBA measurements. The overall submillimetre spectrum of this object is now well constrained, and two different conclusions can be derived: either one single dust component is fitted to the spectra (without the 60 $`\mu `$m data), and in this case the temperature is higher toward the centre of the galaxy than toward the wings, and the spectral index is less steep (1.4 in the centre and 1.7-2 in the wings); either two components are fitted assuming $`\beta `$=2, and this leads to a warm temperature of 29 K and a cold temperature of 16 K. The interstellar medium masses derived by these two methods are quite different, but both are consistent with a gas to dust mass ratio of 240.

The good spectral coverage of the submillimetre range that we present here makes us rather confident in our determination of the dust emission properties. However, the way to fit the data, as well as the opacity properties, lead to inevitable uncertainties in the determination of the masses. A next step in the accurate knowledge of the continuum emission of this galaxy, as well as of other nearby spirals, would be to obtain the same spectral resolution in the same wavelength range with a better angular resolution. This should allow to better discriminate the number of distinct components, as well as the possible variations of the spectral index, and this should be achieved by forthcoming experiments such as the Herschel satellite with the instruments PACS \[Poglitsch et al. 2001\] and SPIRE \[Griffin et al. 2001\].

## 6 Acknowledgments

We thank very much R. Tuffs for his help. We are indebted to the French space agency Centre National d’Études Spatiales (CNES), which supported the PRONAOS project. We are very grateful to the PRONAOS technical teams at CNRS and CNES, and to the NASA-NSBF balloon-launching facilities group of Fort Sumner (New Mexico). We thank very much the technical and scientific teams of ISO, an ESA project with instruments funded by ESA member states (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) with participation of ISAS and NASA.