The Origin of Planetary System Architectures. I.

The Origin of Planetary System Architectures. I. Multiple Planet Traps in Gaseous Discs

Yasuhiro Hasegawa and Ralph E. Pudritz11footnotemark: 1
Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S 4M1, Canada
Origins Institute, McMaster University, Hamilton, ON L8S 4M1, Canada
E-mail: (YH); (REP)

The structure of planetary systems around their host stars depends on their initial formation conditions. Massive planets will likely be formed as a consequence of rapid migration of planetesimals and low mass cores into specific trapping sites in protoplanetary discs. We present analytical modeling of inhomogeneities in protoplanetary discs around a variety of young stars, - from Herbig Ae/Be to classical T Tauri and down to M stars, - and show how they give rise to planet traps. The positions of these traps define the initial orbital distribution of multiple protoplanets. We investigate both corotation and Lindblad torques, and show that a new trap arises from the (entropy-related) corotation torque. This arises at that disc radius where disc heating changes from viscous to stellar irradiation dominated processes. We demonstrate that up to three traps (heat transitions, ice lines and dead zones) can exist in a single disc, and that they move differently as the disc accretion rate decreases with time. The interaction between the giant planets which grow in such traps may be a crucial ingredient for establishing planetary systems. We also demonstrate that the position of planet traps strongly depends on stellar masses and disc accretion rates. This indicates that host stars establish preferred scales of planetary systems formed around them. We discuss the potential of planet traps induced by ice lines of various molecules such as water and CO, and estimate the maximum and minimum mass of planets which undergo type I migration. We finally apply our analyses to accounting for the initial conditions proposed in the Nice model for the origin of our Solar system.

accretion, accretion discs – turbulence – planets and satellites: formation – planet-disc interactions – protoplanetary discs – (stars:) planetary systems
pagerange: The Origin of Planetary System Architectures. I. Multiple Planet Traps in Gaseous DiscsCpubyear: 2009

1 Introduction

Observations of exoplanets (nearly 1700 if candidates are included) show that the formation of multiple planets around their host stars is relatively common.111See the website This trend is confirmed by both the radial velocity and transit techniques such as the Kepler mission (e.g. Howard et al., 2011). Many previous studies based on N-body simulations, investigated the dynamics of the planetary systems (e.g. Rasio & Ford, 1996). It is well known that these simulations can reproduce the observed distribution of eccentricities of exoplanets very well, if they adopt a specific initial condition that planets are closely packed (e.g. Ford & Rasio, 2008). For our Solar system, the Nice model which requires a specific initial arrangement of Jupiter and Saturn explains the dynamics very well (Morbidelli, 2010, references herein). The obvious question is whether or not such conditions are reasonable. In order to answer the question, one needs to consider the early stage of planet formation in which physical processes are controlled by gas dynamics in discs (e.g. Thommes et al., 2008). The presence of gas is necessary for gas giants to be formed. However, it also causes rapid inward planetary migration (Ward, 1997). Since this radial motion strongly depends on the disc properties such as the gas surface density and the disc temperature, it becomes a huge challenge to systematically investigate what are the realistic initial conditions for the later evolution of planetary systems.

Recently, planet traps have received a lot of attention (Masset et al., 2006; Morbidelli et al., 2008). Barriers to planetary migration arise because of inhomogeneities in discs where the direction of planetary migration switches from inwards to outwards, so that migrating planets are halted. (Equivalently, the net torque exerted on planets is zero there, Hasegawa & Pudritz, 2010a, hereafter HP10, references herein). Planets may acquire most of their mass as they accrete material at these barriers - which we call planet traps in globally evolving discs. We use the term barriers for our local analyses ( 3 and 4) while the term of planet traps are used for our global, unified analyses ( 6).

Planet traps were originally proposed by Masset et al. (2006) in order to solve the well known rapid migration problem wherein planets can be lost to discs within years. Matsumura et al. (2007, hereafter MPT07) first addressed a (possible) link between the planet traps and the diversity of exoplanets by showing that planet traps can move due to the time dependent, viscous evolution of discs (also see Matsumura et al., 2009). Although they focused on dead zones in discs (which are the high density, inner regions, so that turbulence there induced by magnetorotational instability (MRI) is quenched, Gammie, 1996), their basic idea is valid for any type of planet trap. More recently, it has been pointed out that a single disc can have a couple of planet traps (Ida & Lin, 2008; Lyra et al., 2010). Population synthesis models confirm that the planet traps and their movements can be important for the diversity of observed exoplanetary systems (Mordasini et al., 2011). As already noted, planet traps have the possibility of strongly enhancing the growth rate of planetary cores (Sándor et al., 2011) and the formation of giant planets. Thus, planet traps have considerable potential for understanding the formation of planetary systems.

In a series of papers, we systematically investigate how inhomogeneities in discs can trap migrating planets in overdense regions where they undergo most of their growth, and how trapped planets establish their planetary systems in viscously evolving discs. In this paper, we undertake a comprehensive study of the various mechanisms that produce planet traps in discs around a variety of young stars, - from high to low mass. One of our major findings is that a single disc can have up to three planet trap regions. We show that the position of any trap depends upon the disc’s accretion rate, which decreases with time as the disc is used up and star formation is terminated. The decreasing accretion rate forces the planet traps to move inwards at varying rates. This ultimately sets up the condition for the mutual interaction of planets in the traps, which can provide the realistic initial conditions for the evolution of planetary systems.

The plan of this paper is as follows. We describe physical processes that create inhomogeneities in discs and summarise our analytical approach for them in 2. The general reader may then proceed to 6 for discussion of the general results (also see Fig. 7). Armed with physical understanding of the inhomogeneities, we investigate corotation and Lindblad torques which are the driving force of planetary migration in 3 and 4, respectively. For the former case, we show that a new barrier arises due to the change of the main heat source for gas disc from viscosity to stellar irradiation. For both torques, ice lines play some role in trapping planets (Ida & Lin 2008, hereafter IL08; Lyra et al. 2010). In 5, we discuss other possible barriers and estimate the maximum and minimum mass with which planets undergo type I migration. Also, we investigate the possibility of the presence of ice lines due to various molecules. In 6, we integrate our analyses and discuss the roles of these planet traps in the formation of planetary systems. We apply our results to an explanation of the Nice model for the architecture of the Solar system in 7. In 8, we present our conclusions. We summarise important quantities which often appear in this paper in Table 1.

Symbols Meaning
Stellar mass
Stellar radius
Stellar effective temperature
Accretion rate (see equation (3.6))
Gravitational constant
Angular frequency (see equation (4.3))
Keplerian frequency ()
Planetary mass
Planetary orbital radius
Gas volume density
Gas surface density
Gas surface density at
Power-law index of if
Disc temperature ()
Power-law index of
Disc effective temperature
Disc temperature of the mid-plane
Condensation temperature for species at the ice line
Disc temperature of the surface
Sound speed ()
Disc scale height ()
Disc photosphere height
Disc aspect ratio ()
Disc radius of disc inhomogeneities
Gas density modification at
Transition width at ()
Disc radius of the heat transition
Disc radius of ice lines
Disc radius of the outer edge of dead zones
Gas pressure ()
Epicyclic frequency (see equation (4.4))
Wavenumber (see equation (4.2))
Lindblad resonant position ()
Forcing function (see equations (4.5) and (4.12))
Value of a quantity at
Surface density of active regions (see equation (4.24))
Power-law index of (generally )
Characteristic disc radius
Modification factor of due to ice lines ()
Parameter representing density jumps due to ice lines
Parameter representing dust traps due to ice lines
Parameter representing density bumps due to ice lines
Mean strength of turbulence (see equation (4.21))
Strength of turbulence in the active zone ()
Strength of turbulence in the dead zone ()
Kinematic viscosity ()
Adiabatic index (=1.4)
Heat flux in the vertical direction
Viscous dissipation rate per unit mass
the Stefan-Boltzmann constant
the Boltzmann constant
Grazing angle (see equation (3.20))
Optical depth ()

The simplest assumption that is only adopted in 3.2 and 4.2.

Table 1: Important quantities

2 Disc inhomogeneities

We describe physical processes governing the structure of protoplanetary discs and discuss how disc inhomogeneities arises from these processes. We discuss the basic features of our analytical modeling of the resultant disc structures - which will be presented in 3 and 4. Table 2 summarises the disc inhomogeneities, related nomenclature that often appears in the literature, the dominant torque that transports angular momentum in that region of the disc, and the section of the paper that treats the analysis.

Disc inhomogeneity Nomenclature Torque Section
Opacity transitions e.g. Ice lines Corotation 3.2.4
Heat transitions N/A Corotation 3.2.3
Turbulence transitions Dead zones Lindblad 4.3
Opacity and turbulent transitions Ice lines Lindbald 4.4

Disc turbulence is assumed to originate from non-MRI based processes

Disc turbulence is assumed to originate from MRI based processes

Table 2: Summary of disc inhomogeneities

2.1 Physical processes

Heating of protoplanetary discs is one of the most important processes in order to understand their geometrical structure (e.g. Dullemond et al., 2007). Since discs are accreted onto the central stars, release of gravitational energy through disc accretion becomes one of the main heat sources. This energy can be dissipated by viscous stresses, leading to viscous heating. Once discs are heated up, then the absorption efficiency of discs which is regulated by their optical depth establishes the thermal structure of discs. In protoplanetary discs, dust gives the main contribution to opacity for photons with low to intermediate energy while gas is the main absorber of high energy photons. In the inner region of discs, viscous heating is very efficient and leads to high disc temperatures ( 1000 K for classical T Tauri stars (CTTSs), D’Alessio et al., 1998). As a result, both metal dust grains and molecules are destroyed there (e.g. Bell & Lin, 1994). Thus, this high disc temperature and resultant destruction of opacity sources produces opacity transitions. It is interesting that opacity transitions are also produced in a completely different situation, that is, of low disc temperatures wherein opacity is enhanced by ”freeze-out” processes. Ice lines are one of the most famous opacity transitions created by this process. Ice lines are defined such that the number density of ice-coated dust grains suddenly increases, which is a consequence of low disc temperatures there. In 3.2.4, we treat ice lines as an example of opacity transitions.

It is well known that protoplanetary discs are also heated up by stellar irradiation (Chiang & Goldreich, 1997; Hasegawa & Pudritz, 2010b). Since viscous heating becomes less efficient for larger disc radii, stellar irradiation plays a dominant role in regulating the thermal structure of discs there. Thus, a transition exists for discs wherein the dominant heating mechanism transits from viscous heating to stellar irradiation. We call this a heat transition. (see 3.2.3).

We have focused so far on heating processes which involve with photons that have relatively low energy. High energy photons such as X-rays from the central stars and cosmic rays, are important for understanding the ionization of protoplanetary discs. Since the MRI is the most favoured process to excite turbulence in discs, it is crucial to evaluate the ionization structure of discs (Gammie, 1996). In the inner region of discs, the column density is very high, so that high energy photons cannot penetrate the entire region. In this ”layered” structure, the mid-plane is effectively screened from being ionized while the surface region is efficiently ionized. Thus, the inner region is less turbulent, and characterised by a dead zone. On the other hand, the outer region has lower column density and hence is fully ionized for the entire region including the mid-plane. As a result, the outer region is fully turbulent. Thus, transitions in the amplitude of turbulence exist in discs which we denote as turbulent transitions (see 4.3). Dead zones are the most famous product of the turbulent transitions.

Finally, we introduce interesting transitions which arise from the combination of opacity and turbulent transitions. If the MRI is the mechanism that drives disc turbulence, then ice lines can also be regarded as this type of transition (IL08, see 4.4.1 for the complete discussion). At the ice lines, the sudden increment of ice-coated dust grains is likely to result in the sudden decrement of free electrons there, since such sticky grains can efficiently absorb them (Sano et al., 2000). This process, initiated by the opacity transitions, reduces the ionization level there, and therefore removes the coupling between the magnetic field and the gas - which kills the MRI instability. Hence ice lines create turbulent transitions. In summary, ice lines are regarded as an opacity transition if disc turbulence is excited by non MRI processes while they act as an opacity and turbulent transition if turbulence is excited by the MRI. In 4.4, we investigate ice lines again, but as an example of the opacity and turbulent transitions.

2.2 Our analytical approach

We present analytical modeling of the disc structures affected by the disc inhomogeneities in 3 and 4. In order to make analytical treatments possible, we adopt expressions that are simple enough, but well capture the physics arising from the inhomogeneities. As an example, we adopt a function for representing an inhomogeneity in the surface density (see equation (3.1)). This profile is more general than a power-law and is likely to be applicable for the cases of opacity, heat, and turbulent transitions. For the opacity and heat transitions, the results of Menou & Goodman (2004, hereafter MG04, see their fig1) validate the usage of this function as do the results of MPT07 for the turbulent transitions. Disc structures are approximated as power-laws for regions far away from the inhomogeneities (see equation (3.1)). Thus, we adopt more general profiles for the disc structures. In addition, we make use of the simplest expressions in order to reduce mathematical complexity, and hence different functions are adopted for different transitions.

Bearing these in mind, we will demonstrate that the disc inhomogeneities are the most plausible sites to produce barriers to type I migration by undertaking a comprehensive study of both corotation and Lindblad torques in 3 and 4, respectively.

3 Barriers arising from corotation torques

Corotation torques play an important role in planetary migration in regions with high viscosity (), because the corotation torque there is (partially) unsaturated (i.e. effective, see Appendix B in Hasegawa & Pudritz 2011, references herein). The corotation torque is known to act as a barrier for two situations. The first situation is when with (Masset et al., 2006; Paardekooper & Papaloizou, 2009). This case can be generally established only at the inner edge of discs. In typical discs, the inner edge is located around a few hundredths of au, which corresponds to the semi-major axis of Hot Jupiters. Thus, the positive gradients of surface density unlikely play a dominant role in explaining the diversity of the observed planetary systems. The second situation arises for adiabatic discs that have inhomogeneities (Lyra et al., 2010). This is more promising since the main sites of planetary formation in protoplanetary discs are generally optically thick and reasonably approximated to be adiabatic. In this section, we investigate only the second case for the above reason.

3.1 Conditions for outward migration

We derive conditions for the disc structures that are required for planets to migrate outwards. These conditions are crucial for the subsequent subsection where we discuss how disc inhomogeneities work as a barrier.

3.1.1 Disc models

We discuss our disc models (also see Table 1). For the surface density affected by disc inhomogeneities, we adopt the following equation;


where is the initial density profile, characterise the density distortion which is a consequence of disc inhomogeneities, is a orbital radius of the disc inhomogeneities, and is the width of the transition (also see Table 1). This analytical modeling is motivated by the results of MG04 who first found the significant effects of opacity transitions on planetary migration by solving the detailed, 1D disc structure equations. Based on their results (see their fig. 1), any temperature distortion created by opacity or heat transitions, is likely to result in a surface density structure which is well expressed by equation (3.1). For the initial profile (), we examine two cases both of which are well discussed in the literature: and . The former power-law index is known as minimum mass solar nebula (MMSN) while the latter one is a steady state solution to disc accretion. For the disc temperature, we adopt simple, power-law structures () and derive the critical value of that results in outward migration below.

3.1.2 Torque formula

We adopt the torque formula derived by Paardekooper et al. (2010) in which discs are assumed to be 2D. In the formula, the total torque is comprised of the linear Lindblad torque and non-linear corotation torques, known as horseshoe drags. If discs are (locally) isothermal, both Lindblad and corotation torques dictate that standard rapid inward migration will occur. More specifically, the vortensity-related horseshoe drag that is only an active (corotation) torque in this case, cannot exceed the Lindblad torque for power-law discs (Paardekooper & Papaloizou, 2009). For discs with more general profiles, the vortensity-related horseshoe drag can exceed the Lindblad torque, but still results in inward migration. Therefore, we assume discs to be adiabatic, which is reasonable for the inner part of discs ( 80 au) around CTTSs (Chiang & Goldreich, 1997). Adopting the torque formula in adiabatic discs (see Paardekooper et al., 2010), the direction of migration may be given as


where is an ”effective” power-law index for which can have a more general profile, is the adiabatic index and we took the softening length for simplicity (also see Table 1). The terms in the first brackets arise from the Lindblad torque, the terms in the second from the vortensity-related horseshoe drag, and the terms in the third from the entropy-related horseshoe drag. As Paardekooper & Mellema (2006) first discovered and subsequent works interpreted it, the entropy-related horseshoe drag which occurs only in adiabatic discs, is scaled by radial, entropy gradients and can drive planets into outward migration (see Appendix B in Hasegawa & Pudritz 2011 for a summary, references herein).

The original formula derived by Paardekooper et al. (2010) is valid exclusively in power-law discs. In order to take into account discs with more general profiles such as equation (3.1), we add a vortensity correction factor in the second term which is defined as




is the vorticity of gas, also known as one of the Oort constants. This is because gradients of vortensities () are very sensitive to disc structure.222Lindblad torques also need a similar correction factor for discs with general profiles. Indeed, we derive an analytical relation for Lindblad torques using equation (3.1) in 4.3 (see equation (4.20)). As shown in Appendix A, however, they switch the direction of migration only if for . In general, the density distortion created by opacity and heat transitions is likely to be less than that. In addition, it is more consistent to use the above terms than equation (4.20) in this torque formulation. Therefore, we adopt the original form for the Lindblad torques. We note that for pure power-law behaviours. In fact, gradients of vortensities are the core of the vortensity-related corotation torque and regulate the transfer of angular momentum there (see Appendix B in Hasegawa & Pudritz 2011). Therefore, inclusion of the factor is likely to be crucial for properly evaluating the importance of the vortensity-related corotation torque relative to the others.

The condition for outward migration is, therefore,


It is fruitful to first consider pure power-law discs where . Setting , the required temperature profile is for MMSN discs () while for , the temperature profile is needed. We stress that such steep temperature profiles can be only achieved by viscous heating in optically thick discs, which is discussed more in the next subsection.

Let us now investigate how the required deviates from the power-law predictions due to the factor . In general, we find that the structure of is very complicated. Therefore, we present the detail discussion of in Appendix A and briefly summarise the two most important effects on here. The first is that the effects of are very local and are likely to be confined within the transition region . This is clearly shown in Fig. 1 (see the vertical dotted line on the bottom panel for representing the transition region). In this figure, we set that is the most likely value for the opacity and heat transitions. For the bottom panel, we plot the critical value of derived from equation (3.5) as a function of , where the characteristic disc radius is set au unless otherwise stated. The top panel shows the behaviour of . For both panels, the solid line denotes the results for the disc with the general profiles (see equation (3.1)) while the dashed line is for the pure power-law discs. We show the results only for the case that , because the results of are qualitatively similar. The second important point is that, in the local region where the effects of become crucial, a required reaches values that are unattainable by any physical process in protoplanetary discs. This means that outward migration cannot happen due to large vortensity-related corotation torques that results in inward migration. In summary, the inclusion of reduces the possible region wherein planets can migrate outwards, but this reduction is well confined in a local region ( an order of ) centered at .

Finally, we note that this torque formula does not take into account any effect of saturation. The problem of saturation, which is central to the problem of corotation torque, is very complicated and strongly dependent on the flow pattern at the horseshoe orbit (Masset & Casoli 2010; also see Appendix B in Hasegawa & Pudritz 2011). Paardekooper et al. (2011); Masset & Casoli (2010) attempted to derive analytical formulae in which the effects of saturation are included. However, both formulae depend on thermal diffusivity in discs that is totally unknown in protoplanetary discs. Hence, we adopted the unsaturated torque formula. This implies that the above required temperature profiles may be the minimum value. Steeper profiles may be needed if (partial) saturation effects are taken into account.

Figure 1: The effects of the vortensity correction factor on and on the top and bottom panel, respectively. For both panels, the solid lines denote the results for discs with general profiles (see equation (3.1)) while the dashed lines are for pure power-law discs. Only the results of are shown because those of the MMSN case are qualitatively similar. We set the transition width parameter . For both panels, the effects of are well confined in the transition region () (see the vertical dotted lines on the bottom panel). In addition, inclusion of reduces the possibility of outward migration there.

3.2 Heat and opacity transition barriers

The entropy-related horseshoe drag can result in outward migration only in discs with steep temperature profiles, as discussed above. Only viscous heating can establish such profiles. On the other hand, stellar irradiation can heat up protoplanetary discs as well. By deriving temperature profiles for both heating processes, we show that the heating transition from viscosity to stellar irradiation activates a new barrier. In addition, we examine the effects of opacity transitions on the temperature structure. We exclusively focus on ice lines as an example of the opacity transitions. For both transitions, we identify the positions of barriers.

3.2.1 Disc models

We adopt simple, power-law profiles for both the surface density and disc temperature ( and ). This is supported by the argument done in the above subsection. Since the effects of occurs only in a local region centered at the position of a transition, none of our findings discussed below is affected. In addition, the resultant deviation for our estimate of the position of barriers is only 10 per cents. Thus, it is reasonable to use power-law discs here.

We take the value of surface density at (labeled by ) as 10 times larger than the standard values of discs for all cases. Many previous studies showed that gas giants can be formed only if the disc mass increases by this amount relative to the MMSN models (e.g. Alibert et al., 2005). Table 3 tabulates the values of around stars with various masses.

In addition, we assume stationary accretion discs that have accretion rates which can be modeled as


using the famous -prescription (Shakura & Sunyaev, 1973).

Herbig Ae/Be stars CTTSs M stars
(/ year)
() 2.5 0.5 0.1
() 2 2.5 0.4
(K) 10000 4000 2850
(g cm)

The value of is 10 times larger than the standard values.

Table 3: Typical quantities for stars with various masses

3.2.2 Viscous heating: Origin of outward migration

We first discuss viscous heating. Assuming local thermal equilibrium in geometrically thin discs, the energy equation can be written as


where is the heat flux in the vertical direction and is the viscous energy dissipation rate per unit mass (e.g. Ruden & Lin, 1986, also see Table 1). In Keplerian discs, the viscous dissipation rate becomes


Integrating equation (3.7) gives


If discs are assumed to be optically thin and isothermal, the flux with the Stefan-Boltzmann constant . Thus, the resultant effective temperature profile becomes


for discs which are assumed to reach a steady state. In this case which is valid in the late stage of disc evolution, accretion rates () is constant over the entire disc. This famous law was first derived by Lynden-Bell & Pringle (1974) and is exactly identical to the temperature profile for flat discs which are heated by stellar irradiation (Adams et al., 1987; Chiang & Goldreich, 1997). In addition, it is well known that the resultant spectral energy distributions (SEDs) do not reproduce the observed ones - flared discs are required.

If the accretion rate () changes with disc radius, then the self-consistent temperature profile is


where is the surface density of discs. For the MMSN models (), while for discs with , . It is obvious that the temperature profiles derived from the isothermal assumption are not steep enough for the corotation torque to provide a barrier.

The isothermal assumption can break down, especially in the main site of planetary formation . The region of discs is reasonably considered to be optically thick and the energy is mainly transported by radiation (D’Alessio et al., 1998). In this case, the flux is generally described by


where is the opacity (e.g. Ruden & Lin, 1986, also see Table 1). Assuming the flux is constant over , integration of equation (3.12) gives the relation between the temperatures of surface and midplane regions;


where is the optical depth and and is the temperature of surface and midplane region, respectively (also see Nakamoto & Nakagawa, 1994, Table 1). If only viscous heating is taken into account, the temperature of the surface region should be much smaller than that of the midplane. As a result, the temperature structure is governed by


When radiative transfer equation is treated explicitly, is replaced by (Hubeny, 1990; Kley & Crida, 2008);


If discs have a constant (see equation (3.6)), then the temperature profile becomes


assuming to be independent of and . For the MMSN models (), while for discs with , . Again, these temperature profiles are not steep enough.

However, Kley & Crida (2008) performed numerical simulations by solving a more complicated energy equation, and showed that the resultant temperature profile goes to for discs with in steady state. We can gain a similar profile if we adopt (Bell & Lin, 1994). In this case, . Therefore, the assumption that is independent of and , likely underestimates the temperature structure. If is not constant, then the self-consistent temperature profile becomes


assuming to be independent of and . For the MMSN models (), while for discs with , . These profiles are exactly what is required (equation (3.2)) in order for planets to migrate outwards. If is adopted, , which is more preferred for outward migration.

Thus, outward migration due to the entropy-related horseshoe drag is expected for viscously heated, optically thick discs. If protoplanetary discs were homogeneous in opacity and viscosity were the only physical process heating them, the results would suggest that large region of discs ought to be devoid of planets. This does not occur, however. We examine two kinds of inhomogeneity of discs below. One of them is stellar irradiation, which gives a new barrier. The other is the ice line which arises as a consequence of an opacity transition (Lyra et al., 2010).

3.2.3 Stellar irradiation: Origin of heat transition barriers

Stellar irradiation is well known to be the main heat source for regions in the disc beyond au (in CTTSs, D’Alessio et al., 1998). Inside of that radius, viscous heating dominates over stellar irradiation. Since the temperature slope controlled by stellar irradiation is much shallower than that of viscous heating, planets which migrate inwards can be halted at the turning point where viscous heating begins to take over. We call this barrier a heat transition barrier.

In order to calculate the position of the turning point, we adopt radiative disc models of Chiang & Goldreich (1997). In this model, two kinds of disc temperature are calculated. One of them represents the surface layer which is directly heated by the central star, called the super-heated layer. The other is for the midplane layer which is heated by the super-heated layer. Assuming radiative equilibrium discs, the temperature of the super-heated layer is given as


where K for discs around the CTTSs. Since the super-heated layer radiates equal amounts of energy inwards and outwards, the temperature of the midplane layer which is approximated to be optically thick is written as


where is the grazing angle at which photons emitted from the star strike the discs (also see Table 1). In general, is written as


where is the height of the visible photosphere above the mid-plane. Thus, the temperature of the midplane is strongly dependent on the geometrical shape of discs if stellar irradiation is taken into account.

For flat discs in which the aspect ratio is the constant, the temperature of the midplane layer becomes


Again, this is the famous law. For flared discs which are needed for reproducing the observed SEDs (Kenyon & Hartmann, 1987), the second term of the right hand side in equation (3.20) becomes dominant. Assuming that dust is well mixed with the gas in discs,




is the mean molecular weight of the gas and is the Boltzmann constant (also see Table 1). Assuming is constant, the self-consistent temperature of the midplane is




It is obvious that this temperature profile is much shallower than that due to viscous heating.

Thus, the radius at which temperature profile switches is (equating equations (3.14) and (3.24))


where we adopt with (Bell & Lin, 1994). The usage of this form of that is only valid for the region outside of ice lines is reasonable, which is discussed below.

Fig. 2 shows the heat transition radius for discs around stars with various masses. For the solid lines, we adopt the values of in Table 3 with (also see Table 1). For comparison purposes, the dashed lines denote the case with (which is the standard surface density in the literature). The transition radius generally increases with increasing . Otherwise, it becomes a decreasing function of (see the dashed line for the case of M stars). Table 4 summarises the values of for both and with . It is obvious that massive discs around Herbig Ae/Be stars extend to the order of several hundred au while lowest mass discs around M stars shrink to the order of 1 au. For the intermediate disc mass around CTTSs, is the order of a few ten au.

Figure 2: The heat transition radius as a function of the power-law index . For all panels, the solid line denotes the case using the values in Table 3 while only decreases by a factor of 10 for the dashed line, which is the value used more often in the literature. In general, the radius increase with increasing , since the quantity in the brackets of equation (3.26) is larger than unity. If it is less than unity (see the dashed line for the case of M stars), the radius decreases.
Herbig Ae/Be stars CTTSs M stars
(au) for 159 16 1.8
(au) for 827 40 2.3
Table 4: The typical values of for

3.2.4 Ice lines: Origin of ice line barriers

We now discuss the effects of ice lines on temperature profiles which provide another barrier (Lyra et al., 2010). We note that our following approach is applicable to any other opacity transitions. We discuss them more in 5.2. At the ice line, the number density of icy dust grains suddenly increases due to the low disc temperature. This sudden increment of dust provides an opacity transition. As a result, the temperature profile around the ice lines becomes shallow, resulting in inward migration while the profile well inside of the ice line becomes steep enough for planets to migrate outwards. Thus, planets migrating inwards are halted around the ice lines. This is known as the ice line barrier. This assumes that corotation torques are active -i.e. non-saturated, which is not necessarily clear at the ice lines (see 4.4.1 and 5.2).

At first, we examine which heat source controls the location of ice lines, . In order to proceed, we adopt equation (3.14) for viscous heating and equation (3.24) for stellar irradiation. Equating to the condensation temperature for HO ice, K (Jang-Condell & Sasselov, 2004), equation (3.14) becomes


where with (Bell & Lin, 1994), and equation (3.24) becomes


In general, the exponent index is negative, so that . This shows that as the surface density () decreases, the position of the ice line barriers moves inwards. We note that we examine the effects of a water-ice line here although the above argument is applicable to ice lines of any material by changing the relevant condensation temperature for species , . We will discuss them more in 5.3.

Fig. 3 shows the above two equations as a function of for discs around stars with various masses. We set , since this choice of minimises the importance of viscous heating in discs (see Fig. 2). The black, solid and dashed lines denote equation (3.27) with and , respectively while the dotted line is for equation (3.28). In all cases, defined by viscous heating is located at greater distances from the star than that by stellar irradiation. Furthermore, the transition radius defined by equation (3.26) is larger than any of these two radii (see the gray, thick lines). Thus, we can conclude that is determined by viscous heating for discs with a wide range of . This agrees with numerical work by Min et al. (2011) who showed the same results by numerically solving the full wavelength dependent, radiative transfer equation by means of a Monte Carlo method in 3D discs. In their simulations, viscous heating is explicitly included as well as stellar irradiation. In addition, this supports the findings of Lyra et al. (2010) who adopt disc models which are only heated by viscous heating. Furthermore, Fig. 3 shows that the radius of the heat transition is located outside of the ice line. This validates our choice of in equation (3.26).

We now investigate temperature profiles created by viscous heating around ice lines by adopting appropriate for each region (Bell & Lin, 1994). For the region inside of the ice lines, with . As a result, equation (3.14) becomes


for discs with general (see equation (3.6)). for the MMSN disc models () while for the disc models with . These steep profiles obviously indicate that planets migrate outwards. For the region around the ice lines, with . Consequently, equation (3.14) becomes


for discs with . for the MMSN disc models () while for the disc models with . Under these shallow profiles, the entropy-related horseshoe drag cannot exceed the Linblad torque, and hence planets migrate inwards. Thus, a planet that migrates inwards beyond the ice line will be halted near the ice line due to this opacity transition.

Figure 3: The location of the water-ice line as a function of the turbulence parameter, . For every panel, the solid and dashed, black lines denote the radius derived from viscous heating using and , respectively (see equation (3.27) and Table 3). The dotted lines are from stellar irradiation (see equation (3.28)). For comparison purposes, the heat transition radius is denoted by the two gray lines (see equation 3.26). For all panels, viscous heating defines the location of the water-ice lines. In addition, the heat transition radius is larger than that of the water-ice lines.

3.3 Comparison

We compare the heat transition barriers to the ice line barriers. We especially focus on the relative location of each barrier. As discussed above, Fig. 3 shows the relative location of these two barriers (see the black and gray lines). The ice line barrier (equation (3.27)) is located inside of the heat transition barrier (equation (3.26)) for discs with a wide range of . However, we again emphasise that the above analyses for these two barriers can be valid only for the disc region with (). In the region with a low value of such as dead zones (), any corotation torque will be saturated (i.e. zero), so that these two barriers are never activated there. We also note that for such low values of , ice lines can become barriers as a consequence of the Lindblad torque rather than the corotation torque. In addition, if the adiabatic approximation breaks down, which can arise for the late stage of disc evolution or the outer part of discs, these two mechanisms cannot be effective. We now examine the mechanisms of barriers activated by Lindblad torque.

4 Barriers arising from Linblad torques

Let us now examine regions wherein any torques arising from corotation resonances are negligible, so that the direction of migration is determined only by the Lindblad torque. Such regions occur for discs with a shallow temperature profile or a low value of in dead zones (, HP10; Paardekooper et al., 2011; Masset & Casoli, 2010; Hasegawa & Pudritz, 2011).

4.1 Lindblad torque

4.1.1 Basic equations

We adopt analytical formulae of the Lindblad torques described in Hasegawa & Pudritz (2011, see their Appendix A for the complete expression). In the formulae, the standard Lindblad torque derived by Ward (1997) is modified, so that the vertical thickness of discs is taken into account without rigorously solving the 3D Euler equations. Thus, we treat the Lindblad torque density of a layer at the height from the mid-plane which is written as


where for the outer(inner) resonances, the wavenumber is treated as a continuous variable;


is the angular velocity of a planet with , the angular velocity is


the epicyclic frequency is


, and the forcing function is




is the resonant position, and (also see Table 1). We have adopted the standard approximation for (Goldreich & Tremaine, 1980; Jang-Condell & Sasselov, 2005). In addition, we have added a factor in the second term of following MG04.

Finally, the total torque exerted by the planets on discs is found by summing over all layers;


Note that the integration of equation (4.1) in terms of with leads to the famous analytical formulae of the Lindblad torque density in 2D discs (Ward, 1997). Also, note that the tidal torques shown above are exerted by the planets on discs. Thus, the positive sign of in equation (4.1) means that the planets exert the tidal torques on discs, resulting in the loss of angular momentum from the planets and vice versa.

The validity of the above approach for mimicking the reduction of the tidal torques in 3D discs (Tanaka et al., 2002) is confirmed by Lubow & Ogilvie (1998) who showed that, in thermally stratified discs, 2D modes carry off more than 95 per cent of angular momentum and vertical modes are not dominant for transferring angular momentum between planets and their disc (also see Ward, 1988; Artymowicz, 1993; Jang-Condell & Sasselov, 2005).

4.1.2 Basic assumptions

We assume discs to be geometrically thin () and Keplerian. Under this assumption, we have


Thus, the torque density becomes






Thus, the behavior of the torque density depends on . In the following subsections, we expand every quantity ( and ) as a series in , by assuming simple power-law disc structures in 4.2 and non-power-law disc structures in 4.3 and 4.4. This expansion is assured, since the torque takes the maximum value at (Ward, 1997), and consequently, allows us to derive simple relations governing the direction of migration. Thus, we perform local analyses as done in 2.

4.2 Simple power-law structures

We derive an analytical relation for the direction of planetary migration. Here, we assume simple power-law structures for the background density and temperature in discs, that is, and . As a result, the relation depends only on the exponent of and .

We can further simplify the mathematics by taking the limit (since the effect of gas pressure becomes negligible in the Keplerian discs; see Appendix B). Equation (4.11) then becomes, to first order in ,




Thus, the total torque is


where the integration range for is approximated to extend over the height of the photosphere at the location of the planet since the density above the photosphere is about two orders of magnitude less than that at the mid-plane.

Consequently, the sign of the net torque depends on the sign expression;


The torque becomes positive (inward migration) when and negative (outward migration) when . When the MMSN disc models () are adopted, is needed for outward migration while, when disc models have , is required. In general, is a decreasing function of . Therefore, some specific features in discs are required for the Lindblad torque to excite barriers. In the following subsections, we focus on dead zones and ice lines.

4.3 Dead zone barriers

Dead zones are known to play very interesting roles in planetary formation and migration (e.g. Matsumura & Pudritz, 2006). They can excite two barriers: density barriers (MPT07) and thermal barriers (HP10). These two barriers are interesting products of turbulent transitions. In this subsection, we derive the analytical relation which predicts when the direction of migration is reversed from inwards to outwards for these cases.

4.3.1 Thermal barriers

Thermal barriers were uncovered by HP10 who first investigated numerically the effects of dead zones on the temperature structure of discs (also see Hasegawa & Pudritz, 2010b). In their paper, the temperature structure of discs with dust settling and a dead zone was simulated by solving the full wavelength dependent, radiative transfer equation by means of a Monte Carlo method. Stellar irradiation is assumed to be the main heat source. HP10 demonstrated that planets migrate outward in a region where a positive temperature gradient is established (see their figs 1 and 3). This temperature gradient is a result of the back heating of the dead zone by a thermally hot dusty wall, and is well represented by with . The hot dusty wall is produced by the enhanced dust settling in the dead zone and the resultant enhanced absorption of stellar irradiation by the wall. Since they adopted with , the profile of the positive temperature gradient is identical to the profile required by equation (4.17) in order for planets to migrate outwards. Thus, our relation predicts the results of the detail numerical simulations very well.

4.3.2 Density jumps

Density barriers were found by MPT07 who undertook a pioneering study on the effects of dead zones on planetary migration. They are created by the formation of a density jump at the outer edge of dead zones which arises as a consequence of time-dependent, viscous evolution of discs (also see Matsumura et al., 2009, for a more complete study). This density jump develops because of the huge difference in the strength of turbulence () between the active and dead zones of the disc.

Here, we analytically derive an analogous relation for density barriers. We adopt equation (3.1) as the background surface density. In this case, is the density difference between the active and dead zones, is a orbital radius of the outer boundary of the dead zone, and is the width of the transition (also see Table 1). Except for this, we adopt the same approximation as above. As a result, the relation depends on the structure of the dead zones as well as the exponent of and . Since planets are slowed down or stopped near , we can approximate .

We note that a slightly more accurate shape of the boundary in equation (3.1) can be expressed in terms of error functions. However, the difference between the error and functions is less than 1 per cent with the proper choice of coefficients (Basset, 1998). Furthermore, functions are fully analytical. For these reasons, we adopt functions.

Adopting the same approximation as above, we find


Note that we expanded in terms of and finally expressed it in terms of assuming


This is reasonable because the resonant positions are pushed away from the planet by a distance due to the gas pressure (Artymowicz, 1993; Ward, 1997, also see equation (4.10)).

Inserting the above surface density into equation (4.11), the sign of the net torque consequently becomes


Again, we can check the validity of equation (4.20) by comparing it with the numerical simulations. In this case, we compare the results of MPT07. Changing parameters and in equation (3.1), they investigated what values of and result in outward migration (see their Appendix B). Table 5 summarises their experiments with our predictions (given in the brackets). One observes immediately that equation (4.20) well explains the results of their numerical simulations. One might wonder if the corotation torque is significant around the density jump, but this is not the case. MPT07 clarified that planets migrate outwards due to the reversal of the balance of Lindblad torques, not corotation torques. This arises because the Lindblad resonant positions are located further away from planets than the corotation ones. Therefore, migrating planets first encounter the inner Lindblad torque at the density jump and this is sufficient to reverse the direction of migration.

In summary, both our simple analytical relations (equations (4.17) and (4.20)) well reproduce the results of the detailed numerical simulations of dead zones as planet traps. This indicates that our assumptions and treatments are reasonable to capture the physics arising in the more complicated numerical simulations.

=5 =10 =100
in (0.07) in (-0.28) out (-6.6)
out (-0.8) out (-2.2) out (-28)
out (-2) out (-4.9) out (-57)

=0.35, =-3/2, and =-1/2. in(out) means the inward(outward) migration calculated by MPT07. The number in the brackets is calculated from equation (4.20).

Table 5: Summary of parameter study

4.4 Ice line barriers

We finally examine ice line barriers. As discussed in 3.2.4, the disc radius of the ice lines is determined by viscous heating for a wide range of (see Fig. 3). In this subsection, we discuss how ice lines establish layered structures, and derive analytical relations for the resultant barrier. Thus, we investigate ice lines as an example of opacity and turbulent transitions. Again, we specify a ice line due to water, although all our analyses are applicable to ice lines of any material. The discussion is presented in 5.3.

4.4.1 Layered structures at ice lines

As mentioned before, the number density of icy grains suddenly increases at the ice lines. If discs are considered to be turbulent due to the MRI, this sudden increment of dust can result in layered structures: the MRI active, surface layer and the MRI dead, inner layer (Gammie, 1996). This can be understood as follows. At the ice lines, the number density of free electrons in the disc also suddenly drops because they are absorbed by such dust grains. Since free electrons are the main contributor of coupling with magnetic fields threading the disc (Sano et al., 2000), the surface density of the MRI active region is strongly diminished, and consequently the surface density of the MRI dead region is enhanced. As a result, a density bump appears which acts as a barrier (see Fig. 4). Thus, the mean value of around the ice lines becomes (see equation (4.21)), and hence, the ice lines can be regarded as a self-regulated, localised dead zone. In summary, the subsequent process initiated by icy grains can reduce the value of at the ice lines. The disc radius at which it appears is determined by viscous heating (see equation (3.27)).

We recall from 3.2.4 that ice lines become a barrier due to the corotation torque (also see IL08). However, any corotation torque may saturate at the ice lines because the turbulence dies out there - a conclusion which is valid for any disc in which turbulence is excited by the MRI. If the MRI is the most general source of disc turbulence, we expect that the ice lines may be natural barriers arising from the Lindblad torques rather than corotation torques.

4.4.2 Disc models

We adopt a parameterised treatment of dead zones in a disc described by IL08. In the region with layered structures, the effective can be generally written as


where is the surface density of the active layer, and and are the strength of turbulence in the active and dead layers, respectively (Kretke & Lin, 2007; Matsumura et al., 2009, also see Table 1). Assuming stationary accretion discs (see equation (3.6)) with equation (4.21), the surface density is given as


where cannot exceed . Thus, equation (4.22) is useful for representing the surface density within dead zones where . Also, this equation implies that the structure of in the dead zones strongly depends on . However, the structure of is not well constrained by theory, since it is very sensitive to the distribution of dust grains as well as the chemical models (Sano et al., 2000; Ilgner & Nelson, 2006). Therefore, we adopt simple prescriptions proposed by Kretke & Lin (2007) and IL08 as follows. For the active region where , the surface density simply becomes


if there is no ice line. The effects of the ice lines on are described below.

4.4.3 Density jumps without ice lines

At first, we examine density barriers produced by dead zones (see equation (4.22)) to compare the our previous analysis done in 4.3. For the surface density of the active layer, we assume


following Kretke & Lin (2007). We note that and are very sensitive to the dust distribution (also see Table 1). However, we set constant and here. A more detailed analysis is presented below where the effects of ice lines on are included. In discs with layered structures, it is generally considered that .

Using the above equation, equation (4.22) is approximately written as

where we performed Taylor expansion of in terms of , and finally expressed it by , assuming . Thus, the direction of migration is determined by


Since the outer edge of dead zones where a density jump develops is established around the region with , we can assume


We emphasise that equation (4.27) gives the position of planets halted by the density barrier which defines the outer edge of the dead zones (). In addition, this equation is another expression that the necessary condition () is safely satisfied. As a result, the sign of the net torque is written as


Compared with equation (4.17), the layered structure with a density jump gives an additional term which controls the direction of migration (since ). Thus, planets migrate inwards when sgn while they migrate outwards when . In standard disc models, and . Our relation therefore dictates outward migration, which is also consistent with our previous analysis done in 4.3.2.

4.4.4 Density jumps with ice lines

We incorporate the effects of ice lines on the surface density of the active layer, following IL08 (see Fig. 4). In this case, we need to consider two separate cases, depending on the ratio of to , where is the disc radius of an ice line. For the case that , we can use equations (4.22) and (4.24) as the surface density while we need to modify equation (4.23) for the case that .

We first examine the case that , that is, an ice line is located inside of a dead zone. In this case, in equation (4.24) is generally given as


where and parameterise the effects of an ice line at which the surface density of the active layer suddenly drops. As noted by Kretke & Lin (2007), the surface density distortion induced by ice lines produces a radial, positive pressure gradient there, and consequently dust can be trapped there. Otherwise, it migrates inwards due to the so-called head winds arising from the gas motion (Weidenschilling, 1977). Since the maximum migration speed of dust is times faster than the viscous evolution of gas, dust is quickly piled up at the ice lines. This effect can be included in which is generally written as


We note that inclusion of the dust effects may upgrade the ice line barriers to planet traps although we keep calling them barriers.

Fig. 4 shows the typical (as controlled by equations (4.29) and (4.30)) and the role of the parameters and in .

Consequently, we can approximate the surface density (see equations (4.22), (4.24), (4.29), and (4.30)) as