text:
# Low–Mass Proto–Planet Migration in T-Tauri 𝛼–Disks

## 1 Introduction

Since the first discovery of a planet orbiting a nearby Sun-like star about a decade ago (Mayor & Queloz 1995), radial velocity surveys of a few thousands of our closest stellar neighbors have uncovered more than a hundred such planets (Marcy & Butler 1998; Cumming, Marcy & Butler 1999; Marcy, Cochran & Mayor 2000). Transit searches, which have confirmed the gaseous giant nature of these planetary companions (Charbonneau et al. 2000), promise to uncover a large number of additional extrasolar planets in the future (Udalski et al. 2002; Konacki et al. 2003).

In addition to the large eccentricities of many of the planets discovered to date, one of the most striking features of this new population is the existence of a subset of extrasolar planets orbiting very close to their parent star (see, e.g., the extrasolar planet almanac111http://exoplanets.org/almanacframe.html and encyclopedia222http://www.obspm.fr/encycl/encycl.html for a census). It is generally accepted that these close-in extrasolar giant planets could not have formed in situ, but instead must have migrated from large distances, where conditions for planet formation are more favorable (see, e.g., Lin, Bodenheimer & Richardson 1996).

There are two physically distinct scenarios considered for the formation of gaseous giant planets (see, e.g., Ruden 1999; Wuchterl, Guillot & Lissauer 2000 for reviews). In the standard core accretion scenario, dust grain accumulation in a proto-planetary disk proceeds to the formation of planetesimals, which later assemble to form proto-planetary cores. Above a mass threshold of a few tens $`M_{}`$, runaway accretion of the surrounding gas onto the proto-planetary cores occurs, leading ultimately to the formation of giant planets made of gas for the most part (Safronov 1969; Mizuno 1980; Pollack et al. 1996; see also Rafikov 2003). Cores of giant planets form more easily beyond the ”snow line,” whose location is typically at AU distances from the parent star (Hayashi 1991; Boss 1996; Sasselov & Lecar 2000). In the alternative disk instability scenario, a rather massive proto-planetary disk is subject to gravitational instabilities leading to the direct formation of (possibly multiple) Jupiter-sized objects (Cameron 1978; Boss 1997; 2000, Armitage & Hansen 1999; Mayer et al. 2002; Lufkin et al. 2003). In the present study, our interests are more focused on the core accretion scenario because the evolution of embedded proto-planets is of general interest, whether or not gravitational instabilities are important for giant planet formation.

In both planet formation scenarios, inward migration is required to explain the existence of the population of close-in extrasolar giant planets. Here again, two physically distinct scenarios for migration have been put forward. Spiral density waves launched by a proto-planet embedded in a gaseous disk apply a negative torque on the proto-planet, which migrates inward as a result of orbital angular momentum losses (Goldreich & Tremaine 1980, henceforth GT80; Lin & Papaloizou 1986; Korycansky & Pollack 1993, henceforth KP93; Artymowicz 1993a,b; Ward 1997a). On the other hand, later in the system’s evolution, it is also possible for a planet embedded in a sufficiently massive planetesimal disk to experience inward migration via repeated planetesimal scattering events (Murray et al. 1998). Our interest in the present study is in the gaseous migration mechanism, partly because it surely affects the precursors of gas giants, at least within the core-instability scenario, whether subsequent planetesimal scattering migration occurs or not.

It has long been recognized that the standard core accretion and gaseous migration scenarios face a serious difficulty when combined. While estimated migration times for massive proto-planetary cores located at $``$ AU distances from their parent star are $`10^5`$ years (see, e.g., Ward 1997a; KP93), estimates for the time required to accrete a gaseous envelope are typically in the range $`10^6`$$`10^7`$ years (see, e.g., Pollack et al. 1996). Short inward migration times imply that cores would be accreted by the central star before they could build up any substantial gaseous envelope, and more generally, they pose a serious threat to the survival of any planetary system in formation (see, e.g., Ward 1997b).

Robust calculations for proto-planetary core migration times are notoriously difficult to achieve, however. While a forming planet with a mass $`>\mathrm{\hspace{0.33em}10}`$$`100M_{}`$ may open a gap and see its subsequent migration tied to the slow viscous evolution of the gaseous disk (“type–II” migration), a lower-mass proto-planetary core migrates in, relative to the gaseous component (“type–I” migration), at a rate which strongly depends on the exact disk conditions in the vicinity of the proto-planet. In that respect, most calculations of “type–I” migration rates to date were made with rather idealized models for the background disk. Both linear torque calculations (Ward 1997a; KP93) and fully non–linear hydrodynamical simulations (Kley, D’Angelo & Henning 2001; Nelson & Benz 2003a,b; D’Angelo, Kley & Henning 2003; Lufkin et al. 2003) usually assume power law profiles for the disk properties, in models often reproducing the minimum mass solar nebula (MMSN) characteristics. This is one of our main motivations to study “type–I” gaseous migration with more realistic proto-planetary disk models. While there are significant uncertainties regarding the nature of angular momentum transport in these disks (Gammie 1996; Glassgold, Najita & Igea 1997; Fromang, Terquem & Balbus 2002; Fleming & Stone 2003; Matsumura & Pudritz 2003) or their overall structure (e.g., Chiang & Goldreich 1997), we believe that a study of how these uncertainties may affect migration time estimates is useful at this time.

Another source of uncertainty in calculations of “type–I” migration rates comes from limitations of spiral density wave theory itself. The theory has been initially developed for infinitely-thin, two-dimensional disks (Goldreich & Tremaine 1980), but three-dimensional effects are important (Tanaka, Takeuchi & Ward 2002). Finite eccentricities and inclinations of embedded proto-planets (Papaloizou & Larwood 2000; Goldreich & Sari 2003; Ogilvie & Lubow 2003) and the potentially magnetized (Terquem 2003) and turbulent nature (Winters, Hawley & Balbus 2003; Papaloizou & Nelson 2003; Nelson & Papaloizou 2003a,b; Papaloizou, Nelson & Snellgrove 2003) of the proto-planetary disk are also ignored in the original theory, while they could all have important consequences. It is probably fair to say that a consensus has yet to be achieved on these various aspects of spiral density wave theory.

Nevertheless, an important motivation for developing a detailed numerical tool to study gaseous migration in proto-planetary disks is the intimate relation, which has become increasingly obvious, between planetary formation and migration on the one hand, and disk evolution on the other hand. Planets form from disk material, migrate by interacting with the disk and at the same time influence strongly the disk long-term evolution by raising torques which affect the disk global angular momentum budget and eventually open gaps in an otherwise smooth gaseous distribution. Therefore, it is likely that a general planet formation theory capable of explaining the diversity of known systems will require following the strongly-coupled evolution of a potentially large ensemble of proto-planets together with the disk in which they formed (see, e.g., Lin & Papaloizou 1986 for early calculations of this type).

In this first investigation of gaseous migration in proto-planetary disks, we focus on “type–I” migration of low-mass proto-planets, ignoring the feedback from Lindblad torques on the disk structure. The radial structure of our disks is not a simple power law but depends upon opacities and other details of accretion. We explore how uncertainties in the structure of proto-planetary disks may affect “type–I” migration rates and how sensitive these rates are to three-dimensional effects. In §2, we describe the characteristics of our T-Tauri disk model, while in §3 the two- and three-dimensional formulations we have adopted for Lindblad torques are summarized. Corotation torques, and our reasons for neglecting them, are discussed briefly in §3.3. We present our results in §4, we discuss possible limitations and extensions of this work in §5 before concluding in §6.

## 2 Disk Models

The present study is exclusively focused on steady-state $`\alpha `$–disk models (see §4.1 for some justification). In the future, however, we intend to study the coupled evolution of planet-disk systems with explicit planetary migration, and this obviously requires being able to follow the evolution of the disk’s properties with time. With this in mind, we have developed a fully time-dependent model for T-Tauri $`\alpha `$–disks. It is largely inspired by the model developed by Hameury et al. (1998) to study the time-dependent evolution of disks around compact objects. The main assumption in our models is that disks are geometrically thin, so that their radial and vertical structures can effectively be treated separately.

### 2.1 Radial Structure and Evolution

Following Hameury et al. (1998), we solve the following radial equations for the conservation of mass, angular momentum, and energy:

$`{\displaystyle \frac{\mathrm{\Sigma }}{t}}`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}(r\mathrm{\Sigma }v_\mathrm{r})+\dot{\mathrm{\Sigma }},`$ (1)
$`j{\displaystyle \frac{\mathrm{\Sigma }}{t}}`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}(r\mathrm{\Sigma }jv_\mathrm{r})+{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}\left({\displaystyle \frac{3}{2}}r^2\mathrm{\Sigma }\nu \mathrm{\Omega }_\mathrm{K}\right)+\dot{J},`$ (2)
$`{\displaystyle \frac{T_\mathrm{c}}{t}}`$ $`=`$ $`{\displaystyle \frac{2(Q^+Q^{}+H)}{C_P\mathrm{\Sigma }}}{\displaystyle \frac{\mathrm{}T_\mathrm{c}}{\mu C_P}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{(rv_\mathrm{r})}{r}}v_\mathrm{r}{\displaystyle \frac{T_\mathrm{c}}{r}},`$ (3)
where $`\mathrm{\Sigma }`$ is the surface column density, $`v_\mathrm{r}`$ is the gas radial velocity in the disk, $`j=(GM_{}r)^{1/2}`$ is the specific angular momentum of material at radius $`r`$ in the disk, $`\mathrm{\Omega }_K=(GM_{}/r^3)^{1/2}`$ is the Keplerian angular velocity, $`M_{}`$ is the mass of the central star (taken to be $`0.5M_{}`$ in the present work), $`T_c`$ is the mid-plane temperature, $`C_P`$ is the heat capacity at constant pressure, $`\mathrm{}`$ is the perfect gas constant, $`\mu `$ is the mean molecular weight and $`\nu `$ is the kinematic viscosity coefficient, for which we adopt a standard $`\alpha `$-prescription (Shakura & Sunyaev 1973). Justifications for the use of an $`\alpha `$–viscosity to describe angular momentum transport in thin accretion disks have been given by Balbus & Papaloizou (1999) in the case of magnetized turbulent transport and by Goodman & Rafikov (2001) in the case of transport by locally-dissipated planetary torques. The gas motion is implicitly described as Keplerian in our formulation of the disk equations. Deviations from strict Keplerianity due to finite radial pressure gradients are calculated separately, when required for the planetary torque calculations, and deviations from Keplerianity due to disk self-gravity are ignored.

The terms $`\dot{\mathrm{\Sigma }}`$ and $`\dot{J}`$ in the first two equations represent additional sources/sinks of mass and angular momentum, respectively, for the disk (both are put to zero in the present study). $`Q^+`$ and $`Q^{}`$ are the local heating and cooling rates per unit surface, respectively. The term $`H`$ accounts for the radial energy flux caused by steep radial gradients, and it is estimated in the framework of the $`\alpha `$ parametrization (see Hameury et al. 1998 for details). This general form of the energy equation accounts for the possibility that the disk is out of thermal equilibrium ($`Q^+Q^{}`$) as can occasionally happen in time-dependent studies (see, e.g., Menou, Hameury & Stehle 1999 for specific examples in dwarf nova disks). For all the steady-state models presented here, however, the radial energy equation reduces to a very good approximation to $`Q^+=Q^{}`$.

Equations (1) and (2) can be combined to obtain the standard viscous diffusion equation for a geometrically thin disk (when $`\dot{\mathrm{\Sigma }}=\dot{J}=0`$; see, e.g., Pringle 1981). The above set of equations is solved on an adaptive numerical grid, to guarantee that regions of the disk characterized by sharp gradients receive a large number of grid points and thus adequate numerical resolution. A number $`N=400`$ grid points has been adopted for all the calculations presented here. The equations are evolved with a fully implicit numerical scheme and a variable time-step. Further details can be found in Hameury et al. (1998). For the present study, steady-state disk models were obtained in the limit $`/t0`$ (in practice, a large enough time-step).

This set of partial differential equations requires a total of six boundary conditions. For a disk with fixed inner and outer radii which is fed, at its outer edge, at a rate $`\dot{M}=\dot{M}_{\mathrm{out}}`$, the boundary conditions at the inner edge ($`n=1`$) and the outer edge ($`n=N`$) of the numerical domain are:

$`\mathrm{\Sigma }|_{n=1}`$ $`=`$ $`ϵ1\mathrm{g}\mathrm{cm}^2,`$
$`r|_{n=1}`$ $`=`$ $`r_{\mathrm{in}}(0.01AU),`$
$`{\displaystyle \frac{T_c}{r}}|_{n=1}`$ $`=`$ $`0,`$
$`\dot{M}|_{n=N}`$ $`=`$ $`\dot{M}_{\mathrm{out}},`$
$`r|_{n=N}`$ $`=`$ $`r_{\mathrm{out}}(200AU),`$
$`{\displaystyle \frac{T_c}{r}}|_{n=N}`$ $`=`$ $`0.`$ (4)
The precise value of $`\mathrm{\Sigma }|_{n=1}`$ is unimportant as long as it is small compared to values in the bulk of the disk.

### 2.2 Vertical Energy Transport

Contrary to Hameury et al. (1998), we do not explicitly solve for radiative-convective energy transport in the vertical direction. Instead, we calculate the local cooling rate $`Q^{}2\sigma T_{\mathrm{eff}}^4`$ (for 2 disk faces) with a one-layer version of the radiative transfer model of Hubeny (1990). According to this model, the disk intrinsic effective temperature, $`T_{\mathrm{eff}}`$, is related to the disk mid-plane temperature, $`T_c`$, by

$$T_{\mathrm{eff}}^4=\frac{4}{3}\frac{\left(T_c^4T_{\mathrm{irr}}^4\right)}{\left(\tau _{\mathrm{ext}}/2+1/\sqrt{3}+1/(3\tau _{\mathrm{abs}})\right)},$$
(5)
where the mid-plane optical depth $`\tau =\kappa \mathrm{\Sigma }/2`$, $`\kappa _{\mathrm{ext}}`$ is the Rosseland-mean extinction opacity (absorption+scattering), $`\kappa _{\mathrm{abs}}`$ is the Rosseland-mean absorption opacity and $`T_{\mathrm{irr}}>0`$ in the presence of external irradiation. Note that in this formulation, $`T_{\mathrm{eff}}`$ accounts only for the disk intrinsic emission (i.e. from viscous dissipation only), while the disk also re-emits all the irradiation flux absorbed locally (as measured by $`T_{\mathrm{irr}}`$). A detailed discussion of this radiative transfer model and of how well it does in the two important optically-thick and optically-thin limits can be found in Hubeny (1990). While Hubeny’s model requires in principle a Planck-mean opacity to be used for $`\tau _{\mathrm{abs}}`$ (in the optically-thin regime), we only use Rosseland-mean opacities here (see §2.3 for a justification). We do not allow for different opacities when considering the disk’s own radiation and the hotter radiation from the central star. Adopting this simple prescription for vertical radiative transfer is an oversimplification of the problem, which is justified a posteriori by satisfactory comparisons to previous studies treating the vertical energy transport problem in more details. Our one-layer approximation neglects vertical transport by convection, whose role is unclear, especially in the context of magnetized accretion disks, which are already subject to a dynamical instability leading to turbulent transport, by construction. We also note that D’Alessio et al. (1998) have estimated that convective transport should be unimportant in all but the very innermost regions of a typical T-Tauri disk.

Viscous dissipation is prescribed in the usual form, with a heating rate per unit area

$$Q^+=\frac{9}{4}\nu \mathrm{\Sigma }\mathrm{\Omega }_K^2,$$
(6)
where the kinematic viscosity, $`\nu =\alpha c_sh`$, is expressed in terms of a dimensionless viscosity parameter, $`\alpha `$ ($`c_s`$ is the local adiabatic sound speed and $`h`$ is the local scale height).

Irradiation by the central star is an important source of heating for the disk at large distances. Rather than solving self-consistently for the magnitude of irradiation based on a detailed model of the disk shape and structure, as in Dubus et al. (1999), we parametrize it with a simple geometrical model. Assuming that the irradiation source is point-like (a reasonable assumption far away from the stellar source), the irradiation flux is given by

$$\sigma T_{\mathrm{irr}}^4=\frac{L_{}(1ϵ)}{4\pi r^2}\frac{h}{r}\left(\frac{d\mathrm{log}h}{d\mathrm{log}r}1\right),$$
(7)
where $`\sigma `$ is the Stefan-Boltzmann constant, $`L_{}`$ is the stellar luminosity, $`ϵ`$ is the disk albedo and $`T_{\mathrm{irr}}`$ should go to zero if the last term in parenthesis becomes negative (disk self-shadowing; see, e.g., Dubus et al. 1999). For $`ϵ=0.5`$, an aspect ratio $`h/r0.1`$, a geometric factor $`0.1`$ (last term in parenthesis in Eq. ) and $`L_{}L_{}`$, this yields an irradiation temperature

$$T_{\mathrm{irr}}90K\left(\frac{r}{1\mathrm{AU}}\right)^{1/2},$$
(8)
which is directly used in equation (5). This prescription artificially exposes disk regions with $`d\mathrm{log}h/d\mathrm{log}r<1`$, that would normally be self-shadowed, to irradiation from the central star, but without serious consequences as we shall see below. Although an improved treatment of disk irradiation would lead to better accuracy, this simple prescription should be sufficient for our purposes in view of other large uncertainties inherent to structural models of proto-planetary disks.

### 2.3 Opacities and Equation of State

We assume that dust (when present) is fully mixed with the gas. Analytic expressions for the Rosseland-mean extinction opacities are taken from Bell et al. (1997; see also Henning & Stognienko 1996). Pollack et al. (1994) have found that Planck-mean and Rosseland-mean extinction opacities are comparable for mixed dust. This justifies our use of Rosseland-mean opacities in Eq. (5) even in the optically-thin regime. Following Pollack et al. (1994; see their Fig. 4b), we also adopt $`\kappa _{\mathrm{abs}}=\kappa _{\mathrm{ext}}`$, since we have found that this is adequate for the range of temperatures at which typical T-Tauri disks are optically-thin.

The equation of state adopted is that of a perfect gas. The mean molecular weight and adiabatic index are typically $`\mu _{H_2}2.4`$ and $`\mathrm{\Gamma }7/5`$ throughout the disk, since hydrogen is everywhere predominantly in molecular form. Other thermodynamical quantities are self-consistently calculated as a function of temperature and density at the disk mid-plane when needed (see Hameury et al. 1998 for details).

## 3 Planetary Torques

In this study, we focus on the low-mass perturber limit and we assume that the feedback on the disk due to the presence of embedded proto-planets can be neglected. We neglect corotation torques in our calculations, and assume for simplicity that all the planets are on circular orbits, without inclination relative to the disk mid-plane. Linear spiral density wave theory has been elaborated by a number of authors (GT80; KP93; Artymowicz 1993a,b; Ward 1997a) and we rely heavily on these works in our formulation of the problem. Thus, we calculate the “torque density” representing waves excited at Lindblad resonances, which are located at radii $`r_m`$ where the disk’s angular velocity $`\mathrm{\Omega }(r)`$ and epicyclic frequency $`\kappa (r)`$ are related to the angular velocity $`\mathrm{\Omega }_p`$ of the planet by $`m\left[\mathrm{\Omega }(r_m)\mathrm{\Omega }_p\right]=\kappa (r_m)`$ for azimuthal harmonics $`m`$.

### 3.1 Torque Density

Following Ward (1997a), the torque density (torque per unit radius) of waves excited in the disk by a planet of mass $`M_p=\mu M_{}`$ orbiting a star of mass $`M_{}`$ at radius $`r_p`$ is

$$\left[\frac{dT}{dr}\left(r\right)\right]_{LR}=\text{sign}(rr_p)\frac{2\mu ^2\mathrm{\Sigma }r_p^4\mathrm{\Omega }_p^4}{r(1+4\xi ^2)\kappa ^2}m^4\psi ^2.$$
(9)
The surface density, $`\mathrm{\Sigma }`$, and epicyclic frequency, $`\kappa `$ (not to be confused with opacity), are evaluated at the radius of interest, $`r`$, and $`\xi mc_s/r\kappa `$, where $`c_s`$ is the adiabatic sound speed. This expression incorporates the torque cutoff *via* the functional dependence of $`\xi (r)`$, $`\psi (r)`$ and $`m(r)`$ (see below).

Note that

$$1+\xi ^2=\frac{(\mathrm{\Omega }\mathrm{\Omega }_p)^2}{(\mathrm{\Omega }\mathrm{\Omega }_p)^2c_s^2/r^2}.$$
(10)
The angular velocity and epicyclic frequency in the gaseous disk are slightly different from their Keplerian values due to finite radial pressure gradients:

$`\mathrm{\Omega }^2(r)`$ $`=`$ $`{\displaystyle \frac{GM_{}}{r^3}}+{\displaystyle \frac{1}{r\rho }}{\displaystyle \frac{}{r}}(\rho c_{\mathrm{iso}}^2),`$ (11)
$`\kappa ^2`$ $``$ $`{\displaystyle \frac{1}{r^3}}{\displaystyle \frac{}{r}}(r^4\mathrm{\Omega }^2),`$
where $`c_{\mathrm{iso}}`$ is the isothermal sound speed. We must include the pressure terms in the disk rotation curve when we compute the torque density, even if we neglect them in the equations for the disk evolution (§2.1). The angular velocity of the satellite, $`\mathrm{\Omega }_p`$, has the Keplerian value $`\sqrt{GM_{}/r_p^3}`$, and the small difference $`\mathrm{\Omega }_p\mathrm{\Omega }(r_p)`$ is the most important cause of the imbalance between the torques exerted on the interior and exterior parts of the disk, and hence of the planet’s migration. Because of the extra radial derivative, $`\kappa `$ is even more sensitive to pressure effects.333There is some danger that $`\kappa ^2`$ may be formally negative near the edges of gaps, which are not considered in the present study.

The dimensionless azimuthal wavenumber $`m`$ is treated as a continuous function of radius. It is computed from the WKBJ dispersion relation by setting $`km/r`$ (since the radial wavenumber vanishes at the Lindblad resonances where the waves are excited):

$$m(r)=\left[\frac{\kappa ^2}{(\mathrm{\Omega }\mathrm{\Omega }_p)^2c_s^2/r^2}\right]^{1/2}.$$
(12)
Where equation (12) predicts imaginary $`m`$, $`dT/dr0`$. Disk self-gravity is neglected in equation (12).

The dimensionless satellite forcing function is

$$\psi \frac{\pi }{2}\left[\frac{1}{m}\left|\frac{db_{1/2}^m}{d\beta }\right|+2\frac{\mathrm{\Omega }}{\kappa }\sqrt{1+\xi ^2}b_{1/2}^m(\beta )\right],$$
(13)
where $`\beta r/r_p`$. The factor $`\mathrm{\Omega }/\kappa `$ does not appear in the original formulation of Ward (1997a) but it is needed so that equation (13) reduces to equation (18) of Goldreich & Tremaine (1980), in the limit $`\xi 1`$. The standard definition of the Laplace coefficient is

$$b_{1/2}^m(\beta )=\frac{2}{\pi }\underset{0}{\overset{\pi }{}}\frac{\mathrm{cos}m\theta d\theta }{\sqrt{12\beta \mathrm{cos}\theta +\beta ^2}}.$$
(14)
This is exact for an infinitely-thin two-dimensional disk, as studied, for instance, by Ward (1997a), and this is the formulation we adopt for 2D Lindblad torques.

To account for three-dimensional effects, we replace the above expression with the following approximation, which includes a correction for the thickness, $`h`$, of the disk

$$b_{1/2}^m(\beta )\frac{2}{\pi \sqrt{\beta }}K_0\left(m\sqrt{\beta 2+\beta ^1+h^2/rr_p}\right),h\frac{c_{\mathrm{iso}}}{\mathrm{\Omega }},$$
(15)
where $`K_0`$ is the modified Bessel function of the second kind of order $`0`$: $`K_0(z)\mathrm{ln}z`$ for $`0<z1`$ and $`K_0(z)e^z\sqrt{\pi /2z}`$ for $`z1`$. Following Goldreich & Tremaine (1980), we have made use of the integral representation

$$K_0(z)=\underset{0}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}(zt)}{\sqrt{t^2+1}}𝑑t$$
(16)
(Abramowitz & Stegun §9.6.21). The appearance of a term proportional to $`h^2`$ in equation (15) comes from expressing the gravitational potential of the perturbing planet as

$$\mathrm{\Phi }_p=\frac{GM_p}{\sqrt{|𝐫𝐫_𝐩|^2+h^2}},$$
(17)
while $`h`$ is set exactly to zero in the two-dimensional theory.

With this expression, we account directly in the perturbing potential for the fact that the disk material subject to the gravitational influence of the proto-planet is actually distributed vertically over a length scale typically $`h`$, the disk scale height. It is clear that this formulation accounts only approximately for three-dimensional effects. In a detailed study, Tanaka et al. (2002) find that the exact vertical distribution of mass in the disk enters explicitly a three-dimensional formulation of spiral density wave theory (see their Eq. , which is derived under a specific vertically isothermal assumption). As we emphasized earlier, the structure of proto-planetary disks, including their vertical structure, is not well known (e.g., Fromang et al. 2002; Fleming & Stone 2003; see also Chiang & Goldreich 1997). In this context, it is unlikely that three-dimensional torques can be calculated very accurately. This is our main motivation for adopting a simple formulation of three-dimensional effects, that should nonetheless capture their most salient features. We note that Papaloizou & Larwood (2000) have adopted a similar formulation, with a vertical “softening length” varying from $`40\%`$ to $`100\%`$ of a disk scale-height, $`h`$ (see also Ward 1989 for a different treatment of vertical thickness). It is unclear what value of the softening length is most appropriate or best reproduces the results of Tanaka et al. (2002) in general, but most cases of interest should be covered by a treatment intermediate between our two- and three-dimensional prescriptions.

The expression for $`\psi `$ includes a derivative with respect to $`\beta =r/r_p`$, and we include the implicit dependence on $`r`$ via $`h^2(r)/r`$ when evaluating this derivative. However, $`m`$ is treated as a constant during the differentiation. To evaluate the term $`db_{1/2}^m(\beta )/d\beta `$ in equation (13), note that $`dK_0(z)/dz=K_1(z)`$ (Abramowitz & Stegun §9.6.27), where $`K_1`$ is the modified Bessel function of the second kind of order $`1`$. Therefore

$$\frac{db_{1/2}^m}{d\beta }(\beta )=(2\beta )^1b_{1/2}^m(\beta )\frac{mK_1\left(m\sqrt{\beta 2+\beta ^1+h^2/rr_p}\right)}{\pi \sqrt{(\beta 1)^2+h^2/r_p^2}}\left[1\beta ^2+\frac{d}{dr}\left(\frac{h^2}{r}\right)\right].$$
(18)
Since $`h=c_{\mathrm{iso}}/\mathrm{\Omega }`$,

$$\frac{d}{dr}\frac{h^2}{r}=\left(\frac{d}{dr}\mathrm{ln}T_c+\frac{2}{r}\right)\frac{h^2}{r}.$$

When $`h|rr_p|r_p`$, then $`\xi 1`$ and $`\psi `$ reduces approximately to a constant:

$$\psi K_1(\kappa /|2A|)+\frac{2\mathrm{\Omega }}{\kappa }K_0(\kappa /|2A|),$$

where $`2Ad\mathrm{\Omega }/d\mathrm{ln}r=(\kappa ^24\mathrm{\Omega }^2)/2\mathrm{\Omega }`$ is the shear rate. Unfortunately, the torque is dominated by the region $`|rr_p|h`$, where this approximation cannot be used. In fact, as $`|rr_p|2h/3`$, $`m`$ diverges and $`\psi 0`$ exponentially. The torque density vanishes at $`|rr_p|2h/3`$: in this region, the flow past the planet is subsonic so there is no wave drag.

### 3.2 Implementation of Planetary Torques in Disk Models

The torque density has a sharply peaked profile, cutting off at a distance $`h`$ from the proto-planet’s location (see, e.g., Ward 1997a). An accurate calculation of the torque density therefore requires that radial profiles in the background disk model be adequately resolved on length-scales smaller than a disk scale height. This is not trivially achieved because the typical grid spacing in our disk models is usually not very much smaller than a disk scale height.

A first method we explored to calculate the torque density was to make use of the adaptive grid. In principle, by requiring that the grid spacing be inversely proportional to the local torque density, one should reach adequate spatial resolution in the disk model (in the vicinity of a proto-planet) for the torque density to be calculated accurately. In practice, we have found that this solution is difficult to implement and often results in numerical instability.

Instead, we have opted for a sub-grid method. Whenever required for a torque density calculation, the radial profiles of relevant quantities in the disk are fitted to cubic splines in the vicinity of the proto-planet’s location. For the torque density calculation itself, cubic spline interpolations of the disk profiles provide estimates of disk properties at distances as close to the proto-planet as required (whose location does not in general coincide with that of a grid point in the disk model). This implementation effectively decouples the accuracy of the torque density calculation from the effective numerical resolution of the background disk model. Generalization to multiple proto-planets embedded in the disk is straightforward.

### 3.3 Corotation Torques

In addition to the Lindblad torques, there are torques at co-orbital corotation resonances. The latter occur very near the semimajor axis of the planet, $`r_Cr_pO(h^2/r^2)`$. Eccentric planets have non-co-orbital corotation resonances as well (*e.g.,* GT80).

The theory of corotation resonances is in a less satisfactory state than that of Lindblad resonances. GT80 derived an expression for the corotation torque due to an individual azimuthal harmonic. In the co-orbital case and for $`m1`$, their expression reduces to

$$\mathrm{\Gamma }_m^{C,\mathrm{GT}}8mK_0(m|1\beta |)^2\left(\frac{GM_p}{r_p}\right)^2\left[\frac{1}{d\mathrm{\Omega }/dr}\frac{d}{dr}\left(\frac{\mathrm{\Omega }\mathrm{\Sigma }}{\kappa ^2}\right)\right]_{r_C},$$
(19)
where $`\beta r/r_p`$ and the sign pertains to the torque on the planet. Pressure forces and the disk thickness were neglected; normally these effects should be unimportant at $`mr/h`$. The dominant contribution to the total torque on an embedded planet, however, comes from large $`m`$. For the Lindblad resonances, GT80 demonstrated that harmonics $`m>r/h`$ are suppressed by pressure even in a two-dimensional disk. This result has been elaborated by later authors as summarized above. Analogous “torque cutoffs” were derived for corotation resonances by Ward (1989) using asymptotic methods, and by KP93 using numerical integration of the linearized equations. The latter concluded that the net torque (summed over $`m`$) from co-orbital corotation resonances is comparable to the residual of Lindblad torques—a result that is suggested by simply truncating the sum over eq. (19) at $`m=r/h`$. Since all harmonics are proportional to the term in square brackets above, however, the corotation torque is very sensitive to the local gradient in surface density—even in sign—and to the *third* derivative of the rotation curve. Coincidentally, the torque vanishes in a keplerian disk with the canonical profile of the minimum-mass solar nebula ($`\kappa =\mathrm{\Omega }\mathrm{\Sigma }r^{3/2}`$). Unfortunately, the torque cutoffs derived by Ward (1989) and KP93 do not agree quantitatively, even though both groups used linear theory. A subset of KP93’s results were confirmed in the 2D limit of a three-dimensional linear analysis by Tanaka *et al.* (2002), who treated the special (and relatively tractable) case of disks that are both radially and vertically isothermal.

Linear theory itself is problematic. The angular momentum deposited at Lindblad resonances is carried away by a wave flux. But the torque at corotation accumulates in a small amount of disk material unless spread by viscous diffusion or by differences in the radial drift rates of the planet and the disk (Ward 1989; Balmforth & Korycansky 2001; Masset 2001, 2002). Hence the corotation resonances are prone to nonlinear saturation. The larger the gradient of $`\kappa ^2/2\mathrm{\Omega }\mathrm{\Sigma }`$ (potential vorticity) at the resonance, the stronger is the torque in linear theory, but the more sensitive to saturation. In a disk with a smooth profile, saturation is not likely if the planet’s mass is much less than both of the gap-opening criteria cited in §4.4.

In summary, we have neglected corotation torques for lack of a simple and reliable parametrization of their effects.

## 4 Results

### 4.1 T-Tauri $`\alpha `$-Disk Models

Although accretion disks in T-Tauri systems are not expected to be in steady state, a quasi-steady assumption may be a reasonable first approximation for all but their outermost regions. Indeed, assuming the disk is freely spreading, the global disk evolution should be determined by comparatively slow viscous processes in the outer regions, while the inner regions are able to follow this evolution via a succession of quasi-steady states (see, e.g., Cannizzo, Lee & Goodman 1990; Ruden & Pollack 1991). Using this general property, Hartmann et al. (1998) have shown that the observational characteristics of pre-main sequence stars can be successfully interpreted with freely-spreading accretion disk models in which the inner regions go through a succession of quasi-steady states with an accretion rate that is slowly decreasing with time. They infer mass accretion rates in the range $`10^7`$$`10^9M_{}`$ yr-1, with typical values $`10^8M_{}`$ yr-1 for $`10^6`$ year old T-Tauri star systems.

Based on this, our reference T-Tauri disk model is a steady-state model with a constant accretion rate $`\dot{M}=10^8M_{}`$ yr-1 throughout. We adopt $`M_{}=0.5M_{}`$ for the mass of the central star and a default value of the viscosity parameter $`\alpha =0.02`$. Disk irradiation by the central star, when included, follows the prescription given in equation (8). Later on, we will vary various disk parameters to study how changes in the background disk model influence the planetary migration process. We detail here the properties of our reference disk model because it illustrates rather well the basic characteristics of all the disk models we have considered.

Figure 1a shows the radial profiles, from $`0.01`$ to $`200`$ AU, of the surface density (solid line), Rosseland-mean optical thickness (short-dashed line), aspect ratio (long-dashed line) and Toomre $`Q_T`$ parameter for self-gravity (dotted line) in the reference model. For each quantity, a profile with and without disk irradiation is shown. Note that the profiles very near the inner edge of the numerical domain are affected by our choice of boundary condition but they become independent of it beyond $`0.02`$ AU or so. The disk is characterized by an aspect ratio rising from $`>\mathrm{\hspace{0.33em}10}^2`$ close to the central star to $`>\mathrm{\hspace{0.33em}0.1}`$ at AU distances and beyond. Disk self-gravity is at most weak, even for the coldest outer regions, at $`r=200`$ AU. The disk makes a transition from an optically-thick regime to an optically-thin regime at a distance of several AUs from the central star.

Figure 1b shows radial profiles of the midplane temperature (solid line) and internal effective temperature, arising from internal viscous dissipation only (dashed line), for the same reference model with and without irradiation. The transition from an optically thick to an optically thin regime is again easily identified at a distance of several AUs from the central star. In the irradiated model, the outer, optically-thin regions of the disk are vertically isothermal (the midplane temperature being equal to the irradiation temperature).

In the outer regions of the non-irradiated model, the midplane temperature (lower solid line) is significantly in excess of the effective temperature (dashed line) because of the inefficiency of the optically-thin gas at radiating the dissipated energy. It is clear from figures 1a and 1b that the effect of disk irradiation is relatively minor and that it only is important for regions of the disk beyond 1 AU from the central star. Note also that regions of opacity transitions in the disk can already be identified in figure 1 as regions where the slopes of many of the profiles shown change. As we will see later, these slope changes have potentially important effects on the planetary migration process.

Our results for the steady-state structure of a typical T-Tauri $`\alpha `$-disk are in good overall agreement with previous comparable studies (see, e.g., D’Alessio et al. 1998; Papaloizou & Terquem 1999; Fromang et al. 2002). In particular, the comparison to D’Alessio et al. (1998) is quite satisfactory given that these authors treated the disk structure in considerably more details than we did. These successful comparisons validate a posteriori a number of simplifying assumptions that we made in our treatment of the disk structure (see §2).

This success should not, however, eclipse the fact that there remains considerable theoretical uncertainty regarding the structure of proto-planetary disks. Perhaps the most important source of uncertainty comes from our ignorance of the nature of angular momentum transport (“viscosity”) in the bulk of proto-planetary disks (Gammie 1996; Glassgold, Najita & Igea 1997; Fromang et al. 2002; Fleming & Stone 2003; Matsumura & Pudritz 2003). Over the past few years, it has also become clear that the role of a super-heated layer of dust is important to explain the observational characteristics of T-Tauri disks (Chiang & Goldreich 1997), an aspect of the problem that has simply been ignored in our treatment.

### 4.2 Tests of Lindblad Torque Calculations

Tests of our implementation of the torque density calculations serve several purposes. We want to recover previously established results, guarantee numerical convergence and identify the main differences between the classical 2D treatment of spiral density wave theory and our simple 3D treatment.

Figure 2 shows results for our 2D implementation, in a format such that they can be directly compared to similar results obtained by Ward (1997a), for a disk with fixed aspect ratio, $`h/r=0.07`$. The inner, outer and differential Lindblad torques (or, equivalently, integrated torque densities) are normalized to the value

$$T_0=\pi \mathrm{\Sigma }_pr_p^4\mathrm{\Omega }_p^2\left(\frac{M_p}{M_{}}\right)^2\left(\frac{h_p}{r_p}\right)^3,$$
(20)
where the subscript $`p`$ refers to the proto-planet’s location.

The background disk is idealized by assuming that the radial profiles of surface density and midplane temperature are exact power laws, with indexes

$$k\frac{d\mathrm{ln}\mathrm{\Sigma }}{d\mathrm{ln}r},l\frac{d\mathrm{ln}T_c}{d\mathrm{ln}r}.$$
(21)
A comparison between our figure 2 and figure 3 of Ward (1997a) shows that we recover his results on both the magnitude of Lindblad torques and their dependence on $`k`$ and $`l`$ with high accuracy. In particular, we confirm that 2D Lindblad torques are essentially independent of the background profile of surface density (as measured by $`k`$).

Figure 3 shows very different results for our simple 3D implementation of Lindblad torques. The 3D torques are reduced by a factor $`3`$-$`5`$ from the 2D values, and their dependence on the background disk profiles is affected. More precisely, there is now a clear dependence on the background surface density profile (as measured by $`k`$). In fact, the weakness of 3D Lindblad torques and their dependence on the background disk properties are such that a sign reversal of the differential Lindblad torque appears possible for sufficiently steep surface density profiles (large $`k`$) and sufficiently shallow midplane temperature profiles (small $`l`$).

We have confirmed that the results shown in figures 2 and 3 are unaffected by changes of numerical resolution in the background disk models, in the range $`N=100`$ to $`1000`$. These tests show adequate numerical convergence of our torque density calculations based on a sub-grid method (see §3.2).

### 4.3 Migration Rates

Given a proto-planet of mass $`M_p`$ embedded in a disk at a radius $`r_p`$, the net torque, $`\dot{J}_p`$, experienced by the proto-planet is calculated as a function of the disk properties in the vicinity of $`r_p`$. In the present work, we neglect the feedback on the disk structure caused by the proto-planet’s influence. This assumption should be reasonably accurate in the low-mass limit. Rather than following the actual migration with time of the proto-planet in the gaseous disk, we estimate the local inward migration time of the proto-planet as

$$T_{\mathrm{mig}}\frac{J_p}{\dot{J}_p}\frac{M_p\sqrt{GM_{}r_p}}{\dot{J}_p},$$
(22)
where we have assumed a low-mass proto-planet again in the expression for the Keplerian angular momentum, $`J_p`$. By populating the entire disk with a sufficient number of identical mass proto-planets, we obtain a radial profile of local migration times. In what follows, we show how such migration time profiles are affected by variations in the background disk model. We have found that populating the disk with 200 identical proto-planets, logarithmically spaced from $`r_{\mathrm{in}}=0.01`$ AU to $`r_{\mathrm{out}}=200`$ AU was enough to obtain accurate migration time profiles in all cases of interest.

Figure 4 shows radial profiles of inward migration time, $`T_{\mathrm{mig}}`$, for Earth-mass proto-planets embedded in our reference T-Tauri disk model, with $`\dot{M}=10^8M_{}`$ yr-1 and $`\alpha =0.02`$ (see Fig. 1 for detailed disk properties). The solid line shows migration times calculated with the 3D formalism, while the short-dashed line shows the same result according to the standard 2D formalism of spiral density wave theory. As expected, migration times are significantly longer in the 3D formulation. The dotted line shows, for the 3D formulation, that the effects of neglecting disk irradiation on the estimated migration times are significant only at distances $`>\mathrm{\hspace{0.33em}1}`$ AU from the central star, where the disk structure becomes sensitive to irradiation. In addition, the long-dashed line shows the profile of migration times obtained with a modified 3D formulation, assuming that the gas is vertically extended over only half a local disk scale height for the torque calculation ($`h0.5h`$ in Eq. ). This confirms that migration times are sensitive to the exact vertical distribution of gas in the disk.

The most striking features of figure 4 are the sudden variations of $`T_{\mathrm{mig}}`$ with radius, at distances $`<\mathrm{\hspace{0.33em}1}`$ AU from the central star. Due to a stronger dependence on the background disk properties (in particular on the surface density profile), radial variations of $`T_{\mathrm{mig}}`$ are much more pronounced in the 3D formulation of Lindblad torques than in the 2D one, even if they are clearly present in both cases. It is also apparent that these variations of $`T_{\mathrm{mig}}`$ occur at specific locations in the disk. They are easily identified as the regions where changes of opacity regime occur in the disk (see the corresponding slope changes in radial profiles of $`\mathrm{\Sigma }`$, $`T_c`$ and other quantities in Fig. 1). These radial variations of $`T_{\mathrm{mig}}`$ (or, equivalently, of the differential Lindblad torque) have their origin in differences between the inner and outer disk properties when the proto-planet is located in the close vicinity of an opacity transition. In these special regions of the disk, with different inner and outer slopes for $`\mathrm{\Sigma }`$ and $`T_c`$, differential Lindblad torques can take values that are quite different from more typical values adopted in the rest of the disk.

This property of differential Lindblad torques is interesting because it suggests that narrow regions (e.g., at $`0.1`$ and $`0.25`$ AU in Fig. 4) corresponding to local maxima of the migration times may be the preferred sites of mutual interactions between proto-planets. Clearly, this possibility and its consequences can only be investigated by actually following the time-dependent migration of a population of proto-planets, and accounting for their mutual gravitational interactions as they approach stalling regions such as the ones identified in figure 4.

Another notable feature of the migration times shown in figure 4 is their overall magnitude. In agreement with previous studies (see, e.g., Fig. 14 of Ward 1997a), we find that the typical migration time of an Earth-mass proto-planet located at a few AU distances from the central star is of order a few $`10^5`$ years, according to the 2D formalism of spiral density wave theory. As we already mentioned in §1, this relatively short timescale has been recognized as a major difficulty for the standard core accretion scenario for giant planet formation because it may not allow enough time for the build up of a substantial gaseous envelope. On the other hand, the typical migration times obtained with the 3D formulation are typically one order of magnitude longer than the 2D ones, and regions in the vicinity of opacity transitions are characterized by migration times approaching $`10^7`$ years for an Earth-mass proto-planet. We view this result as an indication that a proper treatment of 3D effects in Lindblad torque calculations, together with a detailed background disk model, may be essential to the resolution (even if only partial) of this long standing difficulty for the core accretion scenario.

In what follows, we investigate the effect on migration times of variations in the background disk model. Since the 3D formalism is arguably the most appropriate (proto-planetary disks do have a finite thickness) and it is more sensitive than the 2D one to details in the background disk properties, we chose to show results for migration times calculated with the 3D formulation only. We also restrict ourselves to disk models that include the effects of irradiation by the central star on the disk structure.

Figure 5 shows the effect on migration times of varying the viscosity parameter $`\alpha `$ in the disk model. As $`\alpha `$ is reduced, the disk mass is increased (at a fixed mass accretion rate $`\dot{M}=10^8M_{}`$ yr-1), so is the differential Lindblad torque, so that migration times are reduced. A secondary consequence of varying $`\alpha `$ is to increase the values of the radii at which opacity transitions occur in the disk, so that the regions of maximum migration times are shifted outward as well. The relative amplitude of radial variations of $`T_{\mathrm{mig}}`$ remain comparable even as $`\alpha `$ is varied.

Figure 6 illustrates the effect of varying the mass accretion rate $`\dot{M}`$ in a background disk model with the same viscosity parameter $`\alpha =0.02`$ as in our reference model. This time, the overall magnitude of migration times is only slightly affected (disks with larger $`\dot{M}`$ are somewhat more massive), while the profiles are shifted outward as the mass accretion rate is increased. Finally, we note that we have checked directly the linear dependence of $`T_{\mathrm{mig}}`$ with the proto-planet’s mass, as expected for this migration of “type I” (without disk feedback). Consequently, even for a proto-planet of mass $`M_p=10M_{}`$, which in some specific core accretion models is considered as a typical threshold mass for the onset of runaway envelope accretion, migration times in some regions (at sub-AU distances) of our reference T-Tauri disk model are still in excess of $`10^6`$ years.

### 4.4 Conditions for Gap Formation

One of the assumptions made in our derivation of migration times is that the feedback torque exerted on the disk by a proto-planet is negligibly small. While this assumption should be reasonably accurate in the low-mass limit, we can further quantify it by estimating the mass above which the feedback torque becomes strong enough for the proto-planet to open a gap. Following Rafikov (2002; see also Lin & Papaloizou 1986 and Ward & Hourigan 1989), a proto-planet with a mass in excess of the ”viscous mass”

$$M_\mathrm{V}\left(10.34\alpha \right)^{1/2}\left(\frac{h_p}{r_p}\right)^{5/2}M_{}$$
(23)
generates a torque larger than the disk’s own viscous torque and opens a gap. In addition, an independent condition for gap formation is that the proto-planet’s Hill radius exceed the local disk scale-height, which translates into a ”thermal mass” limit

$$M_\mathrm{T}\left(\frac{h_p}{r_p}\right)^3M_{}.$$
(24)
The local condition for gap formation is presumably given by the larger of these two masses.

Figure 7 shows radial profiles of $`M_\mathrm{V}`$ and $`M_\mathrm{T}`$ in three representative T-Tauri $`\alpha `$–disk models for which we have estimated migration times previously. In each panel, the solid line corresponds to the viscous mass, $`M_\mathrm{V}`$, and the dashed line to the thermal mass, $`M_\mathrm{T}`$. The upper panel shows the conditions for gap formation in our reference T-Tauri disk model ($`\dot{M}=10^8M_{}`$ yr-1 and $`\alpha =0.02`$), while the middle and lower panels show these conditions in a high accretion rate model ($`\dot{M}=10^7M_{}`$ yr-1) and in a low $`\alpha `$ model ($`\alpha =10^3`$), respectively.

In each case, one easily identifies the strong dependence of $`M_\mathrm{V}`$ and $`M_\mathrm{T}`$ on the local disk aspect ratio, $`h/r`$, and how the location of the plateau at $``$AU distances varies somewhat from model to model. For both the reference and high $`\dot{M}`$ models, the viscous mass, $`M_\mathrm{V}`$, determines everywhere the conditions for gap formation, while the thermal mass, $`M_\mathrm{T}`$, is relevant everywhere in the low $`\alpha `$ model. Gap formation conditions in models with intermediate values of $`\alpha `$ would be determined by both the viscous and thermal masses, depending on location in the disk.

Except in the very innermost disk regions, effective gap opening conditions correspond to masses in excess of $`10M_{}`$ and reach values larger than 100 Earth masses typically beyond $`0.1`$-$`0.3`$ AU. In addition, it should be noted that the viscous mass threshold (Eq. ) was determined from a 2D formulation of spiral density wave theory. Since a more realistic treatment reveals that 3D Lindblad torques are significantly weaker than 2D ones, a derivation of $`M_\mathrm{V}`$ accounting for the finite thickness of proto-planetary disks should lead to even larger values of the viscous mass for gap formation. Therefore, provided one is interested in proto-planets with masses significantly smaller than the values shown in figure 7, estimates of local migration times ignoring disk feedback should be reasonably accurate, except for the effects of corotation torques, and except insofar as the collective action of many small planets may modify the disk.

## 5 Discussion

Our detailed calculations of Lindblad torques for low-mass proto-planets embedded in steady-state T-Tauri $`\alpha `$–disks reveal strong sensitivities to radial properties of the background disk and its finite thickness. While we have found these results to be robust within the set of calculations we performed, several important limitations should also be emphasized at this point.

We have ignored corotation torques for lack of an adequate parametrization of their effects. Although Lindblad torques are usually thought to dominate migration rates (e.g. KP93), the corotation resonances may be even more sensitive to local conditions and to feedback than the Lindblad ones, as explained in §3.3. Very likely, corotation torques would alter our quantitative results at order unity but would reinforce our main qualitative conclusions, namely, that the migration rate is strongly modulated by opacity transitions and other radially localized features of $`\alpha `$ disks. This will have to be tested in future work.

It is unlikely that our steady-state constant–$`\alpha `$ disk models are accurate representations of T-Tauri disks in general. Indeed, the nature of angular momentum transport in proto-planetary disks is not well understood, even though it determines their structure to a large extent. Gammie (1996) observed that MHD turbulence driven by the magneto-rotational instability (Balbus & Hawley 1991; 1998) should operate in the innermost regions ($`r<\mathrm{\hspace{0.33em}0.1}`$ AU) of T-Tauri disks under standard thermal ionization conditions, but probably not in their colder, outer regions. He proposed instead that MHD-turbulent accretion at large radii proceeds only in surface layers that are sufficiently ionized by cosmic rays. Glassgold et al. (1997) later refined this layered accretion picture, by taking into account the ionization due to X-rays from the central star. Since then, it has become clear that the radial extent of a layered accretion region is sensitive to details in the disk models (Fromang et al. 2002; Matsumura & Pudritz 2003), with weakly-turbulent properties that were not previously anticipated (Fleming & Stone 2003).

Absent any contribution to angular momentum transport from hydrodynamical processes alone (see Balbus 2003 for a discussion), proto-planetary disks may also have to rely on planetary torques for their ”viscous” evolution (Larson 1990). Goodman & Rafikov (2001) studied the propagation and dissipation of spiral density waves in gaseous disks, within the framework of two-dimensional theory. They showed that dissipation from wave breaking is quasi-local and they estimated a magnitude of angular momentum transport $`\alpha 10^3`$$`10^4`$ for a uniform population of Earth-mass proto-planets embedded in a minimum mass solar nebula. It remains to be seen what efficiencies of angular momentum transport can be reached self-consistently in more realistic proto-planetary disk models that include explicitly the integrated torque feedback from a population of embedded proto-planets.

Clearly, these large uncertainties on the structure of proto-planetary disks will have important effects on migration time estimates. We believe, however, that the identification of opacity transitions in the disk as regions where migration rates vary rapidly with radius is not unique to our steady-state $`\alpha `$-disk models. In the specific layered accretion model of Gammie (1996), for instance, opacity transitions cause variations in the slopes of radial solutions, much as in our models, and this may be a general property of models based on an $`\alpha `$-formalism.

What is less clear perhaps is the magnitude of this effect on radial profiles of migration times. In particular, our models have assumed a direct dependence of the opacity on the local values of density and temperature in the disk. In the presence of substantial radial mixing, for instance due to disk turbulence, these transitions could be made more gradual and the amplitude of the rapid radial variations of migration times we have found may be reduced accordingly. The issue of the turbulent nature of the disk is particularly important in view of recent results from fully turbulent MHD simulations for embedded proto-planets (Nelson & Papaloizou 2003b; see also Laughlin, Steinacker & Adams 2003). These simulations show how turbulent density perturbations can dominate over the specific spiral density pattern associated with an embedded proto-planet and lead to a largely random migration process. It remains to be seen whether long-term migration trends will emerge from such simulations and whether they will be sensitive to global radial disk properties such as those exhibited by our $`\alpha `$–disk models. If not, it may be that our results on radial variations of migration times will be mostly relevant to the ”dead zones” of layered accretion regions, where the turbulence is expected to be comparatively much weaker (Fleming & Stone 2003).

Not surprisingly, many of the uncertainties we have just highlighted provide good motivations for additional work on the subject. It would be interesting to include in our disk models a treatment of layered accretion similar to that described by Gammie (1997) and study proto-planet migration in that context. Another important extension of the present work would be to include the feedback on the disk due to planetary torques. This would allow us to study more accurately the conditions for gap formation, as well as the global efficiency of angular momentum transport in the disk associated with a population of embedded proto-planets. Ultimately, we may thus be able to study with some confidence the coupled evolution of a proto-planetary disk and its population of growing planets.

## 6 Conclusion

We have presented detailed calculations of migration times for low-mass proto-planets embedded in steady-state T-Tauri $`\alpha `$-disks, based on an approximate treatment of two- and three-dimensional Lindblad torques. We have emphasized the strong sensitivity of local migration rates to details in the background disk structure and we have argued that regions of opacity transitions may in general be characterized by sudden radial variations of the migration rates. Localized regions with the slowest migration rates may thus be preferred sites for gravitational interactions between proto-planets. In some of our disk models, we have found regions with migration times in excess of $`10^6`$ yrs for a $`10M_{}`$ proto-planet, and within the framework of the core accretion scenario, this could allow for significant mass buildup of a gaseous envelope before the proto-planet is accreted by the central star.

## Acknowledgments

KM was supported by the Celerity Foundation and the European Community RTN Program “The Origin of Planetary Systems” (PLANETS, contract number HPRN-CT-2002–00308). JG was supported in part by NASA Origins grant NAG5-1164.


format: text

text:
# Single Proton Knock-Out Reactions from 24,25,26F

## I Introduction

The structure of neutron-rich nuclei in the $`p/sd`$-shell region has attracted significant interest recently Bro01 . The emergence of the $`N`$ = 16 shell and the simultaneous disappearance of the $`N`$ = 20 shell sparked the interest especially in the oxygen and fluorine isotopes Oza00 ; Uts01 ; Ots01 . Interaction-, reaction-, and charge-changing cross section Oza01  as well as fragmentation cross sections Lei02  have been measured and analyzed. However, a more detailed knowledge of properties of nuclei in this region is essential to understand the nuclear structure along the dripline.

In view of the emergence of the $`N`$ = 16 shell it has been suggested that the normal level structure is already influenced in 23O Kan02 . One would expect that the ground state of 23O is an $`s_{1/2}`$ level with strong single particle character. However, in order to explain the measured momentum distribution it was necessary to modify the core (22O) structure Kan02 . This interpretation is controversial because an analysis of the data based on the many-body shell model and eikonal theory did not confirm the need for any changes of the level structure in 23O Bro03 ; Kan03 .

Recently, it has been shown that one-nucleon knock-out reactions are a useful tool to extract spectroscopic factors Mad01 ; Han01 ; Bro02 . The single proton stripping reactions (24F,23O) and (25F,24O) should be a sensitive tool to study the ground state properties of 23O and 24O. Although experimentally not confirmed the ground state of 24F is most likely a $`3^+`$ state composed of $`\pi d_{5/2}`$ coupled to a $`\nu s_{1/2}`$ Ree99  and the ground-state of 25F is supposed to be the $`\pi d_{5/2}`$ state. Thus the reactions (24F,23O) and (25F,24O) should have strong couplings to the ground state of 23O and 24O, respectively.

We have measured single-proton knock-out cross section of 24F, 25F, and 26F to populate states of 22O, 23O, and 24O. In most cases it its necessary to measure the reaction in coincidence with $`\gamma `$-rays to cleanly identify the populated state Nav98 ; Han01 ; Baz03 . However, the search for excited states in the last bound oxygen isotopes 23O and 24O have been negative Sta02  so $`\gamma `$-ray coincidence measurements are not necessary in the present cases to identify the final state. The knock-out reaction of the heaviest fluorine isotope feasible up-to-date (26F,25O) leads to the particle-unstable nucleus 25O Lan85 .

## II Experimental Setup

The experiments were performed at the Coupled Cyclotron Facility of the National Superconducting Cyclotron Laboratory at Michigan State University Mar01 ; Bau01 ; Mor03 . A 110 MeV/nucleon 48Ca19+ beam was fragmented on a beryllium target. The A1900 fragment separator contained an acrylic wedge degrader at the intermediate image. The momentum acceptance of the A1900 was limited to 1%. Details of the production of the secondary beams of 24F, 25F and 26F are listed in Table 1.

At the final focus of the separator a stack of silicon detectors was located as shown in Figure 1. The stack consisted of three 500 $`\mu `$m thick silicon surface barrier detectors followed by three 5000 $`\mu `$m thick lithium-drifted silicon-diodes. A secondary 146 mg/cm2 thick 12C target was located following the first silicon detector (#0). The setup was equivalent to that used in the measurement of proton knock-out reactions of 8B and 9C End03 . The fluorine and oxygen isotopes were stopped in the first thick silicon detector (#3).

The incoming fluorine isotopes were uniquely identified by the energy loss in the first detector (#0) and the time of flight measured with respect to the cyclotron RF.

The oxygen isotopes were identified by their $`\mathrm{\Delta }EE`$ signals in the second (#2) and third (#3) detector after the carbon target with the additional condition that an energy loss corresponding to oxygen isotopes was deposited in the first $`\mathrm{\Delta }E`$ detector (#1) following the target. A mass identification spectrum was generated according to Shi79 :

$$\text{Mass ID}\mathrm{ln}(\alpha \mathrm{\Delta }E)+(\alpha 1)\mathrm{ln}(E+c\mathrm{\Delta }E)\alpha \mathrm{ln}(300)$$

where $`\alpha =ab\mathrm{\Delta }E/T`$ and $`T`$ is the $`\mathrm{\Delta }E`$ detector thickness in $`\mu `$m. The identification spectrum was calibrated with the small percentage of oxygen isotopes which were directly produced in the A1900 production target and transmitted to the focal plane. These isotopes did not contribute to a potential background for the actual knock-out data, because of the requirement that the knock-out events had to be cleanly identified as incoming fluorine isotope in the first silicon detector (#0). The values of the parameters $`a=1.7715`$, $`b=0.2`$, and $`c=0.5`$, were within the range suggested in Ref. Shi79 . The parameter $`a`$ was the most sensitive to a clean mass separation with deviations of less than 0.5% resulting in significant improvements. This analysis method had previously already been applied to single proton knock-out reactions for neutron-rich carbon isotopes Kry93 .

## III Data Analysis and Results

Figure 2 shows the populated oxygen isotopes following breakup reactions from 24F (top), 25F (center) and 26F (bottom). The Mass ID parameter was normalized according to A = (Mass ID)/5. The small overall energy dependence of the calibration results in a small shift for the 26F data so that mass 24 occurs in channel 118 instead of channel 120. In addition to the one-proton knock-out to the ground-state, which is guaranteed because no bound excited states exist in these isotopes, several oxygen isotopes of lower mass are populated. These isotopes can be either populated from sequential decay of excited states from the one-proton stripping reaction or directly from more dissipative reactions (i.e., 1p1n or deuteron, 1p2n or triton, etc.). In the following we limit the analysis to the oxygen isotopes 20$``$24. Although the structures in the mass identification spectrum indicate the population of lighter oxygen isotopes the production mechanism of these isotopes is even more uncertain. No events of 25O were observed in the reaction of incident 26F confirming that 25O is unbound.

The solid lines in the mass identification spectra of Figure 2 correspond to the sum of Gaussian fits to the individual masses. The cross sections for populating these masses can then be extracted from the area of the fits, the target thickness and the number of incoming beam particles. The number of incoming particles is determined from the events in the silicon stack detectors. Any losses due to channeling and other incomplete charge collection effects in the detectors should be similar for the fluorine and the oxygen isotopes and thus cancel. This was a small effect; the difference between the number of incoming fluorine isotopes measured in the silicon stack and the first $`\mathrm{\Delta }`$-E/ToF spectrum was $`<`$4%. No target-out runs were performed so that the contributions from reactions in the detectors before (#0) and after the target (#1) were estimated from the acceptances of the $`\mathrm{\Delta }E`$ gates set in the analysis. Incoming particles were identified as F isotopes event-by-event by their $`\mathrm{\Delta }E`$ and ToF gate as mentioned in section II. The acceptance of the $`\mathrm{\Delta }E`$ gate was set so that isotopes reacting in the first 300$`\mu m`$ of detector #0 were rejected. Oxygen isotopes produced in the last 200$`\mu m`$ could not be discriminated against. Similarly, fluorine isotopes reacting in the first 200$`\mu m`$ of detector #1 could not be excluded from the oxygen identification condition. The reaction cross sections for single-proton knock-out in silicon were calculated with the model described in the following subsection. The contributions from the reaction in the silicon resulted then in a correction of 30$`\pm 5`$%. This is in qualitative agreement with the measured target out contributions measured for the proton knock-out reactions on 8B and 9C End03 .

The final cross sections for the population of oxygen isotopes following the breakup reactions are listed in Table 2. The uncertainties are determined from the Gaussian fit with the area and the width as free parameters. In addition, the uncertainty of the contribution from the background from reactions in the detectors is added in quadrature.

### III.1 Direct One-Proton Knock-Out Reactions to 23,24,25O

The dominant configuration of the valence proton in the projectiles of 24,25,26F is $`\pi d_{5/2}`$ according to the shell model and the actual spectroscopic factors to the ground states calculated in the universal $`sd`$ (USD) interaction Bro88 , given in Table 3, are very close to unity. The listed single-particle cross sections for the removal of this $`d_{5/2}`$ proton were calculated in the eikonal theory Mad01 ; Tos99 ; Han03 . Assuming that the nuclei 23O and 24O have no bound excited levels Sta02 , the measured inclusive cross section should correspond to the product of the single proton knock-out to the ground state and the spectroscopic factor. The single particle ($`\sigma _{sp}`$) and the calculated single proton knock-out cross sections ($`\sigma _{calc}`$) are compared to the experimental values ($`\sigma _{exp}`$) in Table 3.

For this calculation the proton separation energies were taken to be 14.4, 15.0 and 16.1 MeV for 24,25,26F, respectively. The rms matter radii were deduced from experimental interaction cross-sections Oza01  as 3.20(4) (23O) and 3.19(13) fm (24O) for the core and 2.32 fm for the 12C target. An uncertainty of 0.1 fm in the radii translates into an uncertainty of 0.9 mb for the single-particle cross-section. The Woods-Saxon parameters for the calculation of the wave function for the relative proton-core motion were taken as the “standard set” of ref. Bro02 . This work deals with neutron and proton states with separation energies up to 19 MeV and strongly suggests that the eikonal calculation furnishes absolute spectroscopic factors, i.e. relating to physical occupancies of the shell-model states. These, however, turned out to be systematically reduced by a factor 0.5$``$0.6 relative to a many-body shell model based on effective interactions Bro02 ; Han03 . This effect has been attributed to contributions from the nucleon-nucleon interaction not contained in the standard shell model. Viewed in this light, and with experimental errors and theoretical uncertainties taken into account, the resulting theoretical cross sections agree well with the measured values.

The strong population of the 23O ground state in the (24F,23O) reaction yields additional evidence that the usual shell-model calculations give a good description of nuclei in this region Bro03 . In fact, even the simplest picture, that of the ground state of 23O described as a $`[0d_{5/2}^6]1s_{1/2}`$ configuration, gives a near-quantitative description, both of the present data and of the 23O neutron knock-out reactions Kan02 ; Sau00 . The spin of 24F, although most likely 3-, is not essential for the argument, given that the proton in any case will be in the $`0d_{5/2}`$ state.

For the population of the ground state of 25O in the proton knock-out from 26F only an upper limit can be extracted because 25O is unbound and decays to 24O. The ground state of 24O can also be populated from the decay of excited states in 25O. Thus, the measured cross section of 4.1$`\pm 1.4`$ mb correspond to a sum of the contributions from the decay of the ground-state and excited states of 25O (see also Figure 3(c) and the following section).

### III.2 Population of Other Oxygen Isotopes

The cross section for the population of other oxygen isotopes in these breakup reactions can only be estimated because several production mechanisms can contribute. In addition to the single-proton knock-out followed by subsequent neutron emission, direct population of the ground state as well as bound excited states have to be taken into account.

Within the eikonal approach one can calculate the single-particle distribution due to contributions from the knock-out of protons from the $`p`$-shell. The spectroscopic strength distribution as calculated within the shell model using the universal $`sd`$ (USD) interaction Bro88  is shown on the left side of Figure 3. The single-particle cross section calculated in the eikonal model depends on the effective proton separation energy which is the sum of the proton separation energy and the excitation energy of the specific state. The proton separation energies for 24F, 25F, and 26F are 14.4 MeV, 15.0 MeV, and 16.1 MeV, respectively, Figure 4 shows that the proton-removal cross sections are very similar for the three fluorine isotopes. In order to extract the total proton-removal cross sections from the $`p`$-shell the spectroscopic distributions were folded with the excitation-energy dependent stripping cross section of Figure 4. The values for the $`p`$-shell proton-removal cross sections are 25.7 mb, 25.2 mb and 25.0 mb for 24F, 25F, and 26F, respectively, and are listed (p,$`\sigma _{calc}`$) in Figure 3.

These populations are located above the one-, two-, three-, and even the four-neutron evaporation channel. The decay branches to the individual neutron channels depend on the details of the excited states in these nuclei. As mentioned above other processes also contribute to the population of these channels, so that calculating these decay branches is not useful. The only quantitative statement that can be extracted from the present data is that the calculated $`p`$-shell proton-removal cross section should not exceed the summed cross section of the open neutron channels.

For the reaction (24F,23O) the 1-n to 3-n channels have to be taken into account. The threshold for the 4-n channel (19O) at 21 MeV is located above the population for the removal of protons from the $`p`$-shell. The sum of the cross sections for these three neutron decay channels (population of 22-20O) as listed in Table 2 is 39.7$`\pm `$3.0 mb, which has to be compared to a total proton-removal cross section from the $`p`$-shell of 25.7 mb. These cross sections are also indicated in Figure 3.

The sum of the population of the 1-n to 3-n channel (23-21O) in the reaction (25F,24O) is 32.4$`\pm `$2.2 mb. This value is only slightly larger than the calculated proton-removal cross section of 25.2 mb. However, the threshold for the 4-n channel (20O) is located at 17.1 MeV, which is below the highest energies of the distribution of $`p`$-shell proton removal so that it also could contribute.

The reaction (26F,25O) has to be treated differently because 25O is unstable and thus the contribution from $`p`$-shell removal can not be distinguished from the removal of the $`d_{5/2}`$ proton. The total calculated single proton-removal cross section from $`d`$\- and $`p`$-wave knock-out is 32.6 mb. The open channels as indicated in Figure 3 are the 1-n through 4-n decay, which corresponds to a cross section of 38.5$`\pm `$5.0 mb. It might be even possible that the 5-n channel (20O) with a separation energy of 17 MeV can contribute.

For all three reactions studied, the predicted $`p`$-shell removal cross sections are consistent with the measured populations of the open neutron decay channels.

It is interesting to note that the widely used empirical fragmentation model EPAX Sum90 ; Sum00  yields cross sections of the order of 10$``$20 mb for the population of these oxygen isotopes in agreement with the values listed in Table 2.

The recent measurements of fragmentation of unstable neutron-rich oxygen isotopes also included single proton knock-out cross sections Lei02 . These measurements were performed at significantly higher energies ($`>`$500 MeV/nucleon) and yielded cross sections of $``$20 mb. These measurements were not performed as close to the dripline as the current experiment and the populated nitrogen isotopes all had bound excited states, so a direct comparison to the present analysis is not possible. Only the reaction of the most neutron-rich projectile 21O populates a nucleus with no known bound excited states (20N). Unfortunately no cross section is quoted for this reaction Lei02 .

## IV Conclusions

The results of the single-proton knock-out reactions were analyzed within the a theoretical approach combining many-body shell model and eikonal theory. The cross section for the reaction (24F,23O) is consistent with a single particle description of these nuclei, while the cross sections for the reaction (25F,24O) shows a suppression of $``$50%. However, experiments with higher precision are necessary in order to infer a change of the structure. The knock-out reactions on 24F and 25F are special because of the non-existence of bound excited states in 23O and 24O, respectively. The present method promises to be a powerful tool to extract spectroscopic information for nuclei at the very edge of the driplines because they typically do not have any bound excited states. Other nuclei will require more sophisticated experimental arrangements, including $`\gamma `$-ray coincidences.

## V Acknowledgments

This work has been supported by the National Science Foundation grant number PHY01-10253 and the PHY00-70911.


text:
# 1 Abstract

## 1 Abstract

Random matrix ensembles (RME) of quantal statistical Hamiltonian operators, e.g. Gaussian random matrix ensembles (GRME) and Ginibre random matrix ensembles (Ginibre RME), had been applied in literature in study of following quantal statistical systems: molecular systems, nuclear systems, disordered materials, random Ising spin systems, quantal chaotic systems, and two-dimensional electron systems (Wigner-Dyson electrostatic analogy). Measures of quantum chaos and quantum integrability with respect to eigenergies of quantal systems are defined and calculated. Information functional is defined as negentropy (opposite of entropy or minus entropy). Entropy is neginformation (opposite of information or minus information. The distribution functions for the random matrix ensembles are derived from the maximum entropy principle.

## 2 Introduction

The Random Matrix Theory RMT considers quantum Hamiltonian operators $`H`$ that are random matrix variables in given basis of eigenfunctions. Their matrix elements $`H_{ij}`$ are independent random scalar variables . There were studied among others the following Gaussian Random Matrix ensembles GRME: orthogonal GOE, unitary GUE, symplectic GSE, as well as circular ensembles: orthogonal COE, unitary CUE, and symplectic CSE. The choice of ensemble is based on quantum symmetries ascribed to the Hamiltonian $`H`$. The Hamiltonian $`H`$ acts on quantum space $`V`$ of eigenfunctions. It is assumed that $`V`$ is $`N`$-dimensional Hilbert space $`V=𝐅^N`$, where the real, complex, or quaternion field $`𝐅=𝐑,𝐂,𝐇`$, corresponds to GOE, GUE, or GSE, respectively. If the Hamiltonian matrix $`H`$ is hermitean $`H=H^{}`$, then the probability density function of $`H`$ reads:

$`f_{}(H)=𝒞_{H\beta }\mathrm{exp}[\beta {\displaystyle \frac{1}{2}}\mathrm{Tr}(H^2)],`$ (1)
$`𝒞_{H\beta }=({\displaystyle \frac{\beta }{2\pi }})^{𝒩_{H\beta }/2},`$
$`𝒩_{H\beta }=N+{\displaystyle \frac{1}{2}}N(N1)\beta ,`$
$`{\displaystyle f_{}(H)𝑑H}=1,`$
$`dH={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{ji}{\overset{N}{}}}{\displaystyle \underset{\gamma =0}{\overset{D1}{}}}dH_{ij}^{(\gamma )},`$
$`H_{ij}=(H_{ij}^{(0)},\mathrm{},H_{ij}^{(D1)})𝐅,`$
where the parameter $`\beta `$ assume values $`\beta =1,2,4,`$ for GOE($`N`$), GUE($`N`$), GSE($`N`$), respectively, and $`𝒩_{H\beta }`$ is number of independent matrix elements of hermitean Hamiltonian $`H`$. The Hamiltonian $`H`$ belongs to Lie group of hermitean $`N\times N𝐅`$-matrices, and the matrix Haar’s measure $`dH`$ is invariant under transformations from the unitary group U($`N`$, F). The eigenenergies $`E_i,i=1,\mathrm{},N`$, of $`H`$, are real-valued random variables $`E_i=E_i^{}`$. It was Eugene Wigner who firstly dealt with eigenenergy level repulsion phenomenon studying nuclear spectra . RMT is applicable now in many branches of physics: nuclear physics (slow neutron resonances, highly excited complex nuclei), condensed phase physics (fine metallic particles, random Ising model \[spin glasses\]), quantum chaos (quantum billiards, quantum dots), disordered mesoscopic systems (transport phenomena), quantum chromodynamics, quantum gravity, field theory.

## 3 The Ginibre ensembles

Jean Ginibre considered another example of GRME dropping the assumption of hermiticity of Hamiltonians thus defining generic $`𝐅`$-valued Hamiltonian $`K`$ . Hence, $`K`$ belong to general linear Lie group GL($`N`$, F), and the matrix Haar’s measure $`dK`$ is invariant under transformations form that group. The distribution of $`K`$ is given by:

$`f_𝒦(K)=𝒞_{K\beta }\mathrm{exp}[\beta {\displaystyle \frac{1}{2}}\mathrm{Tr}(K^{}K)],`$ (2)
$`𝒞_{K\beta }=({\displaystyle \frac{\beta }{2\pi }})^{𝒩_{K\beta }/2},`$
$`𝒩_{K\beta }=N^2\beta ,`$
$`{\displaystyle f_𝒦(K)𝑑K}=1,`$
$`dK={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\gamma =0}{\overset{D1}{}}}dK_{ij}^{(\gamma )},`$
$`K_{ij}=(K_{ij}^{(0)},\mathrm{},K_{ij}^{(D1)})𝐅,`$
where $`\beta =1,2,4`$, stands for real, complex, and quaternion Ginibre ensembles, respectively. Therefore, the eigenenergies $`Z_i`$ of quantum system ascribed to Ginibre ensemble are complex-valued random variables. The eigenenergies $`Z_i,i=1,\mathrm{},N`$, of nonhermitean Hamiltonian $`K`$ are not real-valued random variables $`Z_iZ_i^{}`$. Jean Ginibre postulated the following joint probability density function of random vector of complex eigenvalues $`Z_1,\mathrm{},Z_N`$ for $`N\times N`$ Hamiltonian matrices $`K`$ for $`\beta =2`$ :

$`P(z_1,\mathrm{},z_N)=`$ (3)
$`={\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{\pi j!}}{\displaystyle \underset{i<j}{\overset{N}{}}}|z_iz_j|^2\mathrm{exp}({\displaystyle \underset{j=1}{\overset{N}{}}}|z_j|^2),`$
where $`z_i`$ are complex-valued sample points ($`z_i𝐂`$).

We emphasize here Wigner and Dyson’s electrostatic analogy. A Coulomb gas of $`N`$ unit charges moving on complex plane (Gauss’s plane) C is considered. The vectors of positions of charges are $`z_i`$ and potential energy of the system is:

$$U(z_1,\mathrm{},z_N)=\underset{i<j}{}\mathrm{ln}|z_iz_j|+\frac{1}{2}\underset{i}{}|z_i^2|.$$
(4)
If gas is in thermodynamical equilibrium at temperature $`T=\frac{1}{2k_B}`$ ($`\beta =\frac{1}{k_BT}=2`$, $`k_B`$ is Boltzmann’s constant), then probability density function of vectors of positions is $`P(z_1,\mathrm{},z_N)`$ Eq. (3). Therefore, complex eigenenergies $`Z_i`$ of quantum system are analogous to vectors of positions of charges of Coulomb gas. Moreover, complex-valued spacings $`\mathrm{\Delta }^1Z_i`$ of complex eigenenergies of quantum system:

$$\mathrm{\Delta }^1Z_i=Z_{i+1}Z_i,i=1,\mathrm{},(N1),$$
(5)
are analogous to vectors of relative positions of electric charges. Finally, complex-valued second differences $`\mathrm{\Delta }^2Z_i`$ of complex eigenenergies:

$$\mathrm{\Delta }^2Z_i=Z_{i+2}2Z_{i+1}+Z_i,i=1,\mathrm{},(N2),$$
(6)
are analogous to vectors of relative positions of vectors of relative positions of electric charges.

The eigenenergies $`Z_i=Z(i)`$ can be treated as values of function $`Z`$ of discrete parameter $`i=1,\mathrm{},N`$. The ”Jacobian” of $`Z_i`$ reads:

$$\mathrm{Jac}Z_i=\frac{Z_i}{i}\frac{\mathrm{\Delta }^1Z_i}{\mathrm{\Delta }^1i}=\mathrm{\Delta }^1Z_i.$$
(7)
We readily have, that the spacing is an discrete analog of Jacobian, since the indexing parameter $`i`$ belongs to discrete space of indices $`iI=\{1,\mathrm{},N\}`$. Therefore, the first derivative with respect to $`i`$ reduces to the first differential quotient. The Hessian is a Jacobian applied to Jacobian. We immediately have the formula for discrete ”Hessian” for the eigenenergies $`Z_i`$:

$$\mathrm{Hess}Z_i=\frac{^2Z_i}{i^2}\frac{\mathrm{\Delta }^2Z_i}{\mathrm{\Delta }^1i^2}=\mathrm{\Delta }^2Z_i.$$
(8)
Thus, the second difference of $`Z`$ is discrete analog of Hessian of $`Z`$. One emphasizes that both ”Jacobian” and ”Hessian” work on discrete index space $`I`$ of indices $`i`$. The spacing is also a discrete analog of energy slope whereas the second difference corresponds to energy curvature with respect to external parameter $`\lambda `$ describing parametric “evolution” of energy levels . The finite differences of order higher than two are discrete analogs of compositions of ”Jacobians” with ”Hessians” of $`Z`$.

The eigenenergies $`E_i,iI`$, of the hermitean Hamiltonian $`H`$ are ordered increasingly real-valued random variables. They are values of discrete function $`E_i=E(i)`$. The first difference of adjacent eigenenergies is:

$$\mathrm{\Delta }^1E_i=E_{i+1}E_i,i=1,\mathrm{},(N1),$$
(9)
are analogous to vectors of relative positions of electric charges of one-dimensional Coulomb gas. It is simply the spacing of two adjacent energies. Real-valued second differences $`\mathrm{\Delta }^2E_i`$ of eigenenergies:

$$\mathrm{\Delta }^2E_i=E_{i+2}2E_{i+1}+E_i,i=1,\mathrm{},(N2),$$
(10)
are analogous to vectors of relative positions of vectors of relative positions of charges of one-dimensional Coulomb gas. The $`\mathrm{\Delta }^2Z_i`$ have their real parts $`\mathrm{Re}\mathrm{\Delta }^2Z_i`$, and imaginary parts $`\mathrm{Im}\mathrm{\Delta }^2Z_i`$, as well as radii (moduli) $`|\mathrm{\Delta }^2Z_i|`$, and main arguments (angles) $`\mathrm{Arg}\mathrm{\Delta }^2Z_i`$. $`\mathrm{\Delta }^2Z_i`$ are extensions of real-valued second differences:

$$\mathrm{\Delta }^2E_i=E_{i+2}2E_{i+1}+E_i,i=1,\mathrm{},(N2),$$
(11)
of adjacent ordered increasingly real-valued eigenenergies $`E_i`$ of Hamiltonian $`H`$ defined for GOE, GUE, GSE, and Poisson ensemble PE (where Poisson ensemble is composed of uncorrelated randomly distributed eigenenergies) . The Jacobian and Hessian operators of energy function $`E(i)=E_i`$ for these ensembles read:

$$\mathrm{Jac}E_i=\frac{E_i}{i}\frac{\mathrm{\Delta }^1E_i}{\mathrm{\Delta }^1i}=\mathrm{\Delta }^1E_i,$$
(12)
and

$$\mathrm{Hess}E_i=\frac{^2E_i}{i^2}\frac{\mathrm{\Delta }^2E_i}{\mathrm{\Delta }^1i^2}=\mathrm{\Delta }^2E_i.$$
(13)
The treatment of first and second differences of eigenenergies as discrete analogs of Jacobians and Hessians allows one to consider these eigenenergies as a magnitudes with statistical properties studied in discrete space of indices. The labelling index $`i`$ of the eigenenergies is an additional variable of ”motion”, hence the space of indices $`I`$ augments the space of dynamics of random magnitudes.

## 4 The Maximum Entropy Principle

In order to derive the probability distribution in matrix space we apply the maximum entropy principle:

$$\mathrm{max}\{S_\beta (f_𝒳):1=1,_𝒳=U_\beta \},$$
(14)
which yields:

$$\mathrm{max}\{S_\beta (f_𝒳):f_𝒳(X)𝑑X=1,_𝒳(X)f_𝒳(X)𝑑X=U_\beta \},$$
(15)
where $`X=H`$ or $`X=K`$ for Gaussian or Ginibre ensembles, respectively, and $`_𝒳(X)=\frac{1}{2}\mathrm{Tr}(X^{}X)`$. The maximization of entropy $`S_\beta (f_𝒳)=(k_B\mathrm{ln}f_𝒳(X))f_𝒳(X)𝑑X`$ under two additional constraints of normalization of the probability density function, and of equality of its first momentum and intrinsic energy, is equivalent to the minimization of the following functional $`(f_𝒳)`$ with the use of Lagrange multipliers $`\alpha _1,\beta _1`$:

$`\mathrm{min}\{(f_𝒳)\},`$ (16)
$`(f_𝒳)={\displaystyle (k_B\mathrm{ln}f_𝒳(X))f_𝒳(X)𝑑X}+\alpha _1{\displaystyle f_𝒳(X)𝑑X}`$
$`+\beta _1{\displaystyle _𝒳(X)f_𝒳(X)𝑑X}.`$
It follows, that the first variational derivative of $`(f_𝒳)`$ must vanish:

$$\frac{\delta (f_𝒳)}{\delta f_𝒳}=0,$$
(17)
which produces:

$$k_B(\mathrm{ln}f_𝒳(X)+1)+\alpha _1+\beta _1_𝒳(X)=0,$$
(18)
and equivalently:

$`f_𝒳(X)=𝒞_{X\beta }\mathrm{exp}[\beta _𝒳(X)]`$ (19)
$`𝒞_{X\beta }=\mathrm{exp}[(\alpha _1+1)k_B^1],\beta =\beta _1k_B^1.`$
The variational principle of maximum entropy does not force additional condition on functional form of $`_𝒳(X)`$. The quantum statistical information functional $`I_\beta `$ is the opposite of entropy:

$$I_\beta (f_𝒳)=S_\beta (f_𝒳)=(+k_B\mathrm{ln}f_𝒳(X))f_𝒳(X)𝑑X.$$
(20)
Information is negentropy, and entropy is neginformation. The maximum entropy principle is equivalent to minimum information priciple.


text:
# Untitled Document

This article has been replaced by gr-qc/0412011