Effects of radiation transfer on the structure of selfgravitating disks, their fragmentation and evolution of the fragments
Abstract
We investigate structure of selfgravitating disks, their fragmentation and evolution of the fragments (the clumps) using both analytic approach and threedimensional radiation hydrodynamics simulations starting from molecular cores. The simulations show that nonlocal radiative transfer determines disk temperature. We find the disk structure is well described by an analytical model of quasisteady selfgravitating disk with radial radiative transfer. Because the radiative process is not local and radiation from the interstellar medium cannot be ignored, the local radiative cooling would not be balanced with the viscous heating in a massive disk around a low mass star. In our simulations, there are cases in which the disk does not fragment even though it satisfies the fragmentation criterion based on disk cooling time ( and ). This indicates that at least the criterion is not sufficient condition for fragmentation. We determine the parameter range for the host cloud core in which disk fragmentation occurs. In addition, we show that the temperature evolution of the center of the clump is close to that of typical first cores and the minimum initial mass of clumps to be about a few Jupiter mass.
keywords:
star formation – circumstellar disk – methods: hydrodynamics – smoothed particle hydrodynamics – protoplanetary disk – planet formation1 Introduction
Stars form in gravitationally collapsing molecular cloud cores. Since molecular cloud cores have angular momentum (Goodman et al., 1993; Caselli et al., 2002), a circumstellar disk forms around the protostar. According to recent theoretical studies of the gravitational collapse of molecular cloud cores, the protostar is surrounded by a circumstellar disk early in its evolutionary phase in the case without magnetic field, (Bate, 1998; Machida, Inutsuka & Matsumoto, 2010; Tsukamoto & Machida, 2011) and with magnetic field (Inutsuka, Machida & Matsumoto, 2010; Machida, Inutsuka & Matsumoto, 2011).
As Inutsuka, Machida & Matsumoto (2010) pointed out, a circumstellar disk can be gravitationally unstable during its early evolution. When the protostar forms, its mass is only to . In contrast, the firstcore that is the precursor to the circumstellar disk (Saigo & Tomisaka, 2006; Machida, Inutsuka & Matsumoto, 2010) has an initial mass of . Therefore, in the formation epoch of the protostar, the disk has a mass that is significantly greater than that of the protostar. Observations also suggest that massive disks can form in the main accretion phase (Enoch et al., 2009).
In such massive disks, gravitational instability (GI) can develop. The Toomre’s value,
(1) 
is a wellknown indicator for GI (Toomre, 1964). When , the disk is gravitationally unstable against nonaxisymmetric perturbations and develops spiral arms (Laughlin & Bodenheimer, 1994). The spiral arms readjust the mass and angular momentum of the disk, promoting mass accretion onto the protostar. They also heat the disk gas. By readjusting the surface density and the disk heating, the value increases, and the disk is stabilized. This selfregulative nature is an important aspect of GI. To maintain GI, additional physical processes such as radiative cooling or mass accretion from the envelope are necessary. With these effects, the disk can also fragment and gravitationally bound gaseous objects (which we call “clumps”) form.
Disk fragmentation is a candidate mechanism for the formation of wideorbit planets (Vorobyov & Basu, 2010a; Machida, Inutsuka & Matsumoto, 2011; Vorobyov, DeSouza & Basu, 2013). Wideorbit planets are planets separated from the central star by more than 10 AU (Marois et al., 2008; Thalmann et al., 2009; Lagrange et al., 2009; Marois et al., 2010; Lafrenière, Jayawardhana & van Kerkwijk, 2010). Furthermore, it has been suggested that disk fragmentation can explain the formation of brown dwarfs (Stamatellos & Whitworth, 2009a; Stamatellos et al., 2011) and multiple stellar systems (Machida et al., 2008; Kratter et al., 2010).
Effects of radiative cooling on GI and disk fragmentation has been investigated in Gammie (2001) with two dimensional local shearing box simulations. To model radiative cooling, he employed a simplified cooling law (see right hand side of equation (5)). In his simulations, the disk is initially unstable against GI () and GI immediately develops and heats the disk. When radiative cooling is not so strong, it balances the heating by GI, and the disk settles into a quasisteady state. This quasisteady state also realizes in three dimensional global disk simulations (Lodato & Rice, 2004). On the other hand, when radiative cooling is sufficiently strong, such a quasisteady state can not be realized and the disk fragments. In the simulations of Gammie (2001), disk fragmentation occurred when disk cooling time is comparable to the orbital timescale, (or ). The fragmentation condition are also confirmed by three dimensional global disk simulations (Rice et al., 2003) but the simplified cooling laws are also used.
Since Gammie (2001) and Lodato & Rice (2004) employed the simplified cooling law, they implicitly assumed that radiation just acts as the cooling process to decrease the local gas energy (we call this the assumption of local radiative cooling). However, this assumption is not necessarily true. For example, irradiation from the protostar and that from the inner disk region can heat the disk. Thus, in reality, incoming radiation flux, in addition to outgoing radiation flux and local GI heating, should be considered. Furthermore, the interstellar medium has a typical temperature of about 10 K and it is almost impossible to cool the disk gas below 10 K because of radiation flux from the ambient interstellar medium. Therefore, it is not clear that, in realistic situations, the local balance between radiative cooling and viscous heating is achieved or not as in Gammie (2001). Whether the balance is achieved is very important because the structure of quasisteady self gravitating disk can be determined by the energy balance.
Another important issue is applicability of the fragmentation criterion found by Gammie (2001), and , in realistic situations. Gammie pointed out that, in realistic systems, fragmentation realizes when the external irradiation quickly diminishes and when the gas quickly cools. Thus, he regarded the fragmentation criterion can be applicable in very limiting cases.
On the other hand, Rice et al. (2003) interpreted this criterion as “almost isothermal conditions are necessary for fragmentation”. When the cooling timescale of the disk is small, the gas evolves isothermally during GI growing and pressure repulsion becomes weak compared to the adiabatic evolution case (or inefficient cooling case). Rice et al. (2003) suggested that such almost isothermal evolution in nonlinear evolution phase of GI is necessary for fragmentation. According to this interpretation, how the disk becomes or the energy balance within it is not important because whether fragmentation realize or not depends on the thermal behavior of the gas in the nonlinear evolution of GI. The criterion seems to be regarded as a necessary condition in this interpretation.
In the previous works using analytic approach, however, the criterion is used as if a sufficient condition for fragmentation (Rafikov, 2005; Forgan & Rice, 2011; Kratter & MurrayClay, 2011; Forgan & Rice, 2013). As just described above, the interpretation of the criterion is ambiguous and its applicability in the realistic situation is still not clear.
Another mechanism that makes a disk gravitationally unstable is mass loading from the envelope (e.g., Matsumoto & Hanawa, 2003; Inutsuka, Machida & Matsumoto, 2010; Kratter et al., 2010; Vorobyov & Basu, 2010a; Machida, Inutsuka & Matsumoto, 2011; Tsukamoto & Machida, 2011; Bate, 2011; Stamatellos, Whitworth & Hubber, 2011; Tsukamoto, Machida & Inutsuka, 2013; Tsukamoto & Machida, 2013). This mechanism is efficient in the early phase of disk evolution during which the protostar and disk are embedded in a massive envelope. With mass accretion from the envelope, the value of the disk decreases due to the increase of the surface density. When the mass accretion onto the disk is sufficiently high, the disk eventually fragments, and planetary mass clumps form. In this process, whether fragmentation occurs depends on envelope (or cloud core) parameters such as thermal energy or rotational energy. The parameter range in which fragmentation occurs has already been investigated in our previous studies (Tsukamoto & Machida, 2011, 2013); however, in those studies, the effects of radiative transfer were approximated. Radiative transfer could also play important roles during the early phase of disk evolution; therefore, a parameter survey of disk fragmentation with radiation hydrodynamics simulation is necessary.
Initial properties and the evolution of fragments (or clumps) are other important issues. The ultimate fate of clumps depends on their orbital and internal evolutions. For example, if the radial migration timescale is very short, clumps quickly accrete onto the protostar and disappear. If the radial migration timescale is sufficiently long and the clumps survive in the disk maintaining its mass in the range of planetary masses (), the clumps evolve into the wide orbit planets (Vorobyov & Basu, 2010a; Machida, Inutsuka & Matsumoto, 2011). If the mass of the clumps quickly increases, disk fragmentation can be regarded as the formation process for brown dwarfs and binary or multiple stellar systems (Stamatellos & Whitworth, 2009a).
In spite of its importance, only a few studies have focused on the orbital and internal evolution of clumps in circumstellar disks. Baruteau, Meru & Paardekooper (2011) investigated the orbital evolution of massive planets in a gravitationally unstable disk under the assumption of local radiative cooling. They showed that massive planets rapidly migrate inward in a type I migration timescale (Tanaka, Takeuchi & Ward, 2002). By adopting an analytical approach or threedimensional smoothed particle hydrodynamics (SPH) simulation, Nayakshin (2010) and Galvagni et al. (2012) investigated the internal evolution and the collapse timescale of an isolated clump. The gravitational collapse (the second collapse) occurs in the clump when the central temperature reaches K, at which the dissociation of molecular hydrogen begins, and the gas pressure can no longer support the clump against its self gravity. These authors suggested that the timescale for the second collapse after clump formation is in the range of several thousand years to several years. However, they ignored further mass accretion from the disk onto the clumps.
Tsukamoto, Machida & Inutsuka (2013) investigated the orbital and internal evolution of clumps permitting realistic gas accretion from the disk onto clumps with threedimensional radiation hydrodynamics simulation. They showed that although most of the clumps rapidly migrate inward and finally fall onto the protostar, a few clumps can remain in the disk. They also showed that clumps are convectionally stable (). The central density and temperature of a surviving clump rapidly increase, and the clump undergoes a second collapse within – years after its formation. However, in this works, only one simulation was performed. The evolution process of clumps may change under different initial conditions. For example, the central entropy of the clump would become small when fragmentation occurs in the outer relatively cold region of the disk. Since the central entropy of the clump determines its initial mass and radius of the clump, fragmentation in different disk environments changes the nature of the clumps. It is also expected that the mass accretion from the disk onto clumps becomes small when fragmentation occurs in the outer region where the disk surface density is small at .
In this study, we investigate the structure of selfgravitating disks, fragmentation of the disk, and evolution of clumps using both an analytic approach and radiation hydrodynamics simulations starting from molecular cloud cores. This paper is organized as follows. In §2, we analytically derive the structure of the selfgravitating steady disk using various energy balance relations. It would be insightful to understand the general structure of selfgravitating disks. The numerical method and initial conditions for the simulations are described in §3. Results and discussions are divided into two parts (§4 and §5). In §4, we mainly investigate the structure of a selfgravitating disk that formed in the radiation hydrodynamics simulations. The simulation results for disk evolution are presented in §4.1. In this section, we show that radiation can transfer energy within the disk and can heat the outer region of the disk. Therefore, the assumption of local radiative cooling is not valid in the disk. In §4.2, we discuss the consequences of the results obtained from our numerical simulations and their implications in interpreting previous work. In §5, we investigate disk fragmentation and the evolution of clumps. In §5.1, simulation results are shown. The parameter range of the molecular cloud core in which fragmentation occurs is determined, and the orbital and internal evolutions of clumps are investigated. In §5.2, we discuss the results reported in §5.1 and derive the minimum initial mass of clumps. Finally, we summarize our results and discuss future perspectives in §6.
2 Structures of selfgravitating steady disks with various energy balance equations
The analytical description of the radial structure of a selfgravitating steady disk is useful for interpreting the simulation results shown in this study and in previous studies. In this section, we derive radial profiles for physical quantities in a selfgravitating circumstellar disk using various energy balance equations. We assume that the disk (1) is steady (), (2) can be described by the viscous accretion disk model (Shakura & Sunyaev, 1973) and (3) is in the quasi steadystate against gravitational instability ().
A steady viscous accretion disk should satisfy the following equations,
(2) 
where, , and is the kinematic viscosity. The for a selfgravitating disk becomes
(3) 
Here and in the following, we assume the epicycle frequency scales the same as the angular velocity . We assume that physical quantities obey the single power law in radius, ,, , and . Then, (2) and (3) lead to
(4) 
From these equations, we can determine the profile of a steady selfgravitating disk by specifying a rotation profile and an energy balance equation.
Equations (4) shows that only a globally isothermal disk can achieve a globally constant . In other words, if the disk has a radial temperature profile, then it also inevitably has a radial profile for . Therefore, all selfgravitating disk models that have a radially constant and a radially dependent temperature profile violate the assumption of steady state.
2.1 Disk structure with local cooling law
First, we investigate the steadystate structure of a disk with local cooling law (Gammie, 2001), which is often used in global disk simulations in the context of gravitational instabilities or disk fragmentation (e.g., Rice, Lodato & Armitage, 2005; Lodato & Rice, 2005; Paardekooper, Baruteau & Meru, 2011; Meru & Bate, 2012).
At the steady state with local cooling law, the energy balance is given by relation,
(5) 
This can be rewritten as
(6) 
where, we adopted as in most previous studies and equation (6) gives a disk with constant . Note that there is no reason to expect that is constant.
Combining (2), (3), and (6), we obtain
(7) 
Here, we assumed Keplerian disk, . This cooling law gives a globally isothermal disk and a surface density of . Indeed, Fig. 1 in Baruteau, Meru & Paardekooper (2011) shows that the disk is almost globally isothermal and that the surface density has a profile of in the quasisteady state of their simulations. Their results are consistent with our estimate.
Because the disk is isothermal, we can easily calculate the disk temperature in the quasisteady state (we call this “the equilibrium temperature”) using physical quantities at a radius. For example, we can calculate the temperature of the disk used in Rice, Lodato & Armitage (2005) and Meru & Bate (2012). Since dimensionless disk parameters were used in those studies, we need to convert them into dimensional parameters to derive dimensional disk quantities. If we regard the disk model used in those studies as a disk with , the disk cutoff radius as AU and AU as suggested by Meru & Bate (2012), then the equilibrium temperature can be calculated by substituting the value of the physical quantities into
(8) 
Here we have used , and is the mean molecular weight, is the mass of a hydrogen atom, is the Boltzmann constant, and is the ratio of heat capacities. This equilibrium temperature is too small for a disk temperature around a lowmass star. For example, the temperature of disk at AU is about K if the luminosity of the central star is . Note also that the disk is very compact and condensed ( within AU). This disk has a typical Jeans mass of
(9) 
(see, top panel of Fig. 7). Therefore, the mass of the fragments formed in their simulations would be about . Note that, because the Jeans mass is very small, the resolution requirement for such a disk is very severe and particles are required for an SPH simulation to resolve the Jeans mass and satisfy the resolution requirement suggested by Bate & Burkert (1997) (see, §4.2.4 for more discussion).
2.2 Disk structure with the local energy balance equation
Next, we investigate the disk structure when local viscous heating balances local radiative cooling. This assumption is also often used in the context of disk fragmentation (e.g., Rafikov, 2005; Kratter & MurrayClay, 2011). The local energy balance equation between radiative cooling and viscous heating is given by
(10) 
where is the optical depth and we assume . Here we assume that the disk is optically thick. When the gas temperature below K, the opacity is well described by
(11) 
where is a constant (Semenov et al., 2003). As we will see below, the disks formed in our simulation satisfy the condition, K for AU.
Using (2), (3), and (10), we have
(12) 
where we again assume the disk is Keplerian, . With this energy balance equation, physical quantities have very steep radial profiles. Especially, the radial dependence of the temperature () is significantly steeper than the profile expected from the passively irradiated disk model, (see, Kusaka, Nakano & Hayashi, 1970; Chiang & Goldreich, 1997). Whether such steep profiles actually realize in realistic situations is unclear, and we must investigate the disk profile with threedimensional radiation hydrodynamics simulations.
The exponents in the power laws for the disk profile are determined by (2), (3), the energy balance equation (10), and the rotation profile . Thus, there is no degree of freedom available to change them. Therefore, it is impossible to realize a radially constant in the Keplerian disk with local energy balance equation.
2.3 Disk structure with a given temperature profile
As we will see below, the radiation flux from the inner region can heat the outer region of the disk. Furthermore, in a realistic disk around a protostar, stellar irradiation may have an important effect on the disk temperature. Therefore, it is useful to investigate the steadystate structure of the disk with a given temperature profile.
In a disk whose temperature is passively determined, the temperature profile can be written as
(13) 
instead of using an energy balance equation, such as equation (6) or (10). The value of depends on the disk model. For example, in the irradiated disk model, (see, Kusaka, Nakano & Hayashi, 1970; Chiang & Goldreich, 1997). Using (2), (3), , and , we can derive the structure of the irradiated steady selfgravitating disk,
(14) 
This structure is expected when stellar irradiation directly reaches the disk photosphere. If the radiation is extinct before it reaches the disk photosphere, the temperature profile becomes steeper than that given in (14). With these results in mind, we investigate the disk structure observed in threedimensional radiation hydrodynamical simulations.
3 Numerical Method and Initial Conditions
In the simulations performed for this study, we solved the radiation hydrodynamics equations with the fluxlimited diffusion (FLD) approximation SPH schemes according to Whitehouse & Bate (2004) and Whitehouse, Bate & Monaghan (2005). We adopted the tabulated equation of state used in Tomida et al. (2013), in which the internal degrees of freedom and chemical reactions of seven species are included. The hydrogen and helium mass fractions were assumed to be and , respectively. We used the dust opacity table provided by Semenov et al. (2003) and the gas opacity table from Ferguson et al. (2005).
When the gas density exceeds the threshold density, , a sink particle was introduced. Around this density, the gas temperature reaches the dissociation temperature of molecular hydrogen ( K) and the second collapse begins in the firstcore or in the clumps. Therefore, we can follow the thermal evolution just prior to the second collapse. The sink radius was set to AU.
We adopted an isothermal, uniform and rigidly rotating molecular cloud core as the initial condition. Starting the simulations from molecular cloud cores has several merits for investigating the disk evolution. First, we can investigate the formation and evolution of the disk in a selfconsistent manner. Second, the boundary conditions for radiative transfer can be placed far from the disk photosphere. As Rogers & Wadsley (2011) pointed out, the treatment of the boundary conditions significantly affects disk fragmentation in radiation hydrodynamics simulations. Thus, an appropriate treatment for boundary conditions is crucial to investigate both the temperature structure of the disk and disk fragmentation. For the boundary conditions, in this work, we fixed the gas temperature to be 10 K if the gas density is less than . Thus, the boundary was far from the disk photosphere and did not affect the temperature structure of the disk.
The initial mass and temperature of the cores were 1 and 10 K; therefore, the free parameters of the core were the radius and the angular velocity . The parameters of the initial cloud cores are listed in Table 1. In this table, the values of the initial condition used Tsukamoto, Machida & Inutsuka (2013) (model TMI13) is also shown to allow us to investigate the boundary between fragmentation and nofragmentation models. In the table, and are the ratios of thermal to gravitational energy ( ) and rotational to gravitational energy ( ), respectively. The values for cloud core parameters adopted in this study are consistent with results obtained by recent threedimensional magnetohydrodynamics (MHD) simulations that investigated the evolution of the molecular cloud and involved core formation (Inoue & Inutsuka, 2012). The initial cores were modeled with about 520,000 SPH particles.
Model  R (AU)  Fragmentation  

1  0.58  0.01  No  
2  0.58  0.02  Yes  
3  0.48  0.007  No  
4  0.48  0.01  Yes  
5  0.38  0.007  No  
TMI13  0.38  0.01  Yes  

4 Quasisteady structures of self gravitating disks with radiative transfer
4.1 Simulation results
Overview of disk evolution
We calculated the early evolution of the disk until about 10,000 years after the formation of protostar. Figure 1 shows the evolution of surface density at the center of the cloud core for model 1 in Table 1. In Figs. 1a to 1e, we can clearly see the bimodal structure, i.e., the central pressure supported core (the firstcore) and the disk around it. The figure indicates that, before the second collapse, the disk forms around the central firstcore, and spiral arms develop. These features are commonly seen in the unmagnetized or weakly magnetized molecular cloud cores (Bate, 1998; Tsukamoto & Machida, 2011; Bate, 2011; Machida, Inutsuka & Matsumoto, 2010, 2011; Tsukamoto, Machida & Inutsuka, 2013; Bate, Tricco & Price, 2014). The central temperature of the firstcore exceeds K and the second core (or protostar) forms between Figs. 1e and 1f. The white dot at the center of Figs. 1f to 1i corresponds to the location of the protostar. The disk radius gradually increases by mass accretion from the envelope, and angular momentum is redistributed by the spiral arms.
Figure 2 shows gas temperatures at the center of the cloud at the same epochs as in Fig. 1. Comparing with Fig. 1, Fig. 2 shows that the structure of the gas temperature is more axisymmetric than that of the surface density. Figures 1 and 2 indicate that there is only a very weak correlation between surface density and temperature. If we assume that the heat source for the disk is viscous heating due to the gravitational instability, it is difficult to explain the axisymmetric temperature distribution because gravitational energy mainly converts to thermal energy around the spiral arms and the heating should be localized around the spiral arms, and thus, the temperature structure should trace the surface density structure. Although we can see weak heating by the spiral arms (e.g., in Fig. 2g), the entire temperature structure is almost axisymmetric. Therefore, there must be other heating mechanisms that causes the axisymmetric temperature structure. As we will see below, the radiation flux within the disk determines the axisymmetric component of the gas temperature.
Vertical disk structures
Figure 3 shows the density contour map on the plane at the epoch in Fig. 1i. The red lines show the scale height of disk, , where is the sound velocity on the midplane, and is the angular velocity at . The cyan lines show the height of the vertical photosphere, defined as . The green lines and arrows show the temperature contour and direction of the radiation flux, respectively. The green lines show that, in the outer region, ( AU), the disk is almost vertically isothermal. This feature also occurs in twodimensional radiation hydrodynamics calculations (Yorke, Bodenheimer & Laughlin, 1993). The green arrows around the midplane are almost parallel to the xaxis. This indicates that the radial component of the radiation flux dominates the vertical component in the outer disk region. Therefore, the radial component of radiative transfer play an important role in determining the temperature structure of the disk. We can see weak heating due to the spiral arms at where the green arrows spread out. Although such spiral heating occasionally occurs, it is transient and does not significantly contribute to the temperature in the outer disk region. Actually, as we will see below, the theoretical disk model in which local heating due to GI balances radiative cooling cannot reproduce the radial profiles of this disk (see thick dashed lines in Fig. 4).
The vertical density structure is well described by the vertically isothermal thin disk profile, , although the disk is slightly compressed by the ram pressure of mass accretion from the envelope. The aspect ratio of the disk is about 0.2 at AU. In the inner region, AU, the height of the disk photosphere is greater than the disk scale height. This structure is different from passively irradiated optically thick disk (Kusaka, Nakano & Hayashi, 1970). This is one reason why our disks have a steeper profile (see left middle panel in Fig. 4) than the passively irradiated disk model .
Radial disk structure
Figure 4 shows the azimuthally averaged radial profiles of the disk. Thin solid, dashed, dotted and dasheddotted lines are for profiles at the same epochs as in Figs. 1f, g, h and i, respectively. In the figure, top left panel shows the angular velocity . Because the disk is massive, the rotation profile is not Keplerian (; thick dotted line). We fitted to (thick dotteddashed line) to obtain the power law of the disk with (2) and (3).
The middle left panel shows values calculated from azimuthally averaged quantities,
(15) 
where means the azimuthal average. The panel shows over almost all regions of the disk. This is a general feature of the selfgravitating disks without fragmentation and the disk is in the quasisteady state against gravitational instabilities.
The bottom left panel shows the midplane gas temperature. The thick dashed line shows the theoretically estimated power law under the assumption that local heating balances local cooling (; thick dashed line). This power law is calculated from (2), (3), (10), and . This power law does not agree with the disk structure because the energy balance is not local, and the radial component of the radiative transfer heats the outer disk region. We fitted the temperature to (thick dasheddotted line). The power law can be derived with the diffusion approximation (see §4.2.1 for details).
The right top panel shows the profile of surface density . The thick solid lines show the power law that is predicted from (2), (3), , and . The thick dashed line is the power law predicted from equations (2), (3), , and equation (10). In the outer region with AU, the surface density profile is consistent with the thick solid line. However, in the region of AU, the profile does not agree with the thick solid line because is not constant (see left bottom panel), and the assumption of does not apply in this region. The thickdashed line does not agree with the profile over the entire region.
The right middle panel shows the profile for parameter. parameter at each radius were calculated from
(16) 
where is the component of the viscous stress tensor. The parameter is averaged over where we took years. The stress tensor associated with the gravitational field is given by
(17) 
where and are the radial and azimuthal components of the gravitational force, respectively. The stress tensor associated with velocity field fluctuations or Reynolds stress is calculated by
(18) 
Here, we define a velocity fluctuation by . The total stress tensor is calculated from the sum of these tensors, .
In the profile, corresponds to the epochs of Figs. 1f (solid), g (dashed), and h (dotted). The profile can be described well by the theoretical estimate from (2), (3), , and (; thick solid line). As expected, was not constant. This is a general feature of a nonisothermal selfgravitating disk. Note that the model in which local viscous heating balances local radiating cooling predicts a very steep radial profile (thick dashed line) and fails to explain the profile obtained from the numerical simulation.
The bottom right panel shows the vertical optical depth, . The optical depth was calculated from . In almost all regions of the disk, the vertical optical depth is greater than unity, . Again, the thick solid line describes the radial profile very well.
Cooling time of disk
In Fig. 5, we show the azimuthally averaged cooling time of the disk normalized by orbital frequency, . Here, cooling time is calculated as
(19) 
where is midplane temperature and . We average above in the azimuthal direction. In outer region of the disk, the cooling time satisfies and it is expected that the gas behaves almost isothermally during GI growing. Note that, for our disk, this cooling time is not the timescale on which disk temperature decreases because the temperature of our disk is largely affected by radial radiative transfer. However, we use this cooling time as an indicator how the disk gas behaves against the dynamical compression by GI. When the cooling time is sufficiently small (large), the gas evolves isothermally (adiabatically) during GI growing. As we have seen in Fig. 4, value of the disk is close to 1. Thus, our disk apparently satisfies the fragmentation criterion suggested by Gammie (2001) and Rice et al. (2003). However, the disk does not fragment. This shows that the fragmentation criterion based on disk cooling time is not sufficient condition for fragmentation.
4.2 Discussion
Temperature profile determined by radial radiative transfer in an optically thick disk
The results in §4.1 shows that the temperature profile of the disk formed in our simulation is . This is significantly shallower than the disk in which local heating balances local cooling (see, Fig. 4) and steeper than the passively irradiated disk . In this subsection, we analytically derive the profile .
The right bottom panel in Fig. 4 shows that the vertical optical depth is greater than unity over almost the entire region. Thus, we can use the diffusion approximation to derive the disk temperature profile. Under the assumptions that disk is steady, viscous heating is negligible, and that the radiation temperature is equal to the gas temperature, the energy equation can be expressed as
(20) 
Figure 2 and 3 shows that the temperature distribution can be regarded as axisymmetric and vertically isothermal. Thus, we can assume that . Under the assumption that , (20) becomes
(21) 
Around the midplane, the diffusion approximation is justified, and (21) becomes
(22) 
Using the relation , we have
(23) 
Using (2), (3), and (23), we obtain
(24) 
This agrees very well with the simulation results (see, the thin lines (simulation results) and thick solid lines (theoretical estimate) in Fig. 4).
Importance of nonlocal radiative transfer
As pointed out in §1, the locality of radiative cooling has been assumed in many previous simulations (e.g., Rice et al., 2003; Lodato & Rice, 2005; Rice, Lodato & Armitage, 2005; Rice, Forgan & Armitage, 2012; Paardekooper, Baruteau & Meru, 2011) and in many analytical studies (Rafikov, 2005; Clarke, 2009; Kratter & MurrayClay, 2011) to investigate the nature of GI and the conditions of disk fragmentation with radiative cooling. In contrast, our simulations show that radiation can transfer significant energy in the radial direction even in the absence of the stellar irradiation.
Consequently, the temperature profile becomes in our simulation and is significantly shallower than , which is expected under the assumption that local viscous heating balances local radiative cooling. Thus, the local energy balance between radiation and viscous heating is not satisfied in the disk formed in our simulations. This indicates that the assumption of local radiative cooling is not necessarily valid in a massive disk around a low mass star. The radial profiles of the surface density and the viscous parameter are consistent with a selfgravitating quasisteady disk model with a given rotation profile and . Thus, the disk is in the quasisteady state with an energy balance that is not local.
Therefore, we conclude that the local treatment of radiative cooling is not suitable approximation to investigate massive disk around low mass star.
We only calculated disk evolution just about years after protostar formation and the disk and the protostar were still embedded in a massive envelope (Class 0 phase). One might think that nonlocal radiative transfer does not play the role in later evolutionary phase. However, with twodimensional radiation hydrodynamics simulation, Yorke, Bodenheimer & Laughlin (1993) showed that an almost isothermal temperature distribution in vertical direction also realize in the late evolutionary phase of disk at which 95 % of initial cloud core has already accreted onto the protostar or disk. Thus, we expect that nonlocal radiative transfer also plays an important role for the temperature structure of disks in the later evolutionary phase. Note that a large simulation box in which the boundary is far from disk photosphere is crucial to investigate the effects of nonlocal radiative transfer, otherwise radiation flux is easily affected by the boundary conditions.
We could not find evidence that the nonlocal energy transport of GI suggested by Balbus & Papaloizou (1999) plays an important role because GI heating is itself small compared to radiation energy transfer in the outer region of the disk.
Applicability of fragmentation criterion based on disk cooling time
In Fig. 5, we show that the cooling time of disk corresponds to in outer region and our disk apparently satisfies the fragmentation criterion based on disk cooling time ( and ). Because the cooling time is very short compared to the orbital period, the gas evolves almost isothermally during GI growing in outer region. However, the disk does not fragment.
In our disk, the most unstable wavelength of GI is relatively large and the global spiral arms develop. The spiral arms readjust the surface density and the disk is stabilized. Such a readjustment is also observed in previous works (Laughlin & Bodenheimer, 1994; Kratter et al., 2010). Especially, Kratter et al. (2010) shows that the disk fragmentation is suppressed by the readjustment even with isothermal equation of state. Because such a stabilizing process also play the role in the realistic disk, we conclude that at least the fragmentation criterion is not sufficient condition for fragmentation in realistic disk around low mass star.
A significant difference between our disk and disks used in the previous works with local cooling law is the most unstable wavelength of GI. The disks used in previous studies using local cooling law (e.g., Rice et al., 2003; Meru & Bate, 2012) have very small value of the most unstable wavelength and the spiral arms formed in these simulations are not geometrically thick and have large azimuthal mode numbers . It is expected that the efficiency of the readjustment promoted by such spiral arms is small compared to that promoted by global spiral arms which may form in realistic disk.
Is the fragmentation criterion based on local cooling a necessary condition for disk fragmentation ? We think this is not true, either. Although the efficiency of the radiative cooling may affect the evolution of condensations formed by GI, cooling criterion is not necessary for the fragmentation of disk. This can be understood from the previous works that employed the simplified equation of state. For example, Vorobyov & Basu (2006), Inutsuka, Machida & Matsumoto (2010), Tsukamoto & Machida (2011), and Machida & Matsumoto (2011) employed barotropic equation of state and the gas evolves adiabatically when (the condensation exceeds this density at very early phase of its evolution). Even with such a stiff (or adiabatic) equation of state, fragmentation is observed.
Therefore, we conclude that, in general, the fragmentation criterion based on disk cooling time ( and ) is not necessary nor sufficient conditions for fragmentation.
One may think that our results appear very different from the previous results which show the disk fragmentation occurs when it satisfies the fragmentation criterion and does not occur when it does not satisfy the criterion (e.g., Rice, Lodato & Armitage, 2005; Stamatellos & Whitworth, 2008; Stamatellos & Whitworth, 2009a). However, we emphasize that there is no contradiction. Logically speaking, there are four possible cases,

The disk satisfies ( and ) and the disk fragments,

The disk does not satisfies ( and ) and the disk does not fragment,

The disk satisfies ( and ) but the disk does not fragment,

The disk does not satisfies ( and ) but the disk fragments.
The previous works mentioned above show the cases of (i) and (ii). On the other hand, we point out there exist the cases of (iii) and (iv). Therefore, there is no contradiction.
Note, however, that the finite number of examples of (i) and (ii) are not enough to prove the statement that “ALL the disk fragments when it satisfies and ” or “ALL the fragmenting disk must satisfy and ” because we cannot prove nonexistence of the exceptions with finite number of the simulations. On the other hand, just one counterexample is enough to reject these statements. In this paper, we show there exists (iii) case and this is a counterexample for the statement, “ALL the disk fragments when it satisfies and ”. Therefore, the fragmentation criterion is not sufficient condition. On the other hand, the previous works we mentioned above show that there exists (iv) case. This is a counterexample for the statement “ALL the fragmenting disk must satisfy and ”. Therefore, the fragmentation criterion is not necessary condition. Have we emphasize that the discussions above is based on the azimuthally averaged quantities (see §5.1.2).
Resolution consideration for disk fragmentation simulations
Resolution consideration of the simulations performed in this paper
Numerical resolution is considered to be an important factor in simulations of disk fragmentation. Recently, Meru & Bate (2012) argued that the condition for disk fragmentation using local cooling law strongly depends on numerical resolution. Even though our simulation is based on a more realistic radiative transfer method and the disk structure is completely different from theirs, it is important to check whether the numerical resolution in our simulation was sufficient. As Bate & Burkert (1997) discussed, the SPH method cannot correctly simulate the fragmentation phenomenon if the minimum resolvable mass is greater than the Jeans mass , where and are the number of neighbors and mass of the SPH particles, respectively. Here, according to Bate & Burkert (1997), the Jeans mass is given by
(25) 
where is the gas constant. On the other hand, Nelson (2006) suggested that Toomre mass,
(26) 
is more appropriate mass scale for fragmentation in disk and it should be resolved with . Furthermore he also suggested that the disk vertical scale height should be resolved at least , where is the smoothing length. In Fig. 6, we plotted azimuthally averaged Jeans mass (top left), (top right), (bottom left) and (bottom right) for the same epochs as in Fig. 1f (solid), g (dashed), h (dotted), and i (dasheddotted).
The Jeans mass of our disk is about over almost the entire region of the disk; this value is consistent with the initial clump mass that formed in model 2. As shown in top right panel, the Jeans mass is resolved by larger than particles. In our simulations, we adopted the number of neighbors to be ; hence, the mass resolution is about 20 times higher than that required by Bate & Burkert (1997).
The Toomre mass is resolved by larger than particles and this is about 30 times higher than the resolution criterion suggested by Nelson (2006). The bottom right panel shows that the vertical scale height is resolved by larger than and our simulation also satisfies the resolution requirement for vertical scale height. Therefore, we conclude that the numerical resolution in our simulation was sufficient to resolve fragmentation in the disk.
Resolution consideration of the simulations with local cooling law
Why does the disk fragmentation criterion with local cooling law strongly depend on numerical resolution ? To answer this question, we investigate the resolution requirement of the disk used in Rice, Lodato & Armitage (2005) and Meru & Bate (2012) with the quasisteady state structure. In §2, we investigated the steady state of the selfgravitating disk with local cooling law. As we pointed out, the temperatures of the disks used in Rice, Lodato & Armitage (2005); Meru & Bate (2012) are very small in the quasisteady state. Therefore, we can expect that the requirement on numerical resolution for their disks which settled into quasisteady states is very severe. In this subsection, we only consider the resolution requirement suggested by Bate & Burkert (1997).
In Fig. 7, we show the Jeans mass (top) and the required particle number to resolve the Jeans mass, (bottom) as functions of the radius for a disk that has dimensionless parameters of and is in the quasisteady state . Here according to Rice, Forgan & Armitage (2012), we approximated for comparison and we assumed . The Jeans mass is . This corresponds to about 0.01 if we regard the central star mass as . This value of the Jeans mass is very small compared to that in our simulation and to the mass of wide orbit planets ().
The bottom panel in Fig. 7 shows that a particle number of is required to resolve the Jeans mass. Note that, with , the Jeans mass is not resolved over the entire disk region. Meru & Bate (2012) showed that the convergence of the calculation around , and their results are consistent with our estimate. Our estimate is more severe than that of Rice, Forgan & Armitage (2012) in the outer disk region. The difference is because Rice, Forgan & Armitage (2012) assumed instead of assuming the steady state as in (2). Therefore, they derived and found a larger in the outer disk region. Note, however, that such a disk is not realized with local cooling law because the surface density profile also changes because of the angular momentum transfer to realize the quasisteady state structure (see Fig. 1 of Baruteau, Meru & Paardekooper, 2011).
Note also that the resolution requirement is estimated using the initial total disk mass and initial cutoff radii, and . As the disk evolves, gas accretes onto the central star, and the disk radius increases. As a result, the surface density decreases, and the gas cools further to maintain . In a disk with , the Jeans mass is proportional to
(27) 
Therefore, as the surface density decreases, also decreases as and the resolution requirement becomes more severe as the disk evolves. We can ignore the change in due to mass evolution of the central star because .
In an isolated disk with the local cooling law and without a temperature floor, the disk cools infinitely to maintain , and the Jeans mass becomes infinitely small as the disk evolves (or as the surface density decreases). We must carefully monitor the Jeans mass during the simulations.
On the other hand, in a realistic system, there exists an appropriate range of temperature or of the Jeans mass, depending on the system of interest. We emphasize that the nature of the gravitational instability depends not only on the value of but also on the most unstable wavelength of GI, (), and different length scales give strikingly different outcomes for the nonlinear growth phase. For example, the widths of spiral arms and masses of fragments differ for different length scales (compare the spiral pattern in Fig. 1 with Fig. 1 in Rice, Lodato & Armitage, 2005).
We suggest that future selfgravitating disk simulations should use a disk model that has a realistic Jeans length. For example, the irradiated disk model (e.g., Cai et al., 2008) or a locally isothermal model seems to be more suitable for numerical simulations because we can limit the value of the most unstable wavelength to a realistic range. To construct appropriate initial conditions, the profiles derived in §2 would be useful.
5 Disk fragmentation and evolution of clumps
In this section, we investigate fragmentation of a disk and the evolution of clumps after their formation. Understanding the evolution of the clump is crucial to consider the possible final outcomes (wide orbit planets, brown dwarfs or binary/multiple systems) that could result from disk fragmentation.
5.1 Simulation results
Parameter survey for disk fragmentation
As shown in Table 1, we calculated the evolution of the cloud core for five models. In models 1, 3, and 5, fragmentation did not occur in the early phase of disk evolution. On the other hand, in models 2, 4, and TMI13, fragmentation did take place. The computation results are summarized in Fig. 8. In that figure, circles (triangle) indicate models in which fragmentation did not (did) occur. The figure shows how fragmentation correlates with , the mass accretion rate onto the disk, and , the centrifugal radius of the fluid element. Qualitatively, larger mass accretion rates and larger centrifugal radii lead to fragmentation. The boundary that divides fragmentation and nofragmentation models is consistent with that determined by simulations with the barotropic approximation (see., Tsukamoto & Machida, 2011). This result suggests that, in the early evolution of a disk, fragmentation rarely depends on details of radiative transfer or thermodynamics of the gas. This is because disk fragmentation in the early evolutionary phase is mainly determined by mass accretion rate onto the disk and the angular momentum of mass accretion (Matsumoto & Hanawa, 2003; Saigo & Tomisaka, 2006; Vorobyov & Basu, 2010a, b; Inutsuka, Machida & Matsumoto, 2010; Tsukamoto & Machida, 2011)
Difference of the conditions between fragmenting and nonfragmenting disk
In this subsection, we investigate the difference of the conditions of fragmenting and nonfragmenting disk. In Fig. 9, we show the value and cooling time of model 1 (red lines; nonfragmenting disk) and 2 (green lines; fragmenting disk). The epoch of red lines is the same as in Fig. 1h. The epoch of green lines is the same as in Fig. 10a and just prior to disk fragmentation. The solid lines of Fig. 9 show that the azimuthally averaged value and the azimuthally averaged cooling time normalized by orbital period. They show both disks are gravitationally unstable () and their cooling time is sufficiently small (). From the azimuthally averaged quantities, we can not find any significant difference between two disks. But one fragments and the other does not. What is the difference between them ?
The dashed lines in Fig. 9 shows the minimum value, at each radius. This traces the value at the spiral arms. We can see the clear difference between red and green dashed lines. The of fragmenting disk becomes significantly small () in relatively large width ( AU). On the other hand, the profile of in nonfragmenting disk does not show such a dip. We monitored during entire period of evolution of nonfragmenting disk after protostar formation and found that never becomes . Note also that the width of spiral arms of fragmenting disk is smaller than that of nonfragmenting disk (compare Fig. 1h with Fig. 10a).
From these results, we conclude that the value inside the spiral arms (not azimuthally averaged value) and the width of the spiral arms are important quantities for disk fragmentation. The importance of width of spiral arms is also pointed out by Rogers & Wadsley (2012).
Fragmentation of disks and evolution of the clumps
In this subsection, we consider the evolution of the fragments (or clumps) formed by disk fragmentation. We have already investigated the evolution of clumps in Tsukamoto, Machida & Inutsuka (2013). The parameters for the cloud core investigated in Tsukamoto, Machida & Inutsuka (2013) were and (model TMI13 in Table 1). In model TMI13, low and were adopted. This means that the disk formed in model TMI13 is compact. In that case, fragmentation occurred in the region K. Thus, the entropy of the centers of clumps is greater than that of the typical firstcore (see Fig. 4 in Tsukamoto, Machida & Inutsuka, 2013).
When the initial cloud core has a high value of , it is expected that the disk will be more outspread, and clumps will form in lowtemperature regions. Therefore, the entropy of the clumps would become small. The central entropy of clumps is important when we consider the initial mass and evolution of clumps. Smaller entropy leads to smaller clump mass (see eq. (34)) and a smaller clump radius. Thus, investigating clump evolution for larger value of and would be insightful to understand the initial properties and the evolutionary process of clumps in more detail. To do so, we investigate clump evolution in model 2. The cloud core of model 2 has and values greater than those in model TMI13.
Figure 10 shows the evolution of surface density at the center of the cloud core for model 2. Analogous to Fig. 1, Fig. 10a shows bimodal structure (the firstcore and disk around it) in the early phase of disk evolution. The first clump formed in Fig. 10b. We define the epoch of clump formation as the time when its central density exceeds . Within years after the first clump formation, six more clumps formed.
Figure 11 shows the temperature evolution at the center of the cloud core of model 2. As expected, disk fragmentation and clump formation occur in the very lowtemperature region of K. After disk fragmentation, the central temperature of the clump quickly increases because of compressional heating.
Figure 12 shows the orbits of clumps after clump formation (solid lines) and the fluid elements at the centers of the clumps prior to clump formation (dashed lines). Four clumps, shown with red, green, blue, and cyan lines, finally accreted onto the central firstcore or protostar. On the other hand, the remaining three clumps, shown with magenta, purple, and orange lines merged into one clump and became a secondary protostar (see, Fig. 10i).
Figure 13 shows the temperature evolution of the centers of the clumps as a function of density (solid lines). The figure also shows the temperature evolution of fluid elements at the centers of the clumps prior to clump formation (dashed lines). In this figure, a typical temperature evolution with the barotropic equation of state, which is designed to mimic the temperature evolution of the center of the firstcore (Masunaga & Inutsuka, 2000), is also shown for comparison.
Until the density reaches , the temperature is almost isothermal because there is no significant heating by spiral arms. On the other hand, above , the gas evolution becomes adiabatic. The interesting finding is that the thermal evolution of the central temperature is almost consistent with the evolution in the barotropic approximation (black solid line). Because clumps form in the disk, we might expect that fluid elements at the centers of clumps have a good chance of decreasing their entropy by radiative cooling compared with that of the firstcore, which is surrounded by a spherically symmetric massive envelope. However, our results show that centers of clumps do not efficiently cool before they become optically thick. As a result, the evolution of fluid elements follows the line of the barotropic equation of state.
Figure 14 shows how the mass of clumps evolves over time. Here, we define the mass of a clump, , so that it satisfies
(28) 
where and represent the pressure, density, and cumulative mass at , respectively, and is the radius from the clump center. Note that (28) is identical to the Virial theorem when the surface pressure is negligible, and the clump is in hydrostatic equilibrium. Just after clump formation, the mass of each clump is a few Jupiter mass. This mass is slightly smaller than that of the typical clump formed in Tsukamoto, Machida & Inutsuka (2013) because the central entropy of the clump formed in model 2 is smaller than that in Tsukamoto, Machida & Inutsuka (2013). The mass of the clump indicated by orange in Fig. 14 quickly increased at years. This is because the clump encountered a disk spiral arm and gained mass from it. The clumps indicated by purple and orange merged into the clump indicated by magenta. Even with this coalescence, the mass of the clump indicated by magenta did not increase significantly. This is because the clumps indicated by purple and cyan were tidally destroyed and formed circumclump disk around the clump indicated by magenta. Therefore, the mass of the destroyed clumps did not directly accrete onto the clump indicated by magenta.
In the clump indicated by magenta, the second collapse occurred about 3000 years after its formation. This timescale is consistent with our previous results (Tsukamoto, Machida & Inutsuka, 2013). Note that the second collapse occurs at a slightly smaller mass () compared to that in our previous results. This is because the central entropy of the clump is smaller than that of the clumps in Tsukamoto, Machida & Inutsuka (2013). The smaller central entropy induces the second collapse at a smaller clump mass (see (3) in Tsukamoto, Machida & Inutsuka, 2013).
5.2 Discussion
Minimum central entropy of clumps
An important finding in §5.1.3 is that the evolutionary paths of the central temperatures of clumps on the plane are close to that of the firstcore even though the clumps form in the disk. This implies that the clump evolution after the fragmentation can be approximately regarded as sphericallysymmetric. Therefore, in the same manner as in the case of the firstcore, we can describe the evolution of the central entropy of the clump.
Following Masunaga & Inutsuka (1999) and Omukai et al. (2005), the minimum entropy is calculated as follows. We assume the quasiadiabatic evolution begins when the clump becomes opaque. If we assume the clump radius is the Jeans length , the clump becomes opaque when
(29) 
is satisfied.
Because we are interested in the minimum central entropy, we neglect heating by the irradiation. In such a case, the clump evolves according to the energy balance between optically thin radiative cooling and compressional heating until it becomes opaque. The energy balance can be written as
(30) 
where is the freefall timescale. We assume the opacity is given as
(31) 
Using equations (29), (30), and (31), the density and temperature where the adiabatic evolution begins, and are given as
(32)  
(33) 
The minimum value of entropy at the center of the clump is given by these values.
Note that the minimum entropy does not depend on when the ratio of heat capacities is equal to because (Omukai et al., 2005). Therefore, the evolution of the central region of the clump follows the same track in the plane irrespective of different . The property is the same as in the evolution of the protostars with various metallicities shown by Omukai et al. (2005).
Minimum initial mass of the clump
Based on the discussions in §5.2.1, we can evaluate the minimum initial mass of a clump by assuming that the clump can be described by a polytropic sphere (as we have shown in Tsukamoto, Machida & Inutsuka, 2013, the clump can be well described with polytropic sphere with ). the mass of the clump with polytropic index is calculated by
(34) 
where , , and are the central temperature, Boltzmann constant, and mean molecular weight, respectively. In the estimate, the heat capacity at constant volume and the ratio of heat capacities, () are assumed to be constant for simplicity. The estimated initial mass is a few Jupiter mass, which is consistent with the initial mass of the clumps formed in our simulation (see, Fig. 14).
Comparison with previous work
There are few studies about the evolution of the clumps with realistic accretion onto them and sufficient numerical resolution to resolve the central structure of clump (Stamatellos & Whitworth, 2009b; Tsukamoto, Machida & Inutsuka, 2013). In this subsection, we compare our results with those in Stamatellos & Whitworth (2009b).
Stamatellos & Whitworth (2009b) investigated the evolution of the clumps with threedimensional simulation starting from a massive isolated disk. In the disk, fragmentation immediately occurs and the clumps form. They show that the clump undergoes the second collapse when its mass reaches and the timescale for the second collapse is several thousand years. They also points out that thermal evolution of the clump is consistent with the evolution of firstcore (Masunaga & Inutsuka, 2000). The evolution process of the clumps formed in our simulation is consistent with that found in their simulations.
However, there are interesting differences between their results and ours. In their simulations, most of the clumps remained in the disk without falling onto the central protostar. On the other hand, in our simulation, a large fraction of clumps fell onto the central star and disappeared. This difference may come from the fact that the initial disk of Stamatellos & Whitworth (2009b) is very massive () and unstable (). In such a massive disk, fragmentation immediately occurs and many clumps simultaneously form. Thus, many clumps reside in the disk at the same epoch and interact each other. On the other hand, in our simulation, only two or three clumps simultaneously exist in the disk at the same epoch because the disk is not so massive. The difference in the orbital evolution of the clump may come from the difference of the disk. Orbital evolutions of clumps depend more strongly on the condition of disk and further studies are needed on this issue.
Possible evolutionary path to realize the small mass clump
Our simulation results have shown that the central entropy of a clump cannot become significantly smaller than that of the first core. However, there is a possible path that can lead to a clump with an central entropy smaller than that of the first core. The above discussion relies on the clump evolution after disk fragmentation. Thus, the value of in the disk has already become . Once gravitational instabilities turns on, the surface density decreases by mass and angular momentum transfer by the spiral arms, or disk fragmentation occurs. As shown in §5.1.3, the rapid evolution of the clump after fragmentation prevents temperature evolution which leads to smaller central entropy than that of the first core.
However, when is greater than unity, the disk surface density can monotonically increase by infall from the envelope without mass and angular momentum redistribution or fragmentation by GI in the disk. The maximum midplane density of a gravitationally stable Keplerian disk () can be calculated from the condition as
(35) 
where, is the mass of the central star. The maximum midplane density is solely determined by the angular velocity.
This indicates that a high midplane density can be achieved in the inner region of the disk. If the disk with such a high midplane density at the inner region and low temperature somehow fragments, a clump with small central entropy could can be created. For example, if a disk with a temperature of 50 K at 10 AU around fragments, adiabatic contraction of the clump begins from and K. In such a case, the initial mass of clump is . However, the realization of such a disk seems to require the inclusion of magnetic field (Inutsuka, Machida & Matsumoto, 2010; Machida, Inutsuka & Matsumoto, 2011), which is beyond the scope of this paper.
6 Summary and future study
6.1 Summary
In this paper, we investigate the structure of selfgravitating disks, their fragmentation and the evolution of the fragments using an analytic approach and threedimensional radiation hydrodynamics simulations. First, we analytically derive the quasisteady structures of selfgravitating disks with various energy equations. We show that local cooling law, which has been widely used in previous studies (e.g., Rice, Lodato & Armitage, 2005; Lodato & Rice, 2005; Paardekooper, Baruteau & Meru, 2011; Meru & Bate, 2012), describe a globally isothermal disk () in the quasisteady state. In addition, we point out that Jeans mass of the disk with local cooling law and the fiducial disk model ( and ) becomes very small and resolution requirement for the simulation is exceedingly severe. We show that at least particles are required to resolve Jeans mass of the disk. We also point out that such a small Jeans mass (, if we regard the mass of the central star as ) are not attainable in a realistic disk around low mass star since it requires unrealistically low temperature.
We also investigated the quasisteady structure of the disk in which radiative cooling locally balances viscous heating and showed that it has very steep radial profiles (e.g., for Keplerian disk) in the quasisteady state. To investigate whether such a steep profile can be realized (in other words, whether the local balance between radiative cooling and viscous heating can be realized), we conducted threedimensional radiation hydrodynamics simulations. We found that the disk does not have the radial profile expected from the local energy balance and that the disk temperature is nonlocally determined by radiative transfer. The radial profile of the disk temperature was and this scaling can be derived analytically by applying diffusion approximation. Because the temperature of the disk is determined by radial radiative transfer within the disk, radial radiation transport is crucial for outer regions of the disk. Thus, we conclude that the description only with local radiative cooling is not viable in massive disks around low mass stars
The disk formed in our simulation satisfies the fragmentation criterion based on disk cooling time ( and ). However, the fragmentation is not observed in the disk. Therefore, the fragmentation criterion is not sufficient condition for fragmentation. Further studies for more accurate criterion is needed.
Mass accretion from the envelope is another process that possibly drives disk fragmentation. By a parameter survey with radiation hydrodynamics simulations starting from the molecular cloud core, we determined the parameter range of the host cloud core needed for the gravitational fragmentation. The parameter range is in agreement with that obtained in our previous study based on a simplified equation of state (Tsukamoto & Machida, 2011). Therefore, we conclude detailed treatment of radiative transfer is not crucial for disk fragmentation driven by mass accretion from the envelope.
We also investigated the internal evolution of fragments (clumps) formed in extended cold disks ( K in the outer region). Even in such a disk, the central temperature of a clump does not sufficiently cool to have smaller central entropy than that of the first core. Therefore, there is a minimum value of the central entropy. Using this value for the central entropy, we derived the minimum initial mass of the clumps to be about a few Jupiter mass. This is consistent with the initial mass formed in our simulations.
6.2 Future study
In our simulations, the flux limited diffusion (FLD) approximation was adopted. It is well known that the FLD approximation does not behave well in optically thin region. Although the disk formed in our simulation can be regarded as optically thick in almost all regions,
it is possible that the temperature profile changes with more realistic radiative transfer methods because the simulation box includes the optically thin region and the radiation flows along a roundabout path in the FLD approximation. Although the difference is expected to be about a factor of a few (see comparison tests of radiative transfer schemes; e.g., Pascucci et al., 2004; Kuiper & Klessen, 2013) , it is important to perform simulations with more sophisticated radiative transfer method. We plan to perform simulations with more sophisticated radiative transfer schemes to confirm our results for the disk structure of a massive disk around a low mass star.
Another important issue, which is not discussed in this paper, is the effects of magnetic fields. Magnetic fields play important roles in the formation and early evolution of circumstellar disks due to their efficient transfer of angular momentum and formation of outflows (see, Inutsuka, 2012, and references therein). In particular, magnetic fields would change the parameter range of fragmentation that is derived in §5.1.1. It would also be important for orbital and internal evolutions of clumps. The disk is well coupled to magnetic fields because magnetic diffusion is not significant in the range of disk densities. On the other hand, the magnetic diffusion is effective inside clumps because of their high density. Thus, the magnetic field and clumps would be decoupled. In our future study, we will investigate the effects of magnetic fields using numerical methods for smoothed particle magnetohydrodynamics (SPMHD) as described in Iwasaki & Inutsuka (2011) and Tsukamoto, Iwasaki & Inutsuka (2013).
Acknowledgments
We thank K. Iwasaki, T. Matsumoto, H. Kobayashi, M. Ogihara, K. Tanaka and T. Muto for their fruitful discussions. We also thank K. Tomida and Y. Hori to provide their EOS table to us. We also thank the referee, Dr. D. Stamatellos for insightful comments. The snapshots were produced by SPLASH (Price, 2007). The computations were performed on a parallel computer, XC30 system at CfCA of NAOJ and SR16000 at YITP in Kyoto University. Y.T. is financially supported by Research Fellowships of JSPS for Young Scientists.
References
 Balbus S. A., Papaloizou J. C. B., 1999, ApJ, 521, 650
 Baruteau C., Meru F., Paardekooper S.J., 2011, MNRAS, 416, 1971
 Bate M. R., 1998, ApJ, 508, L95
 —, 2011, MNRAS, 417, 2036
 Bate M. R., Burkert A., 1997, MNRAS, 288, 1060
 Bate M. R., Tricco T. S., Price D. J., 2014, MNRAS, 437, 77
 Cai K., Durisen R. H., Boley A. C., Pickett M. K., Mejía A. C., 2008, ApJ, 673, 1138
 Caselli P., Benson P. J., Myers P. C., Tafalla M., 2002, ApJ, 572, 238
 Chiang E. I., Goldreich P., 1997, ApJ, 490, 368
 Clarke C. J., 2009, MNRAS, 396, 1066
 Enoch M. L., Corder S., Dunham M. M., Duchêne G., 2009, ApJ, 707, 103
 Ferguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik J. G., Hauschildt P. H., HeffnerWong A., Tamanai A., 2005, ApJ, 623, 585
 Forgan D., Rice K., 2011, MNRAS, 417, 1928
 —, 2013, MNRAS, 432, 3168
 Galvagni M., Hayfield T., Boley A., Mayer L., Roškar R., Saha P., 2012, MNRAS, 427, 1725
 Gammie C. F., 2001, ApJ, 553, 174
 Goodman A. A., Benson P. J., Fuller G. A., Myers P. C., 1993, ApJ, 406, 528
 Inoue T., Inutsuka S., 2012, ApJ, 759, 35
 Inutsuka S., Machida M. N., Matsumoto T., 2010, ApJ, 718, L58
 Inutsuka S.i., 2012, Progress of Theoretical and Experimental Physics, 2012, 010000
 Iwasaki K., Inutsuka S., 2011, MNRAS, 418, 1668
 Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 708, 1585
 Kratter K. M., MurrayClay R. A., 2011, ApJ, 740, 1
 Kuiper R., Klessen R. S., 2013, A&A, 555, A7
 Kusaka T., Nakano T., Hayashi C., 1970, Progress of Theoretical Physics, 44, 1580
 Lafrenière D., Jayawardhana R., van Kerkwijk M. H., 2010, ApJ, 719, 497
 Lagrange A.M. et al., 2009, A&A, 493, L21
 Laughlin G., Bodenheimer P., 1994, ApJ, 436, 335
 Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630
 —, 2005, MNRAS, 358, 1489
 Machida M. N., Inutsuka S., Matsumoto T., 2010, ApJ, 724, 1006
 —, 2011, ApJ, 729, 42
 Machida M. N., Matsumoto T., 2011, MNRAS, 413, 2767
 Machida M. N., Tomisaka K., Matsumoto T., Inutsuka S., 2008, ApJ, 677, 327
 Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J., Lafrenière D., Doyon R., 2008, Science, 322, 1348
 Marois C., Zuckerman B., Konopacky Q. M., Macintosh B., Barman T., 2010, Nature, 468, 1080
 Masunaga H., Inutsuka S., 1999, ApJ, 510, 822
 —, 2000, ApJ, 531, 350
 Matsumoto T., Hanawa T., 2003, ApJ, 595, 913
 Meru F., Bate M. R., 2012, MNRAS, 427, 2022
 Nayakshin S., 2010, MNRAS, 408, 2381
 Nelson A. F., 2006, MNRAS, 373, 1039
 Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627
 Paardekooper S.J., Baruteau C., Meru F., 2011, MNRAS, 416, L65
 Pascucci I., Wolf S., Steinacker J., Dullemond C. P., Henning T., Niccolini G., Woitke P., Lopez B., 2004, A&A, 417, 793
 Price D. J., 2007, PASA, 24, 159
 Rafikov R. R., 2005, ApJ, 621, L69
 Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS, 339, 1025
 Rice W. K. M., Forgan D. H., Armitage P. J., 2012, MNRAS, 420, 1640
 Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 364, L56
 Rogers P. D., Wadsley J., 2011, MNRAS, 414, 913
 —, 2012, MNRAS, 423, 1896
 Saigo K., Tomisaka K., 2006, ApJ, 645, 381
 Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
 Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
 Stamatellos D., Maury A., Whitworth A., André P., 2011, MNRAS, 413, 1787
 Stamatellos D., Whitworth A. P., 2008, A&A, 480, 879
 —, 2009a, MNRAS, 392, 413
 —, 2009b, MNRAS, 400, 1563
 Stamatellos D., Whitworth A. P., Hubber D. A., 2011, ApJ, 730, 32
 Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257
 Thalmann C. et al., 2009, ApJ, 707, L123
 Tomida K., Tomisaka K., Matsumoto T., Hori Y., Okuzumi S., Machida M. N., Saigo K., 2013, ApJ, 763, 6
 Toomre A., 1964, ApJ, 139, 1217
 Tsukamoto Y., Iwasaki K., Inutsuka S., 2013, MNRAS, 434, 2593
 Tsukamoto Y., Machida M. N., 2011, MNRAS, 416, 591
 —, 2013, MNRAS, 428, 1321
 Tsukamoto Y., Machida M. N., Inutsuka S., 2013, MNRAS, 436, 1667
 Vorobyov E. I., Basu S., 2006, ApJ, 650, 956
 —, 2010a, ApJ, 714, L133
 —, 2010b, ApJ, 719, 1896
 Vorobyov E. I., DeSouza A. L., Basu S., 2013, ApJ, 768, 131
 Whitehouse S. C., Bate M. R., 2004, MNRAS, 353, 1078
 Whitehouse S. C., Bate M. R., Monaghan J. J., 2005, MNRAS, 364, 1367
 Yorke H. W., Bodenheimer P., Laughlin G., 1993, ApJ, 411, 274