# MAGNETICALLY SELF-REGULATED FORMATION OF EARLY PROTOPLANETARY DISKS

Patrick Hennebelle, Benoît Commerçon, Gilles Chabrier, Pierre Marchand

► To cite this version:

Patrick Hennebelle, Benoît Commerçon, Gilles Chabrier, Pierre Marchand. MAGNETICALLY SELF-REGULATED FORMATION OF EARLY PROTOPLANETARY DISKS. The Astrophysical Journal Letters, 2016, 830 (1), pp.L8. 10.3847/2041-8205/830/1/L8. http://hal-02317626

HAL Id: hal-02317626

https://hal.science/hal-02317626v1

Submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MAGNETICALLY SELF-REGULATED FORMATION OF EARLY PROTOPLANETARY DISKS

PATRICK HENNEBELLE¹,², BENOÎT COMMERÇON³, GILLES CHABRIER³,⁴, AND PIERRE MARCHAND³

¹Laboratoire AIM, Paris-Saclay, CEA/IRFU/SAP—CNRS—Université Paris Diderot, F-91191, Gif-sur-Yvette Cedex, France

²Laboratoire de radioastronomie, UMR CNRS 8112, École normale supérieure et Observatoire de Paris, 24 rue Lhomond, F-75231 Paris cedex 07, France

³École normale supérieure de Lyon, CRAL, UMR CNRS 5574, Université de Lyon, F-69364 Lyon Cedex 07, France

⁴School of Physics, University of Exeter, Exeter EX4 4QL, UK

Received 2016 July 8; revised 2016 August 3; accepted 2016 August 8; published 2016 October 6

## ABSTRACT

The formation of protoplanetary disks during the collapse of molecular dense cores is significantly influenced by angular momentum transport, notably by the magnetic torque. In turn, the evolution of the magnetic field is determined by dynamical processes and non-ideal MHD effects such as ambipolar diffusion. Considering simple relations between various timescales characteristic of the magnetized collapse, we derive an expression for the early disk radius, $r \simeq 18$ au ($\eta_{AD}/0.1$ s)$^{2/9}$($B_z/0.1$ G)$^{-4/9}$($M/0.1 M_\odot)^{1/3}$, where M is the total disk plus protostar mass, $\eta_{AD}$ is the ambipolar diffusion coefficient, and $B_z$ is the magnetic field in the inner part of the core. This is significantly smaller than the disks that would form if angular momentum was conserved. The analytical predictions are confronted against a large sample of 3D, non-ideal MHD collapse calculations covering variations of a factor 100 in core mass, a factor 10 in the level of turbulence, a factor 5 in rotation, and magnetic mass-to-flux over critical mass-to-flux ratios 2 and 5. The disk radius estimates are found to agree with the numerical simulations within less than a factor 2. A striking prediction of our analysis is the weak dependence of circumstellar disk radii upon the various relevant quantities, suggesting weak variations among class-0 disk sizes. In some cases, we note the onset of large spiral arms beyond this radius.

**Key words:** diffusion – gravitation – hydrodynamics – magnetohydrodynamics (MHD) – protoplanetary disks – turbulence

## 1. INTRODUCTION

Circumstellar disks are of fundamental importance in astrophysics because they are the birth sites of planet formation. Yet, our current understanding of centrifugally supported disks lack a clear description of how and when they form. The exact role played by magnetic field, in particular, remains an unsettled issue. Various teams have been consistently finding that catastrophic braking may occur when the magnetic field and the rotation axis are aligned (Allen et al. 2003; Galli et al. 2006; Price & Bate 2007; Hennebelle & Fromang 2008; Mellon & Li 2008). In such circumstances, magnetic braking can be so intense that the formation of primordial disks at the class-0 stage can be suppressed even for modest magnetizations. Although recent observations have revealed that large disks are rare, if not absent, at class-0 stage (Maury et al. 2010; Tobin et al. 2015), complete inhibition of disk formation cannot be considered as a plausible scenario because of (i) the assumed aligned configuration (Hennebelle & Ciardi 2009; Joos et al. 2012, 2013; Santos-Lima et al. 2012; Seifried et al. 2012) and (ii) the ideal MHD assumption (Dapp & Basu 2010; Li et al. 2011; Krasnopolsky et al. 2012; Tomida et al. 2010, 2015; Tsukamoto et al. 2015; Masson et al. 2016; Wurster et al. 2016). Indeed, when either misalignment, turbulence, or non-ideal MHD effects are included, disks tend to form more easily, although by no means as large as in pure hydrodynamical calculations (Tomida et al. 2015; Masson et al. 2016).

Naively, a broad distribution of disk properties might be expected, depending for instance on the core mass, the amount of rotation or turbulence in the core, and the strength of the field and its configuration. In this paper, we derive a theoretical framework that suggests the opposite, i.e., that disks at their early stages are remarkably regulated by a combination of

magnetic braking and non-ideal MHD effects, leading to similar sizes. In Sections 2 and 3, we develop simple analytical arguments leading to our suggestion. In Section 4, we compare these analytical estimates of the disk sizes to a series of collapse calculations corresponding to a large variety of initial conditions. Section 5 concludes the paper.

## 2. MAGNETIC FLUX DISTRIBUTION

We first aim at assessing the intensity of magnetic braking. Recent 3D simulations (Tomida et al. 2015; Masson et al. 2016) found that ambipolar diffusion leads to a plateau, i.e., a nearly uniform magnetization, in the inner part of the collapsing cores, with typical values of the order of 0.1 G for $1 M_\odot$ cores, and up to ~0.3 G or so for 100 $M_\odot$ cores. To understand this property we assume stationarity, and reduce the problem to one radial dimension (all quantities in the following are simply written $x \equiv x(r)$). In that case, the Faraday equation reduces to

$$v_r B_z \approx \frac{c^2 \eta_{AD} \partial_r B_z}{4\pi}, \quad (1)$$

where $v_r$ denotes the radial velocity, $B_z$ the field vertical (poloidal) component, and $\eta_{AD}$ the ambipolar diffusivity.

Let us remind that in ideal MHD, flux and mass conservation inside concentric cylinders lead to $B_z \propto \Sigma \simeq 2\rho h$, where $\Sigma$ is the column density and $h$ the typical thickness. Assuming mechanical vertical equilibrium, we get $h \simeq C_s/\sqrt{\pi G\rho}$ and therefore $B_z \propto \sqrt{\rho}$.

Because of the complex dependence of the resistivity $\eta_{AD}$ upon magnetic intensity $B$ and density $n$ (Marchand et al. 2016), Equation (1) cannot be solved analytically. To get a solution of Equation (1), we calculated a table of resistivities for a series of densities and magnetic intensities
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top panel: vertical magnetic field component given by Equation (1), $B_z$, as a function of density $\rho$, for two magnetic intensities (solid lines: 100 µG, dashed lines: 30 µG) at the core boundary and three density and velocity profiles. Middle panel: same as top one except that flux freezing is assumed until $t_{ff} \simeq t_{AD}$. The value of $B_z$ varies only by a factor of a few. Bottom panel: freefall and ambipolar times for the four models displayed in red in the middle and bottom panels. At high densities the freefall time is longer than the ambipolar diffusion time.</figcaption></figure>

from Marchand et al. (2016), from which we get $\eta_{AD}$ for any values by interpolation. To integrate numerically Equation (1), the density and the radial velocity must be specified and we set (e.g., Larson 1969; Shu 1977)

$$ 
\begin{aligned}
\rho(r) &= \delta \frac{C_s^2}{2\pi Gr^2}, \\
v_r(r) &= V_0 (r/r_0)^{-1/2}.
\end{aligned} 
\qquad (2) 
$$

We considered two values of the external field, namely $B_z = 30$ and 100 µG, as well as three densities and velocity field amplitudes, two typical of low-mass cores and one typical of high-mass ones. For the low-mass cores, we took $\delta = 1$ (i.e., the singular isothermal sphere, sis) and $V_0 = 2 \times C_s$ at $r_0 = 10$ au and twice this value (which corresponds to a faster collapse). For the high-mass core, we took $\delta = 10$ and $V_0 = 4 \times C_s$.

Figure 1 displays the results. Since in real cores the flux distribution is due to a combination of flux freezing and

ambipolar diffusion, we explore two cases. First, we solve Equation (1) from the edge to the center of the core (top panel). Second, we assume flux freezing up to the point where the freefall, $t_{ff} = \sqrt{3\pi/(32G\rho)}$, and the ambipolar diffusion times, $t_{AD} = 4\pi/(c^2\eta_{AD}r^2)$, become comparable ($10^{-15}$ g cm$^{-3}$ for 100 µG and $3 \times 10^{-14}$ g cm$^{-3}$ for 30 µG; see bottom panel), then we solve Equation (1) (middle panel). Red lines correspond to $V_0 = 2 \times C_s$ and dark ones to $V_0 = 4 \times C_s$. The dashed lines display the low magnetization cases and solid lines the high magnetization ones. As seen, the value of the magnetic intensity outside the core has a weak influence on the value inside it. A clear transition occurs between a slightly sublinear regime (where $B_z \propto \sim \rho^{2/3}$) to a plateau at $\rho \simeq 10^{-15}$ g cm$^{-3}$. From Figure 5 of Marchand et al. (2016), we see that indeed $\eta_{AD}$ displays two different regimes, which correspond to densities respectively below and above $\simeq 10^{-15}$ g cm$^{-3}$. The bottom panel of Figure 1 shows the freefall time and the ambipolar diffusion time for the four low-mass cases displayed by the red lines in the top and middle panels. Clearly, while the freefall time is shorter than the ambipolar diffusion time in the outer part of the core, the reverse is true in the core inner part (Nakano et al. 2002). The magnetic field at the center of the core weakly depends on the physical conditions and remains remarkably constant. The assumption of stationarity is also well justified as the freefall time is much longer than the ambipolar diffusion one. These features agree quite well with the 3D simulations performed by Masson et al. (2016; their Figure 1) and Tomida et al. (2015). Furthermore, we see that a slowly collapsing low-mass core (red lines) has a smaller central magnetic intensity than a more rapidly collapsing one (black lines). The massive cores, which have both a large inward velocity and a large density, display even higher central magnetic intensities.

Altogether the variations of the magnetic field in the inner part of the core remain limited and weakly depend on the initial conditions.

### 3. THEORETICAL ESTIMATE OF THE DISK RADIUS

To obtain an estimate of the disk radius, we examine the relevant timescales and we estimate the various quantities at the disk centrifugal radius location, i.e., at the disk-envelope boundary. Let us stress that the envelope and the disk are connected by an accretion shock that is quite thin. Therefore, outside the disk, the gas in its vicinity is nearly in freefall (see for example Figures 3 and 4 of Hennebelle & Fromang 2008) and the results of Section 2 can be applied.

#### 3.1. Timescales and Equilibria

The first important timescales are the ones that control the evolution of the azimuthal magnetic field, $B_\phi$, which is responsible for the magnetic braking. On one hand, $B_\phi$ is generated by the differential rotation on a timescale $\tau_{far}$, and on the other hand it is diffused vertically by ambipolar diffusion on a timescale $\tau_{diff}$, with

$$ 
\begin{aligned}
\tau_{far} &\approx \frac{B_\phi h}{B_z v_\phi}, \\
\tau_{diff} &\approx \frac{4\pi h^2 B_z^2 + B_\phi^2}{c^2 \eta_{AD} B_z^2} \approx \frac{4\pi h^2}{c^2 \eta_{AD}},
\end{aligned} 
\qquad (3) 
$$

where h denotes the thickness of the disk.
The second relevant timescales are the magnetic braking one and the rotation time. They are given by

$$ 
\begin{aligned}
\tau_{\text{br}} &\approx \frac{\rho v_{\phi} 4\pi h}{B_z B_{\phi}}, \\
\tau_{\text{rot}} &\approx \frac{2\pi r}{v_{\phi}},
\end{aligned}
\qquad (4) $$

where $r \equiv r_d$ denotes the disk radius.

Then, we assume that the gas in the neighborhood of the disk outer part has a Keplerian velocity (in practice it may be a little lower) and is roughly in vertical hydrostatic equilibrium:

$$ v_{\phi} \approx \sqrt{\frac{G(M_{*} + M_{d})}{r}}, \qquad (5) $$

$$ h \approx \frac{C_s}{\sqrt{4\pi G(\rho + \rho_*)}}, \qquad (6) $$

where $M_d$ is the mass of the disk, $M_*$ the mass of the central star, and $\rho_* = M_*/(4\pi)r_*^{-3}$.

Finally, the density in the envelope is given by

$$ \rho(r) = \delta \frac{C_s^2}{2\pi G r^2} \left( 1 + \frac{1}{2} \left( \frac{v_{\phi}(r)}{C_s} \right)^2 \right). \qquad (7) $$

Apart from $\delta$, which is a coefficient on the order of a few, the first term simply corresponds to the singular isothermal sphere (Shu 1977) while the second one is a correction that must be included when rotation is significant, particularly in the inner part of the envelope close to the disk edge, as discussed in Hennebelle et al. (2004; see their appendix and Figure 2). Note that for massive stars, $\delta$ may be up to about 10 as shown in Figure 3 of Hennebelle et al. (2011).

## 3.2. Dependence of the Disk Radius

The disk properties are the result of the balance between various quantities at the disk-envelope boundary. First of all, as mentioned above, the generation of the toroidal field through differential rotation is offset by the ambipolar diffusion in the vertical direction. From Equation (3), with $\tau_{\text{far}} \approx \tau_{\text{diff}}$, we get

$$ \frac{B_{\phi}}{h v_{\phi}} \approx \frac{4\pi}{c^2 \eta_{\text{AD}}} B_z. \qquad (8) $$

Second of all, the braking and the rotation timescales must be of same order, $\tau_{\text{br}} \approx \tau_{\text{rot}}$, yielding

$$ B_{\phi} \approx \frac{2h\rho}{r} v_{\phi}^2 B_z^{-1}. \qquad (9) $$

Combining Equations (8) and (9) yields

$$ \frac{2\rho}{r} v_{\phi} \approx \frac{4\pi}{c^2 \eta_{\text{AD}}} B_z^2, \qquad (10) $$

while vertical and radial equilibria at the disk outermost limit imply

$$ \frac{\delta G^{1/2} (M_d + M_*)^{3/2}}{2\pi r^{9/2}} \approx \frac{4\pi}{c^2 \eta_{\text{AD}}} B_z^2, \qquad (11) $$

where, for the sake of simplicity, we have assumed $\rho \propto v_{\phi}^2$ in Equation (7).

All these relations lead to

$$ r_{\text{d,AD}} \approx \left( \frac{\delta G^{1/2} c^2 \eta_{\text{AD}}}{8\pi^2} \right)^{2/9} B_z^{-4/9} (M_d + M_*)^{1/3}. \qquad (12) $$

The mass of the star/disk system, $M_d + M_*$, grows as the envelope gets accreted. We take $0.1 M_\odot$ as a fiducial value since we are investigating the class-0 phase.

With these values, Equation (12) can be rewritten:

$$ r_{\text{d,AD}} \approx 18 \text{ au} \\ \times \delta^{2/9} \left( \frac{\eta_{\text{AD}}}{0.1 \text{ s}} \right)^{2/9} \left( \frac{B_z}{0.1 \text{ G}} \right)^{-4/9} \left( \frac{M_d + M_*}{0.1 M_\odot} \right)^{1/3}. \qquad (13) $$

The striking result illustrated by Equation (13) is the weak dependence of the disk radius upon all involved quantities. Note that, in principle, the magnetic resistivity $\eta_{\text{AD}}$ depends on density (see Figure 5 of Marchand et al. 2016), but this dependence is very shallow. We find a more pronounced, although still moderate, dependence of the radius upon $B$ as $\sim B_z^{-0.5}$. In principle this could introduce some variations among the disk radii but, as seen in Section 2, the magnetic field in the inner part of the envelope of the cores is also regulated by ambipolar diffusion. Finally, the radius depends also weakly on the mass. Indeed, as accretion proceeds, the disk is expected to become only about twice larger when the star becomes 10 times more massive, i.e., $M_* = 1 M_\odot$.

We also note that $C_s$ does not enter explicitly in Equation (13), suggesting weak dependence of the disk radius on the velocity field, be it purely thermal or turbulent (through an effective sound speed $C_s^{\text{eff}} = (C_s + \langle v_{\text{rms}}^2 \rangle^{1/2})^{1/2}$), as indeed is found in the simulations (see below). In practice, some dependence on the various supports enters in the coefficient $\delta$ but since it appears at the power 2/9 this leads to weak variations.

It is interesting to compare these trends with the purely hydrodynamical case. Let us consider a spherical cloud of density $\rho_0$ in solid body rotation at a rate $\Omega_0$. When a fluid particle initially at radius $R_0$ reaches centrifugal equilibrium in the disk, its radius becomes

$$ r_{\text{d,hydro}} \approx \frac{\Omega_0^2 R_0^4}{4\pi/3 \rho_0 R_0^3 G} = 3\beta R_0 \\ = 106 \text{ au} \frac{\beta}{0.02} \left( \frac{M}{0.1 M_\odot} \right)^{1/3} \left( \frac{\rho_0}{10^{-18} \text{ g cm}^{-3}} \right)^{-1/3}, \qquad (14) $$

where $\beta = R_0^3 \Omega_0^2 / 3GM$ denotes the core rotational support. Whereas the mass dependence remains the same as above, the radius now strongly (quadratically) depends on the initial rotation rate. As cores have typical $\beta \approx 0.02$ (Goodman et al. 1993; Belloche 2013), purely hydrodynamical disks should be on average significantly (typically 5–6 times) larger than the ones we predict.

# 4. COMPARISON WITH NUMERICAL SIMULATIONS

## 4.1. Initial Conditions

To test the validity of our analytical model, we have performed two series of numerical simulations of non-ideal MHD collapse with ambipolar diffusion with the RAMSES code (Teyssier 2002; Fromang et al. 2006; Masson et al. 2012).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Ratio of the disk radius measured in the simulations over the theoretical estimate (Equation (13)) as a function of the total (disk plus star) mass. The left panel corresponds to the low-mass cores and the right panel to the high-mass ones. As seen, the ratio is on average of the order of 1. The points that deviate significantly ($r_{d,simu}/r_{d,AD} \simeq 2-3$) correspond to the development of prominent spiral patterns which connect to the disk (see text).</figcaption></figure>

<figure><img src="image_4.png" /></figure>

The first types of simulations are identical and/or similar to the ones performed in Masson et al. (2016). They have an initial core mass of $1 M_\odot$, a uniform density profile, and a uniform magnetic field with a mass-to-flux over critical mass-to-flux ratio of 2 or 5. We considered various levels of turbulence, ranging from $M = 0.2$ to $M = 1.2$, different values of $\alpha$ (thermal over gravitational energy) and $\beta$, and different angles $\theta$ between the initial magnetic field and the rotation axis. For the second type of simulation, we considered a massive core of $100 M_\odot$ with a uniform temperature of 20 K. The initial density profile follows $\rho(r) = \rho_c/(1 + (r/r_c)^{-2})$, where $\rho_c \sim 7.7 \times 10^{-18} \text{ g cm}^{-3}$ and the extent of the central plateau is $r_c = 0.02 \text{ pc}$. The initial core radius is $r_0 = 0.2 \text{ pc}$. Radiative transfer is properly accounted for in the simulations, as in Commerçon et al. (2011a), and takes into account the feedback from protostellar luminosity using pre-main sequence evolution models (Hosokawa et al. 2010) attached to sink particles. The coarser grid resolution is $64^3$, and we allow for nine additional levels of refinement, which gives a minimum resolution of 5 au (the sink accretion radius is then of 20 au).

## 4.2. Results

The disk radius is defined according to the criteria described in Joos et al. (2012). We first perform an azimuthal average of the rotation, radial velocity, and sound speed. We then select the rings for which both the radial velocity and the sound speed are smaller than 50% of the rotation velocity. Figure 2 displays the ratio of the disk radius measured in the simulations at different times over the radius inferred from Equation (13), as a function of the total (star+disk) system mass. The left panel corresponds to $1 M_\odot$ mass cores and the right one to $100 M_\odot$ mass cores. The central mass corresponds to either the mass of the first Larson core in the low-mass models or the mass of the sink particle in the high-mass ones. For these latter simulations, we also considered a value of 0.3 G for $B_z$ in Equation (13), as mentioned in Section 2. As seen in the figures, the agreement

between the theoretical predictions and the simulation is globally quite satisfactory. Most of the points lie between 0.5 and 2, indicating that our theoretical estimate $r_{d,AD}$, given by Equation (13), agrees within less than a factor 2 with the numerical results.

For some simulations, notably for the low-mass cores, we see a sudden and steep increase of the radius by a factor 2–3 above some mass. We verified that this occurs when the disk and stellar mass reaches about 30%–50% of the prestellar mass, depending on the various parameters, and the estimated Toomre parameter becomes much smaller than unity. This behavior thus corresponds to the nonlinear development of spiral patterns, which connect to the disk (see Figures 8 and 11 of Masson et al. 2016), making the definition of a disk radius rather ambiguous. The dynamics of these patterns clearly differs from the one of an axisymmetric disk.

## 5. CONCLUSION

In this paper we proposed simple analytical arguments for the formation of early circumstellar disks in collapsing magnetized cores, suggesting that the disks are self-regulated by the magnetic braking and the ambipolar diffusion. The disk radius estimates derived from the theory have been compared to the values obtained from a series of non-ideal MHD simulations, covering a large range of masses, turbulent support, geometrical configurations and magnetic intensities. The comparisons show an agreement between the theoretical and numerical results, within less than a factor 2. The most striking result is the weak dependence of the disk size on the core mass, the intensity of the field, or the level of turbulence in the core, suggesting small variations between class-0 disk sizes under different environments. Clearly, further observations should be able to probe this prediction.

We thank the anonymous referee for a helpful report. This research has received funding from the European Research
Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 247060 and No. 306483). We acknowledge financial support from "Programme National de Physique Stellaire" (PNPS) of CNRS/INSU, France. We thank Rolf Kuiper for providing tabulated PMS evolution tracks.

## REFERENCES

* Allen, A., Li, Z.-Y., & Shu, F. 2003, ApJ, 599, 363

* Belloche, A. 2013, in EAS Publ. Ser. 62, Role and Mechanisms of Angular Momentum Transport During the Formation and Early Evolution of Stars, ed. P. Hennebelle & C. Charbonnel (Les Ulis: EDP Sciences), 25

* Commerçon, B., Hennebelle, P., & Henning, T. 2011a, ApJL, 742, L9

* Dapp, W., & Basu, S. 2010, A&A, 521, 56

* Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

* Galli, D., Lizano, S., Shu, F. H., & Allen, A. 2006, ApJ, 647, 374

* Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528

* Hennebelle, P., & Ciardi, A. 2009, A&A, 506, L29

* Hennebelle, P., Commerçon, B., Joos, M., et al. 2011, A&A, 528, A72

* Hennebelle, P., & Fromang, S. 2008, A&A, 477, 9

* Hennebelle, P., Whitworth, A. P., Cha, S.-H., & Goodwin, S. P. 2004, MNRAS, 348, 687

* Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478

* Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, 128

* Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S. 2013, A&A, 554, 17

* Krasnopolsky, R., Li, Z.-Y., Shang, H., & Zhao, B. 2012, ApJ, 757, 77

* Larson, R. B. 1969, MNRAS, 145, 271

* Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180

* Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

* Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32

* Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24

* Maury, A. J., André, P., Hennebelle, P., et al. 2010, A&A, 512, A40

* Mellon, R. R., & Li, Z. 2008, ApJ, 681, 1356

* Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199

* Price, D. J., & Bate, M. R. 2007, Ap&SS, 311, 75

* Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2012, ApJ, 747, 21

* Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2012, MNRAS, 423, L40

* Shu, F. H. 1977, ApJ, 214, 488

* Teyssier, R. 2002, A&A, 385, 337

* Tobin, J. J., Looney, L. W., Wilner, D. J., et al. 2015, ApJ, 805, 125

* Tomida, K., Machida, M. N., Saigo, K., Tomisaka, K., & Matsumoto, T. 2010, ApJL, 725, L239

* Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117

* Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, MNRAS, 452, 278

* Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037