text:
# Some generic aspects of bosonic excitations in disordered systems

\[

## Abstract

We consider non-interacting bosonic excitations in disordered systems, emphasising generic features of quadratic Hamiltonians in the absence of Goldstone modes. We discuss relationships between such Hamiltonians and the symmetry classes established for fermionic systems. We examine the density $`\rho (\omega )`$ of excitation frequencies $`\omega `$, showing how the universal behavior $`\rho (\omega )\omega ^4`$ for small $`\omega `$ can be obtained both from general arguments and by detailed calculations for one-dimensional models.

\]

An understanding of universal properties of excitations in disordered systems occupies a central place in condensed matter physics. Much of the work in this area has focussed on systems of non-interacting fermions, as models for quasiparticles in disordered conductors or superconductors with interactions treated at the mean-field level. Classification of such systems according to symmetry provides an important starting point, and in many instances the possibilities are the three represented by the Wigner-Dyson random matrix ensembles . A recent development, however, has been the recognition that there exist symmetry classes additional to those of the Wigner-Dyson ensembles. These additional symmetry classes arise in fermionic systems which have either chiral symmetry, as for tight-binding models on bipartite lattices with only off-diagonal disorder , or particle-hole symmetry, as for the Bogoliubov - de Gennes Hamiltonian in disordered superconductors . A characteristic feature of both cases is that eigenstates at positive and negative energies are related in pairs and zero energy emerges as a special point in the spectrum.

It is natural to anticipate that similar mathematical structure may be important for the theory of non-interacting bosonic excitations or classical harmonic modes in disordered systems. Our aim in the following is to examine how far this is the case and what consequences it has. Studies of systems of this type have a long history, with celebrated early work by Dyson  on the dynamics of a disordered chain of masses and springs, and applications which include phonons in disordered solids  and spin waves in random magnets . An obvious parallel between between excitations in these systems and those in fermionic systems belonging to the one of the additional symmetry classes is that bosonic excitations arise in positive and negative frequency pairs, and zero is a special point in the frequency spectrum. A second, more formal parallel is that random magnets, for example, in common with superconductors, may give rise to quadratic Hamiltonians containing terms that annihilate and create particle pairs. There are also clear differences between systems with bosonic and fermionic excitations, which prevent an elementary transcription of established ideas. Most importantly, while the exclusion principle guarantees that a fermionic system with a quadratic Hamiltonian has a ground state, for a corresponding bosonic system constraints must be imposed to ensure that it is stable. As a consequence, even to construct a phenomenological treatment such as random matrix theory for bosonic excitations, it is necessary to keep in mind their origins in an appropriate non-linear problem. In turn, a distinction arises that is specific to bosonic excitations, between those which are Goldstone modes and those which are not. We focus below on the latter and discuss elsewhere systems with continuous symmetry .

Any quadratic bosonic Hamiltonian can be written in the forms

$`H`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{N}{}}}\left[M_{ij}p_ip_j+K_{ij}x_ix_j+C_{ij}(x_ip_j+p_jx_i)\right]`$ (1)
$``$ $`{\displaystyle \underset{i,j=1}{\overset{2N}{}}}_{ij}\xi _i\xi _j{\displaystyle \frac{1}{2}}(𝐚^{}𝐚)\left(\begin{array}{cc}\mathrm{\Gamma }& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& \mathrm{\Gamma }^T\end{array}\right)\left(\begin{array}{c}𝐚\\ 𝐚^{}\end{array}\right).`$ (2)
Here, $`x_i`$ and $`p_i`$ are the coordinates and momenta of the oscillators, $`\xi _i=p_i`$, $`\xi _{N+i}=x_i`$ with $`1iN`$, and the vectors $`𝐚^{}`$, $`𝐚`$ have as entries bosonic creation and annihilation operators $`a_i^{},a_i=(x_i\pm ip_i)/\sqrt{2}`$. The matrix $``$ is real symmetric, $`\mathrm{\Gamma }`$ is Hermitian and $`\mathrm{\Delta }`$ is symmetric. The condition for time-reversal invariance is that $`C=0`$ or, equivalently, that $`\mathrm{\Gamma }`$ and $`\mathrm{\Delta }`$ are real.

We are concerned with spectral properties of models of this kind which are generic when $``$ is random. One anticipates that these will be found at small excitation frequencies $`\omega _i`$, and we concentrate particularly on their average density $`\rho (\omega )=N^1\delta (\omega \omega _i)`$ at small $`\omega `$. Indeed, the final expression of Eq. (1) illustrates the parallel between bosonic models and the Bogoliubov - de Gennes Hamiltonian, from which universal behaviour has been derived for the density of quasiparticle states in disordered, gapless superconductors . In fermionic systems such behaviour arises essentially from the interplay of disorder and symmetry, and does not depend on details of the distribution of Hamiltonian matrix elements. By contrast, for the bosonic systems we discuss, we show that the requirements of stability impose features on the distribution of $``$ which determine the form of $`\rho (\omega )`$ at small $`\omega `$.

A simple argument leading to this conclusion has been given previously . In outline it is as follows. Recall that the Hamiltonian of Eq. (1) is characterised not only by the oscillator frequencies $`\omega _i`$ but also by the eigenvalues $`\kappa _n`$ of the matrix $``$, which we refer to as oscillator stiffnesses; let $`d(\kappa )=(2N)^1\delta (\kappa \kappa _n)`$ be their average density. From a discussion of the curvature distribution at absolute minima of random functions of one variable, it is suggested  that $`d(\kappa )\kappa ^{3/2}`$ for small $`\kappa `$. Then, using the relation $`\omega _i^2=\kappa _1\kappa _2`$ which holds in a single-mode system, one arrives at $`\rho (\omega )d(\omega ^2)\omega \omega ^4`$ for small $`\omega `$. A weakness at both steps in this argument is that one degree of freedom is treated in isolation. We show here how the same behaviour emerges without such a restriction.

In a stable system, stiffnesses are positive and frequencies are real. It is helpful to introduce a description that guarantees this property. To this end, write $``$ as a square of a real matrix $`Q`$, in the form $`=Q^TQ`$ (possible provided all $`E_n0`$). Also, note that frequencies are the eigenvalues of an auxilliary matrix $`^{}=\sigma _2`$, where $`\sigma _2=\sigma _yI_N`$, $`\sigma _y`$ is the usual Pauli matrix and $`I_N`$ is the $`N\times N`$ identity matrix . Since the eigenvalues of $`\sigma _2Q^TQ`$ coincide with those of $`\mathrm{\Omega }=Q\sigma _2Q^T`$, which is Hermitian and antisymmetric, frequencies are real and occur in pairs $`\pm \omega _i`$. While there is in general no simple relation between stiffnesses and frequencies, several special cases represent important exceptions. If (in Eq. (1)) $`M=I_N`$ and $`C=0`$ (as for vibrational problems with all masses equal), then $`\omega _i=\pm \sqrt{\kappa _i}`$, where the non-trivial stiffnesses appearing here are the eigenvalues of $`K`$; whereas if $`M=K`$ and $`C=0`$ (as for a random-bond ferromagnet), then $`\omega _i=\pm \kappa _i`$.

Moving beyond these special cases, a straightforward approach to multimode problems comes from concentrating on disorder realisations in which one stiffness, say $`\kappa _1`$, is much smaller than all others. The probability density for such disorder realisations varies as $`\kappa _1^{3/2}`$, by the arguments of Ref . In the limiting case $`\kappa _1=0`$, two eigenvalues of $`\mathrm{\Omega }`$ vanish (since $`\mathrm{det}()=\mathrm{det}(\mathrm{\Omega })`$, and since the eigenvalues of $`\mathrm{\Omega }`$ are paired), while for $`0<\kappa _1\kappa _n`$, $`n1`$, perturbation theory yields a pair of excitation frequencies $`\omega _i\pm \sqrt{\kappa _1}`$. As a result, $`\rho (\omega )\omega ^4`$ for $`\omega \omega _p`$, where $`\omega _p`$ is the lowest excitation frequency in a typical sample. More generally, consider a macroscopic system in which low frequency modes are localized with a localisation length which remains finite as $`\omega 0`$. Treating each localization volume independently, we obtain $`\rho (\omega )\omega ^4`$ for $`\omega \omega _p`$, where $`\omega _p`$ is in this case the lowest excitation frequency in a typical localization volume.

Turning to detailed calculations, it is useful to construct a formulation in which $`Q`$ appears linearly, by considering the enlarged matrices

$$\stackrel{~}{}=\left(\begin{array}{cc}0& Q\\ Q^T& 0\end{array}\right)\text{and}\stackrel{~}{}^{}=\left(\begin{array}{cc}0& Q\\ \sigma _2Q^T& 0\end{array}\right).$$
(3)
Clearly, the eigenvalues of $``$ and $`^{}`$ are squares of those of $`\stackrel{~}{}`$ and $`\stackrel{~}{}^{}`$ respectively. The $`2\times 2`$ structure of $`\stackrel{~}{}`$ displayed in Eq (3) is interesting partly because it is the defining feature of the chiral symmetry classes, studied previously in a variety of contexts , while matrices with the structure of $`\stackrel{~}{}^{}`$ constitute a new, chiral bosonic problem. Whereas much past work on chiral symmetry classes has taken the elements of $`Q`$ to be independent random variables, drawn from a given distribution, it is clear following our introductory discussion that a central concern in our case is to determine the distribution of $``$ and hence that of $`Q`$. Moreover, to allow further progress, the form of $`Q`$ must be sufficiently simple (for example, local), and it does not seem obvious in advance whether this will be the case for any given problem. In fact, there is considerable flexibility, since $``$ determines $`Q`$ only up to left multiplication by an orthogonal matrix, and we have found convenient, explicit expressions for $`Q`$ in several one-dimensional examples . In addition, the approach proves useful even for problems in which $`Q`$ is known only implicitly, as we now show.

We consider the random field XY model, which has been studied extensively as a description of pinned charge density waves . It has the Hamiltonian

$$H=_0^L\left[\frac{1}{2}\mathrm{\Pi }^2+\frac{1}{2}(_x\varphi )^2+h(\varphi ,x)\right]𝑑x$$
(4)
where $`\mathrm{\Pi }`$ is the momentum conjugate to the angle $`\varphi `$, $`h(\varphi ,x)=h(x)\mathrm{cos}(\varphi (x)\chi (x))`$, and $`h(x)`$ and $`\chi (x)`$ are the random amplitude and phase of an applied field, with $`h(\varphi ,x)=0`$ and $`h(\varphi ,x)h(\varphi ^{},x^{})=\delta (xx^{})h_0\mathrm{cos}(\varphi \varphi ^{})`$. This continuum problem has, in the discrete notation used above, $`M=I_N`$, $`C=0`$ and hence $`\omega =\pm \sqrt{\kappa _i}`$. We are therefore concerned only with $`K`$: for the sake of an obvious analogy, we refer to it as a Hamiltonian, denoting it by $``$ and its eigenvalues by $`E`$ in place of $`\kappa `$.

Stationary configurations of the field $`\varphi `$ satisfy

$$_x^2\varphi +_\varphi h(\varphi ,x)=0$$
(5)
while amplitudes $`\psi `$ of normal mode excitations about the ground-state $`\varphi _0`$, with frequency $`\omega =\pm \sqrt{E}`$, obey

$$\psi _x^2\psi +_\varphi ^2h(\varphi _0,x)\psi =E\psi .$$
(6)
We want to find the density of states for the Hamiltonian of Eq. (6). In this form, this problem has been the subject of many publications . Summarising what is known: a characteristic energy scale, the pinning energy  $`E_p`$, separates two regimes. For $`EE_p`$ the density of states is only weakly affected by the random field and can be computed perturbatively in the field strength. For $`EE_p`$ the density of states is strongly influenced by disorder and believed to vary with $`E`$ as a power law, $`d(E)E^\beta `$, but the value of $`\beta `$ remains controversial. Here we show, in argeement with arguments due to Aleiner and Ruzin  and to Fogler , that $`\beta =3/2`$.

Since Eq. (6) has the form of a Schrödinger equation with a random potential $`_\varphi ^2h(\varphi _0,x)`$, it is tempting to anticipate behaviour familiar from simple choices of potential distribution, such as Gaussian white noise. That would be too naive, however, because – as for any bosonic Hamiltonian – the spectrum of $``$ is positive and so the potential cannot be arbitrarily random. To investigate the implications of this fact, it is useful to follow Feigelman, introducing the notion of a partial energy, $`(\varphi ,y)`$: the ground-state energy of the half-chain with coordinate values $`0xy`$ and boundary condition $`\varphi (y)=\varphi `$. Interpeting Eq. (5) as an equation of motion for a system with space coordinate $`\varphi `$ and time coordinate $`x`$, $`(\varphi ,y)`$ plays the role of an action and satisfies the Hamilton-Jacobi equation

$$_x+\frac{1}{2}\left(_\varphi \right)^2=h(\varphi ,x)$$
(7)
while the classical trajectory $`\varphi _0(x)`$ obeys $`d\varphi _0/dx=_\varphi `$. Now define $`V(x)=_\varphi ^2(\varphi _0(x),x)`$, which we call the chiral potential. As a direct consequence of Eq. (7) it satisfies

$$\frac{d}{dx}V+V^2=_\varphi ^2h(\varphi _0(x),x).$$
(8)
The Hamiltonian of Eq. (6) can therefore be rewritten as $`=Q^TQ`$, where

$$Q=\frac{d}{dx}+V(x).$$
(9)
Thus the bosonic problem of Eq. (6) is equivalent to a one-dimensional chiral problem specified by Eqs. (3) and (9), in which $`V(x)`$ should be determined along with $`\varphi _0(x)`$.

At this point we can draw on the work of Comtet et al, in which one-dimensional chiral problems of this kind have been analyzed in detail. In particular, they show that if the chiral potential $`V(x)`$ has a positive average $`V(x)=E_p^{1/2}`$ (sometimes referred to as staggering ), then the low-lying ($`EE_p`$) states of Eq. (3) are localized with a localization length $`\xi E_p^{1/2}`$. In this regime, the integrated density of states $`N(E)=_0^Ed(E^{})𝑑E^{}`$ can be calculated as the probability of a negative fluctuation of $`V(x)`$ in the interval $`x_1<x<x_2`$ for which $`2_{x_1}^{x_2}𝑑xV(x)2U<\mathrm{log}(E)`$. The probability for a rare event of this kind is expected under many circumstances to vary with $`E`$ as $`\mathrm{exp}[\alpha \mathrm{log}(E)]=E^\alpha `$, but the value of $`\alpha `$ depends on the details of the distribution of $`V(x)`$. In the following, we establish the equivalence implied by our notation, between $`V(x)^2`$ and pinning energy, and calculate the low-lying density of states for the random field XY model by studying rare negative fluctuations of the particular chiral potential that arises in this context.

To appreciate why $`V(x)>0`$, consider the ground state of the half-chain as a function of the boundary condition, $`\varphi `$: it varies smoothly except for jumps at a small number of isolated points within one period. At these points $`(\varphi ,y)`$ is continuous but has an unward cusp ; since $`(\varphi ,y)`$ is periodic, an average of $`V(x)_\varphi ^2`$ computed excluding these cusps (as on $`\varphi _0`$) is naturally positive, and its size can be estimated from the average sum of discontinuities in $`_\varphi `$. Alternatively, one can attempt a more detailed analysis of Eq. (7), which is similar to the KPZ equation , but with two differences. One, the absence of a diffusion term $`D_x^2`$, is unimportant since it is well known that the limit $`D0`$ in the KPZ equation yields solutions which are the global minimum of Eq. (4). The second, involving the nature of noise correlations is more significant: while studies of the KPZ equation deal with noise that is uncorrelated in both $`x`$ and $`\varphi `$, we are concerned with noise that has long-range correlations in $`\varphi `$. This situation is familiar in the context of Burgers turbulence , and indeed, introducing $`u=_\varphi `$ one finds the Burgers equation $`u_x+uu_\varphi =_\varphi h.`$ Cusps in $`(\varphi ,y)`$ as a function of $`\varphi `$ corrspond to the shock waves of Burgers turbulence, and using dimensional analysis the dependence on noise strength $`V(x)h_0^{1/3}`$ can be obtained from Eq. (7), as in the standard picture of the pinning energy.

Turning to fluctuations in $`V(x)`$, a direct attack, via a solution of Eq. (7), is not appropriate since many aspects of Burgers turbulence remain controversial. Instead we treat Eq. (8) as a Langevin equation for $`V(x)`$ with $`^2h(\varphi _0(x),x)`$ playing the role of a random force. If $`\varphi _0(x)`$ is determined by minimising Eq. (4) for fixed $`\varphi (0)`$ and $`\varphi (L)`$, the resulting random force has built-in correlations which ensure that $`V(x)`$ remains bounded. Alternatively, one may solve Eq. (5) for $`\varphi (x)`$ with fixed initial conditions, $`\varphi (0)`$ and $`_x\varphi (0)`$. In this case, the random force is uncorrelated and appropriate values for $`_x\varphi (0)`$ generate ground states corresponding to different choices of $`\varphi (L)`$. Other values of $`_x\varphi (0)`$, however, correspond to maxima of Eq. (4), and still others to higher energy local minima. From the perspective of the Langevin equation, a maximum of Eq. (4) is signalled by escape of $`V(x)`$ towards large negative values. Those noise realisations that result in escape should be eliminated by supplementing Eq. (8) with an absorbing boundary condition at $`V(x)=\mathrm{}`$. The trajectories of $`V(x)`$ that are not absorbed should be weighted in two ways. First, we must eliminate solutions to Eq. (5) which are only local minima of Eq. (4). They are separated from the absolute minimum by maxima and hence have neighboring trajectories on which $`V(x)`$ escapes. The necessary weight is therefore the probability that no near neighbors of a given trajectory are absorbed. Second, we weight trajectories representing the absolute minimum for different $`_x\varphi (0)`$ in such a way as to give a uniform density on $`\varphi (L)`$.

Omitting in the first instance both these weights, the form of the probability distribution $`P(U)`$ at large negative $`U`$ can be found by passing from the Langevin equation to a path integral, fixing $`U`$ with a Lagrange multiplier, and calculating the lowest eigenvalue of the corresponding Fokker-Planck operator, which is non-zero because of the absorbing boundary condition. In this way we obtain $`P(U)e^{2U}`$ as $`U\mathrm{}`$. To calculate the weighting factors, we consider a family of trajectories parameterized by $`s`$, with coordinate $`\varphi (x)+\eta (s,x)`$ and chiral potential $`V(x)+W(s,x)`$, where $`\eta (s,x)`$ and $`W(s,x)`$ are assumed small. Deriving linearized evolution equations by expanding Eqs. (5) and (8), and solving these, we find $`\eta (s,x_2)=e^U\eta (s,x_1)`$ and $`W(s,x_2)=e^{2U}F(s)`$, where

$`F(s)=`$ $`W(s,x_1)`$ (10)
$`+\eta (s,x_1)`$ $`{\displaystyle _{x_1}^{x_2}}\mathrm{exp}\left(3{\displaystyle _{x_1}^x}V(y)𝑑y\right)_\varphi ^3h(\varphi ,x)dx.`$ (11)
In order that no trajectories neighboring $`\varphi (x)`$ are absorbed, we require $`W(s,x_2)1`$ for all $`s`$ and hence, at large negative $`U`$, $`F(s)e^{2U}`$, which occurs with probability $`e^U`$. Of the surviving trajectories, we should retain only those that have a fluctuation in the integrated chiral potential similar to $`U`$. These satisfy the condition $`W(s,x_2)1`$, and hence $`F(s)e^{2U}`$, implying $`|\eta (s,x_1)|e^U`$ and $`|\eta (s,x_2)|e^{2U}`$, thus generating a weight $`e^{2U}`$. Combining these factors, the probability density that $`U`$ is large and negative on $`\varphi _0(x)`$ varies as $`e^{5U}`$, yielding $`N(E)E^{5/2}`$, and hence $`\rho (\omega )\omega ^4`$, as anticipated above. A detailed account of these arguments will be presented elesewhere .

The mapping described here of Eq. (4) into a chiral problem also goes through for a discrete version of the XY spin chain. It gives, not surprisingly, a staggered one-dimensional random hopping model, as studied for example in Refs. . This formulation is especially well suited to numerical calculations using the transfer matrix technique. Fig. 1 shows the result of a computation of $`N(E)`$ for a XY spin chain of $`10^6`$ sites, with Hamiltonian $`H=_i\left\{h_i(\varphi _i)\mathrm{cos}(\varphi _i\varphi _{i1})\right\}`$, where $`h_i(\varphi )=A_i\mathrm{cos}(\varphi \chi _i)`$, and $`[A\mathrm{cos}(\chi ),A\mathrm{sin}(\chi )]`$ is uniformly distributed on a disc of radius $`0.01`$. The power law obtained analyticially is confirmed.

We remark finally that, while the random field XY chain represents a special case in the sense noted above, since frequencies and stiffnesses are simply related, we expect the same form more generally for $`\rho (\omega )`$, as implied by the discussion preceeding Eq. (3). Calculations for the random field Heisenberg spin chain lead to two-channel chiral problems involving $`\stackrel{~}{}`$ and $`\stackrel{~}{}^{}`$. For such problems we obtain  $`d(\kappa )\kappa ^{3/2}`$ and $`\rho (\omega )\omega ^4`$.

We thank A. Altland and M. Zirnbauer for valuable discussions, and acknowledge support by EPSRC under grant GR/J78327.


format: text

text:
# Predictions for the frequency and orbital radii of massive extrasolar planets

## 1 Introduction

Radial velocity surveys of nearby stars show that a significant fraction – at least 8% – have massive planets with orbital radii substantially less than that of Jupiter (Marcy & Butler 2000; Udry et al. 2000; Butler et al. 2001). Most of these planets lie at radii where massive planet formation is theoretically believed to be difficult (Bodenheimer, Hubickyj & Lissauer 2000). The difficulties are most pronounced for those planets at the smallest orbital radii, where the temperatures in the protoplanetary disc would have exceeded those for which ices (and possibly even dust) can exist. These problems can be avoided if planets formed at larger radii, and then migrated inwards to where they are detected today. Mechanisms that could lead to this migration include gravitational interaction with a gaseous and viscous protoplanetary disc (Lin, Bodenheimer & Richardson 1996), planet-planet scattering (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997; Ford, Havlickova & Rasio 2001; Papaloizou & Terquem 2001), or scattering of planetesimals by massive planets (Murray et al. 1998).

In the right circumstances, all of the suggested mechanisms can lead to substantial migration. More quantitative comparisons with the data are therefore warranted. Recently, Trilling, Lunine & Benz (2002) have reported promising results for the protoplanetary disc migration model. They showed that the broad distribution of orbital radii was consistent with the planets becoming stranded, during their inward migration, by the dispersal of the protoplanetary disc. They used the model to estimate both the fraction of stars that must have formed planets, and how many ought to have survived migration.

In this paper, we extend the work of Trilling et al. (2002). We include a physical mechanism for disc dispersal into a model for the formation and migration of massive planets, and make quantitative comparison with the observed distribution of planetary orbital radii. Our disc model, described in Section 2, combines viscous evolution with mass loss at the outer edge, for example due to photoevaporation (Johnstone, Hollenbach & Bally 1998; Clarke, Gendrin & Sotomayor 2001). We then run this model repeatedly, on each occasion adding a single planet to the disc at a specified formation time $`t_{\mathrm{form}}`$. By varying $`t_{\mathrm{form}}`$, subject to the constraint that no planets can be formed once the disc mass is too low, we study in Section 3 how the final orbital radii of planets depend upon the formation time. With the addition of plausible assumptions about how the rate of planet formation varies with time, we then convert the results into a prediction for the radial distribution of planetary orbits. In Section 4 we compare this distribution with the current data, and find that good agreement is obtained. In Section 5 we make use of the model to estimate the fraction of stars targeted by radial velocity surveys that ought to harbour planets, including those with masses too low, or orbital radii too large, to be currently detected. We also predict the number of stars with planets at radii suitable for detection via microlensing, and the fraction of systems where significant planetary migration does not occur. Our conclusions are summarized in Section 6.

## 2 Methods

### 2.1 Fully viscous disc model

The simplest protoplanetary disc model that meets the basic observational constraints on disc mass, disc lifetime, and disc destruction time-scale is a viscous disc with mass loss from the outer regions (Hollenbach, Yorke & Johnstone 2000; Clarke, Gendrin & Sotomayor 2001). We use this model for the majority of the calculations described in this paper, since it has a minimal number of free parameters.

For the viscous disc model runs, we assume that over the relevant range of radii (roughly, between 0.1 AU and 10 AU), the kinematic viscosity $`\nu `$ can be described as a power law in radius $`R`$,

$$\nu =\nu _0\left(\frac{R}{1\mathrm{AU}}\right)^\beta .$$
(1)
The initial conditions for the surface density are a constant accretion rate disc, with a profile which corresponds to the assumed viscosity of equation (1). At radii $`RR_{\mathrm{in}}`$, the inner edge of the disc, this implies,

$$\mathrm{\Sigma }=\frac{\dot{M}}{3\pi \nu }=\mathrm{\Sigma }_0R^\beta $$
(2)
where $`\dot{M}`$ is the accretion rate and $`\mathrm{\Sigma }_0`$ a constant. For our standard disc model, we take $`\beta =3/2`$ (Weidenschilling 1977; Hayashi 1981), and choose $`\nu _0=1.75\times 10^{13}\mathrm{cm}^2\mathrm{s}^1`$ to give a sensible time-scale for disc evolution of a few Myr. Finally, $`\mathrm{\Sigma }_0`$ is set in the initial conditions so that the initial disc mass is $`0.1M_{}`$. The inner and outer boundaries for the calculation are $`R_{\mathrm{in}}=0.066\mathrm{AU}`$, and $`R_{\mathrm{out}}=33.3\mathrm{AU}`$. To extend the predictions to smaller radii, a limited number of runs have also been completed with an inner boundary at $`R_{\mathrm{in}}=0.03\mathrm{AU}`$. A zero torque boundary condition (i.e. a surface density $`\mathrm{\Sigma }=0`$) is applied at the inner boundary, while at the outer boundary we set the radial velocity to zero.

Current observations do not directly constrain the surface density profile of protoplanetary discs at the AU scale. To check how sensitive our results are to changes in this profile, we have also run a model with $`\beta =1/2`$, again normalized to give a sensible time-scale for disc evolution. A flatter profile than the $`R^{3/2}`$ of our standard model is expected if the angular momentum transport within the disc follows the Shakura & Sunyaev (1973) $`\alpha `$-prescription, with a constant $`\alpha `$ (e.g. Bell et al. 1997).

Even a small mass of gas, trapped exterior to the orbit of a massive planet, will eventually soak up a large amount of angular momentum and drive migration. This gas has to be lost in order to obtain a converged final radius for migrating planets. Trilling et al. (2002) make the simplest assumption, and halt migration by removing the disc instantaneously. We adopt a different approach, and assume that the gas is removed by photoevaporation. In this process, the surfaces of the disc are heated by the absorption of ultraviolet radiation, which may originate either from the central star or from an external source (e.g. a massive star in a young cluster). The hot gas can escape as a wind at radii $`R>R_g`$, where the escape velocity is less than the sound speed of the heated material. Observations of the Orion nebula show that this process drives mass loss at significant rates, at least in clusters containing massive stars (Bally, O‘Dell & McCaughrean 2000).

Our implementation of photoevaporation within a time-dependent disc model is based upon the description given by Johnstone, Hollenbach & Bally (1998; see also Shu, Johnstone & Hollenbach 1993; Richling & Yorke 2000). The standard model has $`R_g=5\mathrm{AU}`$, and a radial distribution of mass loss,

$`\dot{\mathrm{\Sigma }}`$ $``$ $`R^1R>R_g`$
$`\dot{\mathrm{\Sigma }}`$ $`=`$ $`0R<R_g.`$ (3)
For a disc with an outer radius $`RR_g`$, this means that the total mass loss scales linearly with $`R`$. The initial mass loss rate (integrated over the disc out to 25 AU) is chosen to be $`\dot{M}=5\times 10^9M_{}\mathrm{yr}^1`$.

Figure 1 shows how the disc mass changes with time, for the standard model with $`\beta =3/2`$ and mass loss from $`R>5\mathrm{A}\mathrm{U}`$. The disc mass declines from $`0.1M_{}`$ to $`10^2M_{}`$ in 3 Myr, mainly due to accretion, before dropping rapidly once the effects of the disc wind take hold. Although the power laws and location of the turnover vary with the adopted parameters, generically similar behavior occurs for a range of disc viscosity and mass loss rates (Clarke, Gendrin & Sotomayor 2001).

### 2.2 Magnetically layered disc model

In order to check the robustness of the results, we investigate how the results would change if we adopted a fundamentally different, but strongly theoretically motivated, model for disc evolution. The magnetically layered model proposed by Gammie (1996) is based upon the observation that cold gas at $`T\stackrel{<}{}10^3\mathrm{K}`$ has a low ionization fraction, which suppresses angular momentum transport via magnetohydrodynamic turbulence (Fleming, Stone & Hawley 2000). If there are no other sources of angular momentum transport, then only the hot inner disc, plus a surface layer ionized by cosmic rays, will be viscous. At radii $`R1\mathrm{A}\mathrm{U}`$, this implies that accretion occurs only through a thin surface layer, below which lies a thick ‘dead zone’ of quiescent gas. Numerical models for the long-term disc evolution (Armitage, Livio & Pringle 2001) find that there is an early phase, lasting for perhaps a Myr, during which accretion occurs in short high accretion rate outbursts. This is followed by a quiescent phase, characterised by low rates of accretion onto the star. This is due to a much reduced (compared to fully viscous disc models) viscosity at radii of the order of 1 AU.

Detailed models for the layered disc are able to provide a reasonable fit to observations of the accretion rate in Classical T Tauri stars (Gammie 1996; but see also Stepinski 1998). For the purpose of studying planetary migration, we adopt an approximate approach, and assume that the most important effect of the layered disc is to reduce the viscosity near the midplane. Specifically, after 1 Myr, we reduce the viscosity where the surface density exceeds the thickness of the magnetically active surface layer. The new viscosity $`\nu _{\mathrm{layer}}`$ is,

$$\nu _{\mathrm{layer}}=\nu \times \mathrm{min}(\frac{100\mathrm{g}\mathrm{c}\mathrm{m}^2}{\mathrm{\Sigma }},1)$$
(4)
where $`\nu `$ is as given in equation (1), and the same viscosity parameters as for the standard disc model are chosen. To obtain a sensible disc lifetime we also increase the mass loss rate to $`\dot{M}=3\times 10^8M_{}\mathrm{yr}^1`$, and allow mass to be lost at all radii by taking $`R_g=R_{\mathrm{in}}`$.

### 2.3 Planet formation assumptions

The location in the disc where massive planet formation is easiest depends upon two competing factors. The characteristic time-scales for planet building in the disc are faster at small radii. However, most of the mass in protoplanetary discs lies at large radii, with a jump in the surface density of solid material beyond the ‘snow line’ where ices can form (Hayashi 1981; Sasselov & Lecar 2000). We assume that the outcome of this competition is that planets form at an orbital radius $`a=5\mathrm{A}\mathrm{U}`$, similar to that of Jupiter in the Solar System. The exact choice of this location is less important than the fact that it lies outside the radii where extrasolar planets are currently observed. In the model which we are testing, this means that migration is necessary to explain the orbital radii of all currently known planets.

In each run of the disc model, we form a planet with mass $`M_p`$, at time $`t_{\mathrm{form}}`$ (where $`t=0`$ is defined as the time when the disc mass is $`0.1M_{}`$). We specify the mass of our planets (typically $`2M_J`$, where $`M_J`$ is the mass of Jupiter) in advance of each run, but test at the time of formation that there is enough mass in the vicinity of the planet to form it consistently. Specifically, we check that the disc mass between $`R=0.6a`$ and $`R=1.6a`$ exceeds the desired planet mass. Provided that this is satisfied, we remove the appropriate amount of gas from the disc around $`R=a`$, and add a planet with mass $`M_p`$ at that location. This instantaneous formation scheme ignores the possibly lengthy time required to assemble the planet’s gaseous envelope (e.g. Lissauer 1993, and references therein). However, Trilling, Lunine & Benz (2002) have shown that including a gradual build-up of mass makes little difference to the outcome of migration.

After running the model a number of times, the basic output is a plot of the final radius of a planet as a function of the time when it formed. To convert this to the number of planets expected at different radii, which is the observable quantity, we also need to specify how the probability of massive planet formation, per unit time, varies with time. We assume that this probability is uniform (or equivalently, that the rate of massive planet formation, averaged over many stars, is constant). This cannot be true over long periods of time. However, we are only interested in the window of formation times which allow a massive planet to survive migration. This window, shown as the shaded region in Figure 1, is short compared to the disc lifetime, and as a result there is only a relatively small change in the physical properties of the disc during this time. For example the disc mass, which drops by two orders of magnitude over the first 3.5 Myr, changes by less than 50%. Whatever factors influence the probability of planet formation, it should therefore be a reasonable first approximation to assume that the probability is constant over the short range of formation times which lead to a surviving massive planet.

### 2.4 Model for planetary migration

Our model for planetary migration is almost identical to that of Trilling et al. (2002), who describe in detail the approximations involved. The coupled evolution of the planet and the protoplanetary disc is described by the equation (Lin & Papaloizou 1986),

$$\frac{\mathrm{\Sigma }}{t}=\frac{1}{R}\frac{}{R}\left[3R^{1/2}\frac{}{R}\left(\nu \mathrm{\Sigma }R^{1/2}\right)\frac{2\mathrm{\Lambda }\mathrm{\Sigma }R^{3/2}}{(GM_{})^{1/2}}\right].$$
(5)
The first term on the right-hand side of this equation describes the diffusive evolution of the surface density due to internal viscous torques (e.g. Pringle 1981). The second term describes how the disc responds to the planetary torque, $`\mathrm{\Lambda }(R,a)`$, where this function is the rate of angular momentum transfer per unit mass from the planet to the disc. We take,

$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{q^2GM_{}}{2R}}\left({\displaystyle \frac{R}{\mathrm{\Delta }_p}}\right)^4R<a`$
$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{q^2GM_{}}{2R}}\left({\displaystyle \frac{a}{\mathrm{\Delta }_p}}\right)^4R>a`$ (6)
where $`q=M_p/M_{}`$,

$$\mathrm{\Delta }_p=\mathrm{max}(H,|Ra|),$$
(7)
and $`H`$ is the scale height of the disc. Guided by detailed protoplanetary disc models (Bell et al. 1997), we adopt $`H=0.05R`$. This form for $`\mathrm{\Lambda }`$ (equation 6) is that used by Lin & Papaloizou (1986), modified to give a symmetric treatment for the disc inside and outside the orbit of the planet.

The transfer of angular momentum leads to orbital migration of the planet at a rate,

$$\frac{\mathrm{d}a}{\mathrm{d}t}=\left(\frac{a}{GM_{}}\right)^{1/2}\left(\frac{4\pi }{M_p}\right)_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}R\mathrm{\Lambda }\mathrm{\Sigma }dR.$$
(8)
Although simple, this formalism for treating the planet-disc interaction has been shown (Trilling et al. 1998) to give comparable results to more sophisticated methods (Takeuchi, Miyama & Lin 1996).

We solve equation (5) using an explicit method on a grid that is uniform in a scaled variable $`R^{1/2}`$, with 300 grid points between $`R_{\mathrm{in}}`$ and $`R_{\mathrm{out}}`$. The stellar mass is $`M_{}=M_{}`$. Typically, the timestep is limited by the radial velocity in the disc gas very near the location of the planet. Evolving the system on this (small) timestep serves no useful purpose once a gap has been opened, so we reduce $`\mathrm{\Lambda }`$, and thereby limit $`|v_r|`$, in the vicinity of the planet. The results of the migration calculations are unaffected by this modification.

## 3 Results

### 3.1 Predicted semi-major axis of stranded planets

Figure 2 shows the final semi-major axis of planets as a function of the formation time, using the standard disc model of Fig. 1. Each point represents a separate run of the model. The fate of a planet depends upon when it was formed, and upon its mass. Planets that are formed too early migrate inward of 0.1 AU, the smallest radius we consider quantitatively in this study. What happens to them subsequently depends upon whether there is a stopping mechanism to halt further migration, but many of them seem likely to be consumed by the star. Planets that start to form too late can migrate outwards, because at late times the gas in the disc at 5 AU is itself moving away from the star. At still later times, it becomes impossible to form a Jupiter mass planet at all, because the dispersing disc has too little surface density (c.f. the discussion in Shu, Johnstone & Hollenbach 1993). In between these limits, there is a window of formation times, lasting for around 0.5 Myr for a $`2M_J`$ planet, which result in the planet being stranded at radii between 0.1 AU and 6 AU. The migration of more massive planets slows down earlier, since their larger angular momentum produces greater resistance to disc-induced migration (Syer & Clarke 1995). Hence, the window of allowable formation times is longer for more massive planets than for less massive ones.

### 3.2 Dependence upon disc model and planet mass

Figure 2 shows how the final planetary radius depends upon the formation time (note that we plot $`a(t_{\mathrm{form}})`$ on a logarithmic scale in this figure). The curve steepens towards small radii / early formation times. This implies that if planet formation occurs with uniform probability per unit time during the allowed window, then fewer planets per logarithmic interval in $`a`$ are expected at small radius (i.e. equal intervals in $`\mathrm{log}(a)`$ correspond to smaller intervals in $`\mathrm{\Delta }t_{\mathrm{form}}`$ at small $`a`$).

With the assumption of a uniform rate of planet formation, the predicted number of planets per logarithmic interval in radius is,

$$\frac{\mathrm{d}N_p}{\mathrm{d}\mathrm{log}(a)}\left(\frac{\mathrm{d}\mathrm{log}(a)}{\mathrm{d}t_{\mathrm{form}}}\right)^1,$$
(9)
which can readily be derived from Fig. 2. Figure 3 shows the results for two different planet masses, and for two variants on the standard disc model. A cubic spline fit has been used to obtain a smooth estimate of the derivative in equation (9), and the curves have been normalised to unity at 1 AU.

The predicted number of planets per logarithmic interval in $`a`$ rises rapidly with increasing radius for all the models considered. For the standard disc model, the predicted number rises by more than a factor of two between 0.1 AU and 1 AU, and by a further factor of two by 3 AU. The good agreement between the curves for $`1M_J`$ and $`2M_J`$ planets shows that the observed mass function for massive planets is predicted to be the same at different radii. This means that we can use the observed mass function at small radii, which is complete down to lower masses, to estimate how many low mass planets must accompany known high mass planets at larger radii. It does not mean that the initial mass function for planets is the same as the observed one. In fact, Fig. 2 makes it clear that high mass planets can form over a longer interval than low mass planets and still survive migration. To end up with the same number of $`1M_J`$ as $`2M_J`$ planets per unit interval in $`\mathrm{log}(a)`$, more low mass planets must have formed. In other words, the (currently unobservable) initial mass function for planets must rise more steeply towards low masses than the observed one.

Figure 3 also shows the predicted relative abundance of planets for two variants of the disc model. Increasing the inner radius beyond which mass is lost to 10 AU does not alter the predicted distribution of planets within 5 AU. Likewise, the predicted distribution is not substantially changed in a model where the power-law slope of the viscosity is $`\beta =1/2`$. Note that for the $`\beta =1/2`$ run we considered planets with mass $`0.5M_J`$, since insufficient gas was available between $`0.6a`$ and $`1.6a`$ to form planets of several Jupiter masses.

### 3.3 Comparison with an instantaneous disc dispersal model

In Figure 4, we show how the model including disc mass loss compares to a simpler model in which the disc is assumed to be lost instantaneously at a time $`t_{\mathrm{disperse}}`$. The disc model used is identical to our standard model with $`\beta =3/2`$, except that the rate of mass loss is set to zero. We compute two models, with $`t_{\mathrm{disperse}}`$ equal to 3 Myr and 4 Myr respectively.

From Fig. 4, it can be seen that the instantaneous disc dispersal model makes very similar predictions for the radial distribution of inwardly migrating planets. Larger differences, however, occur for outward migration. In the instantaneous disc dispersal model, there is no outward migration of planets, because the radial velocity in the disc at 5 AU remains inward at $`t_{\mathrm{disperse}}`$. The model including mass loss, conversely, allows a substantial fraction of planets to migrate outwards. From Fig. 2, we find that for Jupiter mass planets as many as a third may end up in orbits exterior to that in which they formed.

### 3.4 Effect of ongoing planetary accretion

Even after a gap has been opened, numerical simulations show that planetary accretion may continue (Bryden at al. 1999; Kley 1999; Lubow, Seibert, & Artymowicz 1999; Nelson et al. 2000). Significant ongoing accretion – which is not included in our models – would lead to a correlation between planet mass and orbital radius, with more massive planets expected at small orbital radius. There is no evidence for such a correlation in the data, which could either indicate that migration has not occurred, or that ongoing accretion occurred only at a low level.

Lubow, Seibert, & Artymowicz (1999) find that, for low masses, the rate of ongoing accretion can be comparable to the disc accretion rate outside the planet. This is of the order of $`10^9M_{}\mathrm{yr}^1`$ at the epoch when surviving planets form. Migration from 5 AU takes a few hundred thousand years, so we estimate that accretion across the gap could add a few tenths of a Jupiter mass to low mass planets. This is a relatively small fraction of the mass for the $`2M_J`$ planets that are the main focus of this paper. However, it is only marginally consistent with the existence of the lowest mass hot Jupiters, which have $`M_p\mathrm{sin}(i)<0.3M_J`$. This may indicate that these planets formed within the snow line (leading to a shorter migration time).

### 3.5 Magnetically layered protoplanetary discs

Figure 5 shows the results for the runs using the magnetically layered disc model. The decline of the disc mass with time shows the effect of the switch to a low viscosity state. After the initial fully viscous phase, which resembles the standard model, there is a long plateau which arises because the low viscosity in the inner disc produces a bottleneck to accretion onto the star. During this period the disc is slowly destroyed from the outside inwards by the disc wind, which finally disperses the disc completely after 6 Myr. We note that, in a magnetically layered disc model, high rates of mass loss from the disc are needed in order to disperse the disc within a reasonable time-scale. The mass loss rates inferred for discs in parts of the Orion nebula are easily large enough (Johnstone, Hollenbach & Bally 1998), but it is unknown whether high enough rates are possible for isolated discs.

Figure 5 also shows the final location of migrating planets. As before, we assume these planets all formed at 5 AU, and vary the time of formation in a series of runs. The resulting $`a(t_{\mathrm{form}})`$ curve flattens for $`a\stackrel{>}{}1\mathrm{AU}`$. As with the standard disc model, therefore, we predict a broad distribution of final orbital radii, with more planets (per logarithmic interval in $`a`$) ending up at large radii than at small radii. However, there are also important differences. Planets must form earlier (relative to the disc lifetime) in the layered disc if they are to migrate successfully to small orbital radii. This is because little migration can take place once the disc has entered its quiescent, low viscosity phase. For a planet to end up inside 1 AU, it must have formed at an earlier time when the disc was still viscous and able to drive rapid migration.

## 4 Comparison with current data

Before comparing our results to current data, we need to consider possible biases. In an idealized radial velocity search for extrasolar planets, two important selection effects are expected. First, there is almost no sensitivity to planets with orbital periods $`P`$ longer than the duration of the survey (until a whole orbit is seen, the planetary signal may be confused with that from a more massive, distant companion). Second, a planet with mass $`M_p`$ and orbit inclination $`i`$ produces a radial velocity signal with an amplitude proportional to $`M_p\mathrm{sin}(i)a^{1/2}`$. A survey which makes a fixed number of observations, with some given sensitivity, can therefore detect planets above a minimum $`M_p\mathrm{sin}(i)a^{1/2}`$. Although this is a relatively weak radial dependence, it still corresponds to an order of magnitude difference in the minimum detectable mass over the range of radii probed by current surveys.

Real surveys are less easily characterised. Instruments have improved steadily over time, and the intrinsic limits to radial velocity measurements vary on a star-by-star basis. Nevertheless, the selection effect caused by the scaling of the minimum detectable mass with $`a^{1/2}`$ can readily be seen in the data. Figure 6 shows the location of all the known extrasolar planets with $`M_p\mathrm{sin}(i)<10M_J`$, in the $`M_p\mathrm{sin}(i)`$ \- $`a`$ plane. Most of the hot Jupiters at $`a<0.1\mathrm{AU}`$, where surveys are most sensitive, have masses that are actually substantially smaller than that of Jupiter, while the typical detected planet with an orbital radius greater than 1 AU has a mass of several $`M_J`$.

A simple count of the total number of detected planets as a function of radius leads to the cumulative distribution shown in Figure 7. A Kolmogorov-Smirnov test (e.g. Press et al. 1989) shows that this is marginally consistent (KS probability $`P=0.05`$) with a uniform distribution in $`\mathrm{log}(a)`$. However, from Fig. 6, it can be seen that this uncorrected distribution includes low mass planets, at small orbital radius, that would not be detectable beyond 1 AU. The mass function of planets has substantial numbers of low mass planets (Zucker & Mazeh 2001; Tabachnik & Tremaine 2002), so this leads to a significant selection effect (Lineweaver & Grether 2002).

An unbiased estimate of the radial distribution of extrasolar planets can be obtained by considering a subsample of planets which are sufficiently massive that they could be detected at any radius. A conservative definition of such a subsample can be obtained by noting that for the Lick planet search, circa 1999, the minimum detectable mass at 99% confidence is quoted by Cumming, Marcy & Butler (1999) as,

$$M_p\mathrm{sin}(i)\stackrel{>}{}0.7M_J(a/AU)^{1/2}.$$
(10)
Taking this as a guide, a mass cut at $`M_{\mathrm{cut}}=M_p\mathrm{sin}(i)>1.2M_J`$ should yield an approximately complete subsample out to 3 AU. Since the precision of the observations is likely to have improved since 1999, we also consider a subsample with the mass cut placed at $`0.6M_J`$. This corresponds to a radial velocity amplitude at 3 AU of $`10.5\mathrm{ms}^1`$, which is equal to the smallest amplitude signal of any detected planet. As one might expect, this practical detection threshold lies at $`34\sigma `$, where $`\sigma 3\mathrm{ms}^1`$ is the best-case error on a single radial velocity measurement.

As shown in Fig. 7 and Table 1, neither of these subsamples is consistent with the hypothesis that it is drawn from a uniform distribution in $`\mathrm{log}(a)`$, which can be excluded at the 99.9% confidence level. This conclusion is also robust to different choices of the minimum and maximum orbital radii of planets in the subsample.

To compare with the theoretical distribution derived earlier, we henceforth take $`M_{\mathrm{cut}}=0.6M_J`$. This may somewhat underestimate the number of planets at radii of 2-3 AU, but any such bias is smaller than the uncertainties introduced by the relatively small sample size. We then take the current compilation of 76 detected planets (from exoplanets.org), discard those with masses outside the range $`0.6M_J<M_p\mathrm{sin}(i)<10M_J`$, and bin 54 survivors into four logarithmic intervals in semi-major axis, 0.033 AU – 0.1 AU, 0.1 AU – 0.3 AU, 0.3 AU – 0.9 AU, and 0.9 AU – 2.7 AU.

Figure 8 shows the comparison between observation and theory. The disc inside 0.1 AU may well be strongly influenced by the protostellar magnetosphere (Konigl 1991; Najita, Edwards, Basri & Carr 2000), which we have chosen not to model. We therefore concentrate attention on the outer three radial bins, between 0.1 AU and 2.7 AU. The number of known extrasolar planets in this (approximately) complete subsample rises rapidly with increasing semi-major axis. Although the error bars on the data are large, excellent agreement is obtained with the predictions of the theoretical migration model using the standard disc model. This agreement persists when different choices of binning and sample selection are made. In particular, we can obtain an equally good fit (with a different normalization) if we adopt a more conservative mass cut at $`1.2M_J`$.

From the Figure, we also note that in this mass range, the number of hot Jupiters at $`a<0.1\mathrm{AU}`$ is in fact consistent with an extrapolation of the data from larger radii. There is no need for an additional stopping mechanism to explain the abundance of these planets at the smallest radii. There is, however, weak evidence for an over-abundance of lower mass planets (with $`M_p\mathrm{sin}(i)<0.6M_J`$, not plotted in the Figure) in the innermost radial bin.

## 5 Discussion

### 5.1 Frequency of massive planets

The good agreement between theory and observations seen in Figure 8 allows us to use the model to estimate the fraction of stars targeted by radial velocity surveys that possess massive planets, including those with masses too low, or semi-major axis too large, to have been detected so far. Of course, we cannot extrapolate to masses smaller than the $`0.2M_J\mathrm{sin}(i)`$ that is the smallest yet detected, as we have no information at all on the mass function below this threshold. We can, however, estimate how many planets with $`0.5M_J\mathrm{sin}(i)`$ ought to exist at 1 AU (say), beyond the radius where they are currently detectable.

To do so, we note that the radial velocity surveys are most complete (i.e. reach down to the smallest masses) at small radii. Between 0.1 AU and 0.2 AU, equation (10) suggests that planets with masses of $`0.3M_J\mathrm{sin}(i)`$ and above are detectable. 8 planets are known in this radial interval, corresponding to approximately 1% of target stars (since the total haul of 76 planets implies that approximately 8% of stars have planets). The range of planet formation times which result in the planet becoming stranded between 0.1 AU and 0.2 AU is approximately $`2.5\times 10^4\mathrm{yr}`$, which may be compared to the 0.42 Myr window which leaves planets between 0.1 AU and 5 AU. Hence, we estimate that the fraction of target stars possessing planets with masses $`0.3M_\mathrm{J}<M_p\mathrm{sin}(i)<10M_\mathrm{J}`$, and radii $`0.1\mathrm{AU}<a<5\mathrm{A}\mathrm{U}`$, is $``$ 15 %. Roughly half of these have already been detected.

More speculatively, we can also try to extrapolate to larger radii. If we assume that massive planets at 5-8 AU formed near the inner boundary of that region and migrated outwards, then a similar analysis to that above suggests that an additional 5-10% of stars could possess such planets. This conclusion is more tentative than for inward migration, however, because the extent of outward migration varies significantly with planet mass (Fig. 2), and could be very different in the magnetically layered disc model.

### 5.2 Frequency of massive planet formation

In our standard model, planets which survive migration with final orbital radii between 0.1 AU and 6 AU must have formed during a 0.5 Myr allowed window. This window lasts for only $``$ 20% of the disc lifetime. Migration without a stopping mechanism is thus a moderately inefficient process (Trilling et al. 2002) – several planets are likely to have formed for every planet which survives today. This inefficiency means that either a large fraction of stars formed of the order of one planet, or that a smaller fraction of stars formed many planets. Our models do not distinguish between these possibilities, though in principle they make different predictions. For example, the fraction of stars which consume planets (and may show resulting metallicity enhancements) is smaller if only a few stars form many massive planets. The interpretation of current metallicity measurements within such a framework, however, remains somewhat controversial (e.g. Gonzalez 1997; Udry, Mayor & Queloz 2001; Murray et al. 2001; Pinsonneault, DePoy & Coffee 2001; Suchkov & Schultz 2001).

### 5.3 Number of planets detectable via microlensing

Detailed monitoring of microlensing events has so far failed to find anomalies in the light curves characteristic of massive, bound planets (Bennett & Rhie 1996; Albrow et al. 2001; Gaudi et al. 2002). The results exclude the possibility that more than around a third of lensing stars have Jupiter mass planets at radii $`1.5\mathrm{AU}<a<4\mathrm{A}\mathrm{U}`$. Our results are consistent with this limit. For masses greater than $`0.3M_J`$, we estimate that 7% of the stars targeted by radial velocity surveys have planets within this range of radii. If this is also true of the lensing stars (which are typically less massive – around $`0.3M_{}`$), then it suggests that the existing monitoring campaigns are within a factor of $`5`$ of reaching limits where a detection may be expected.

### 5.4 Probability of forming a planet which does not migrate

Unlike in the case of extrasolar planetary systems, the evidence for migration in the Solar System is extremely limited. In situ analysis of the Jovian atmosphere shows that it may contain material that originated in the outer Solar System (Owen et al. 1999), while substantial migration of Uranus and Neptune has been suggested (Thommes, Duncan & Levison 1999). There is little to suggest that Jupiter formed a significant distance away from its current location.

Migration models permit both inward and outward migration, so forming a planet which does not migrate is perfectly possible, though it requires fortuitous timing. If we take ‘no significant migration’ to mean that the planet moves a distance,

$$\frac{\mathrm{\Delta }a}{a}0.1,$$
(11)
then in our standard disc model the window of formation time over which this occurs lasts for $`\mathrm{\Delta }t=6.8\times 10^4\mathrm{yr}`$. This may be compared with the 0.5 Myr range of formation times that result in the planet having a semi-major axis in the range $`0.1\mathrm{AU}<a<6\mathrm{A}\mathrm{U}`$. Making the same assumption as before – that the rate of planet formation remains constant across this larger interval – we estimate that no significant migration occurs in approximately 10-15% of systems in which a planet survives. The no migration outcome which may describe Jupiter in the Solar System would then be an unlikely event, but not a rare one (Lineweaver & Grether 2002). Of course forming several planets, none of which migrated significantly, would be much less likely. However, as already noted, migration of the other giant planets in the Solar System has been seriously considered and does not appear to be excluded on observational grounds.

### 5.5 Eccentricity

Our model for migration ignores the typically substantial eccentricities of extrasolar planets. How serious this omission is depends upon the (unknown) mechanism which leads to the eccentricity. If interactions between a single planet and the disc are responsible (Artymowicz et al. 1991; Papaloizou, Nelson & Masset 2001; Goldreich & Sari 2002), then a description similar to the one we have developed here ought to suffice. Conversely, eccentricity may arise from planet-planet scattering in a multiple planet system (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997; Papaloizou & Terquem 2001). This would be accompanied by order unity changes in the semi-major axis of the surviving planet. The predicted distribution of planetary orbital radii, shown in Figure 8, would be changed if all planets (say) suffered a shrinkage of their orbits by a factor of two subsequent to migration. We have no additional need to invoke such a process to explain the orbital radii of the planets, but more work in this area is warranted.

## 6 Conclusions

We have shown that the radial distribution of massive extrasolar planets, beyond 0.1 AU, is consistent with the planets forming at 5 AU, before migrating inwards through a viscous protoplanetary disc. This means that although theory suggests that massive planets could form closer to the star (Papaloizou & Terquem 1999; Bodenheimer, Hubickyj & Lissauer 2000; Sasselov & Lecar 2000), there is no observational requirement for planets to form at $`a<3\mathrm{A}\mathrm{U}`$. The predicted distribution is a robust feature of disc models which include dispersal of the gas via a disc wind, so the level of comparison is currently limited by the small numbers of detected planets. The continuation of existing radial velocity surveys, along with forthcoming astrometric searches (e.g. Lattanzi et al. 2000), will allow for more stringent tests, provided that the selection effects of the surveys are well understood.

The good agreement between observations and theory allows us to use the model to estimate the fraction of stars that possess massive planets, including those currently undetected on account of their low mass or large semi-major axis. Limiting ourselves to planets with masses $`0.3M_\mathrm{J}<M_p\mathrm{sin}(i)<10M_\mathrm{J}`$, and radii $`0.1\mathrm{AU}<a<5\mathrm{A}\mathrm{U}`$, we estimate that around 15% of stars in the current target sample of extrasolar planet searches possess planets. In some of our models a significant number of planets migrate outwards, and these could form a sizeable additional population. These numbers are comparable to those found by Trilling, Lunine & Benz (2002), and imply that even in existing samples, most of the massive planets remain to be discovered.

Finally, we have used the results to quantify the fraction of systems in which the most massive planet does not suffer significant radial migration. This allows us to address the concern that migration theories could be inconsistent with the detailed knowledge of our own Solar System, in which Jupiter may have formed close to its current location. We find that a no-migration outcome requires fortuitous timing of the epoch of planet formation, but not extraordinary luck. Of the order of 10-15% of planetary systems are expected to have a massive planet which has migrated by less than 10% in orbital radius. Having a most massive planet at the orbital radius of Jupiter is thus expected to occur in around 1-2% of solar-type stars. The orbital radius of Jupiter, at least, is consistent with our understanding of extrasolar planets, though the near-circular orbit may yet prove to be unusual.

## Acknowledgments

PJA thanks the Institute of Astronomy for their hospitality during the course of this work. JEP is grateful to STScI for continuing support under their Visitors Program. SL acknowledges support from NASA grants NAG5-4310 and NAG5-10732. We thank the anonymous referee for a detailed and helpful review.


text:
# Generalized Wigner Function Formulation for Quantum Systems with Open Boundaries

## Abstract

A rigorous microscopic theory for the description of quantum-transport phenomena in systems with open boundaries is proposed. We shall show that the application of the conventional Wigner-function formalism to this problem leads to unphysical results, like injection of coherent electronic states from the contacts. To overcome such basic limitation, we propose a generalization of the standard Wigner-function formulation, able to properly describe the incoherent nature of carrier injection at the device spatial boundaries as well as the interplay between phase coherence and energy relaxation/dephasing within the device active region. The proposed theoretical scheme constitutes the first rigorous derivation of the phenomenological injection model commonly employed in the simulation of open quantum devices.

Present-day technology pushes device dimensions toward limits where the traditional semiclassical transport theory MC  can no longer be applied, and more rigorous quantum-transport approaches are required QT . However, in spite of the quantum-mechanical nature of carrier dynamics in the core region of typical nanostructured devices —like semiconductor superlattices and double-barrier structures— the overall behaviour of such quantum systems is often the result of a complex interplay between phase coherence and energy relaxation/dephasing QCL , the latter being primarily due to the presence of spatial boundaries WF . It follows that a proper treatment of such novel nanoscale devices requires a theoretical modelling able to properly account for both coherent and incoherent —i.e., phase-breaking— processes on the same footing. To this end, a generalization to open systems —i.e., systems with open boundaries— of the well known Semiconductor Bloch Equations (SBE) SBE  has been recently proposed PRL . However, the theoretical analysis presented in PRL  is primarily related to the interplay between phase coherence and energy relaxation within the device active region, and —apart from its abstract formulation— no detailed investigation of the carrier-injection process (from the electrical contacts into the device active region) has been performed so far.

Aim of the present Letter is to provide a rigorous quantum-mechanical description of the coupling dynamics between the device active region and external charge reservoirs, able to account for the semi-phenomenological injection models commonly employed in state-of-the-art simulations of realistic one- and two-dimensional open quantum devices IWCE . Among such simulation strategies it is worth mentioning the approach recently proposed by Fischetti and co-workers Max ; note-Max : By denoting with $`f_\alpha `$ the carrier distribution over the electronic states $`\alpha `$ of the device, the transport equation proposed in Max  is of the form:

$$\frac{d}{dt}f_\alpha =\underset{\alpha ^{}}{}\left(W_{\alpha \alpha ^{}}f_\alpha ^{}W_{\alpha ^{}\alpha }f_\alpha \right)+\frac{f_\alpha ^bf_\alpha }{\tau _\alpha }.$$
(1)
Here, $`f_\alpha ^b`$ denotes the equilibrium carrier distribution in the contacts while $`\tau _\alpha `$ can be regarded as the device transit time for an electron in state $`\alpha `$. As anticipated, in spite of a rigorous treatment of the scattering dynamics, the last (relaxation-time-like) term describes carrier injection/loss on a partially phenomenological level and does not depend on the real position of the device spatial boundaries.

In order to provide a fully microscopic real-space description of the carrier-injection process, we shall start rivisiting the theoretical approach proposed in PRL . By applying to the conventional SBE for a closed system PRL ; note-rho ; note-L

$$\frac{d}{dt}\rho _{\alpha _1\alpha _2}=\underset{\alpha _1^{}\alpha _2^{}}{}L_{\alpha _1\alpha _2,\alpha _1^{}\alpha _2^{},}\rho _{\alpha _1^{}\alpha _2^{}}$$
(2)
the usual Weyl-Wigner transform

$$u_{\alpha _1\alpha _2}(𝐫,𝐤)=𝑑𝐫^{}\varphi _{\alpha _1}\left(𝐫+\frac{𝐫^{}}{2}\right)\frac{e^{i𝐤𝐫^{}}}{(2\pi )^{\frac{3}{2}}}\varphi _{\alpha _2}^{}\left(𝐫\frac{𝐫^{}}{2}\right)$$
(3)
one gets:

$$\frac{d}{dt}f(𝐫,𝐤)=L(𝐫,𝐤;𝐫^{},𝐤^{})f(𝐫^{},𝐤^{})𝑑𝐫^{}𝑑𝐤^{},$$
(4)
where

$$f(𝐫,𝐤)=\underset{\alpha _1\alpha _2}{}u_{\alpha _1\alpha _2}(𝐫,𝐤)\rho _{\alpha _1\alpha _2}$$
(5)
is the well-known Wigner function WF ; Buot  while the quantity $`L(𝐫,𝐤;𝐫^{},𝐤^{})`$ is the Liouville operator in the new phase-space representation $`𝐫,𝐤`$.

For a closed system, the Wigner function $`f`$ is defined for any value of the real-space coordinate $`𝐫`$ and its time evolution is fully dictated by its initial condition. In contrast, for a system with open boundaries $`f`$ is defined only within a given region $`\mathrm{\Omega }`$ of interest and its time evolution is determined by the initial condition at $`t_{}`$ inside such region as well as by the value of its incoming part on the boundary $`𝐫_b`$ of the domain at any time $`t^{}>t_{}`$. Such boundary-condition scheme Modena  —also referred to as “U scheme”— is depicted in Fig. 1. More specifically, in order to properly impose the desired spatial boundary conditions, let us rewrite Eq. (4) as

$$\frac{d}{dt}f(𝐫,𝐤)=\stackrel{~}{L}(𝐫,𝐤;𝐫^{},𝐤^{})f(𝐫^{},𝐤^{})𝑑𝐫^{}𝑑𝐤^{}+S(𝐫,𝐤)$$
(6)
by adding and subtracting a source term $`S(𝐫,𝐤)=v(𝐤)f^b(𝐤)\delta (𝐫𝐫_b)`$ note-S . The latter is compensated by a renormalization $`\mathrm{\Delta }L=\stackrel{~}{L}L`$ of the Liouville operator $`L`$ in (4) note-sbc :

$$\mathrm{\Delta }L(𝐫,𝐤;𝐫^{},𝐤^{})=v(𝐤)\delta (𝐫𝐫_b)\delta (𝐫𝐫^{})\delta (𝐤𝐤^{}).$$
(7)
By applying the inverse of the Weyl-Wigner transform in (3) to Eq. (6) we get:

$$\frac{d}{dt}\rho _{\alpha _1\alpha _2}=\underset{\alpha _1^{}\alpha _2^{}}{}\stackrel{~}{L}_{\alpha _1\alpha _2,\alpha _1^{}\alpha _2^{}}\rho _{\alpha _1^{}\alpha _2^{}}+S_{\alpha _1\alpha _2}.$$
(8)
Equation (8) is the desired generalization to open systems of the SBE in (2) note-oc .

In order to validate the theoretical approach presented so far, we shall focus on a very simple semiconductor nanostructure: a single-barrier equidistant from the device contacts (see Fig. 2). As basis states $`\alpha `$ we adopt the scattering states of the device potential profile; moreover, to better identify the role played by carrier injection, we shall neglect all other sources of energy relaxation/dephasing in the device active region, like carrier-phonon and carrier-carrier scattering. Under these assumptions, Eq. (8) in steady-state conditions reduces to:

$$\frac{i}{\mathrm{}}\left(ϵ_{\alpha _1}ϵ_{\alpha _2}\right)\rho _{\alpha _1\alpha _2}\underset{\alpha _1^{}\alpha _2^{}}{}\mathrm{\Delta }L_{\alpha _1\alpha _2,\alpha _1^{}\alpha _2^{}}\rho _{\alpha _1^{}\alpha _2^{}}=S_{\alpha _1\alpha _2},$$
(9)
where $`ϵ_\alpha `$ denotes the energy of the scattering state $`\alpha `$ and $`\mathrm{\Delta }L_{\alpha _1\alpha _2,\alpha _1^{}\alpha _2^{}}`$ is the renormalization term in (7) written in the density-matrix representation $`\alpha _1\alpha _2`$.

Figure 2 shows results for the single-barrier potential profile when carriers are primarily injected from left. Here, the simulated real-space charge distribution obtained from the phenomenological injection model in Eq. (1) (dashed curve) is compared to that of the microscopic model in (9) (solid curves). As we can see, the two models give completely different results. The phenomenological model gives basically what we expect: since we have significant carrier injection from left only and since the potential barrier is relatively high, the carrier distribution is mainly located on the left side. In contrast, the microscopic model gives an almost symmetric charge distribution. In order to understand the origin of this unphysical result, let us focus on the nature of the source term in (8). Contrary to the phenomenological injection/loss term in (1), the latter is intrinsically non-diagonal, i.e., the injection of a carrier with well-defined wavevector $`𝐤`$ is described by a non-diagonal source contribution $`S_{\alpha _1\alpha _2}`$. In other words, we inject into the device active region a coherent superposition of states $`\alpha _1`$ and $`\alpha _2`$, in clear contrast with the idea of injection from a thermal —i.e., diagonal— charge reservoir. More specifically, in this case the generic scattering state $`\alpha `$ on the left comes out to be an almost equally weighted superposition of $`+k`$ and $`k`$: $`\varphi _k(z)=a_ke^{ikz}+b_ke^{ikz}`$. This, in turn, tells us that the generic plane-wave state $`k`$ injected from the left contact is also an almost equally weighted superposition of the left and right scattering states. This is why the charge distribution (solid curve in Fig. 2) is almost symmetric: any electron injected from left couples to left as well as to right scattering states. The anomaly of the microscopic model is even more pronounced if we look at the carrier distribution in momentum space (see inset in Fig. 2). While for the phenomenological model (dashed curve) we get a positive-definite distribution showing, as expected, the two symmetric wavevector components of the scattering state, the microscopic result is not positive definite; this tells us that our boundary-condition scheme does not ensure that the solution of Eq. (6) is a “good” Wigner function.

The scenario previously discussed is highly non-physical; it can be mainly ascribed to the boundary-condition scheme employed so far, which implies injection of plane-wave electrons \[see source term in Eq. (6)\], regardless of the shape of the device potential profile. This is an intrinsic limitation of the conventional Wigner-function representation $`𝐫,𝐤`$. It is then clear that, in order to overcome this limitation, what we need is a boundary-condition scheme realizing diagonal injection over the scattering states $`\alpha `$ of the device potential profile.

To this end, in this Letter we propose a generalization of the Wigner-function formulation considered so far. The key idea is to extend the Weyl-Wigner transform in (3) from the $`𝐤`$ to a generic basis set $`\{|\beta \}`$ according to note-u :

$$u_{\beta _1\beta _2}^{\alpha _1\alpha _2}(𝐫)=\mathrm{\Omega }𝑑𝐫^{}\varphi _{\alpha _1}\left(𝐫+\frac{𝐫^{}}{2}\right)\chi _{\beta _1}^{}\left(𝐫+\frac{𝐫^{}}{2}\right)\chi _{\beta _2}\left(𝐫\frac{𝐫^{}}{2}\right)\varphi _{\alpha _2}^{}\left(𝐫\frac{𝐫^{}}{2}\right).$$
(10)
In analogy to (5), our generalized Wigner function is given by note-WF-new :

$$f_{\beta _1\beta _2}(𝐫)=\underset{\alpha _1\alpha _2}{}u_{\beta _1\beta _2}^{\alpha _1\alpha _2}(𝐫)\rho _{\alpha _1\alpha _2}.$$
(11)
By adopting as basis states $`|\beta `$ the scattering states of the device potential profile and assuming a diagonal source term of the form

$$S_{\beta _1\beta _2}(𝐫)=v_{\beta _1}f_{\beta _1}^b\delta _{\beta _1\beta _2}\delta (𝐫𝐫_b),$$
(12)
its equation of motion will be

$$\frac{d}{dt}f_{\beta _1\beta _2}(𝐫)=\underset{\beta _1^{}\beta _2^{}}{}𝑑𝐫^{}\stackrel{~}{L}_{\beta _1\beta _2,\beta _1^{}\beta _2^{}}(𝐫,𝐫^{})f_{\beta _1^{}\beta _2^{}}(𝐫^{})+S_{\beta _1\beta _2}(𝐫)$$
(13)
with a renormalization $`\mathrm{\Delta }L_{\beta _1\beta _2,\beta _1^{}\beta _2^{}}(𝐫,𝐫^{})`$ given by

$$v_{\beta _1}\delta _{\beta _1,\beta _2}\delta _{\beta _1\beta _2,\beta _1^{}\beta _2^{}}\delta (𝐫𝐫_b)\delta (𝐫𝐫^{}).$$
(14)
We stress that now the source term in (12) describes diagonal injection over the scattering states (with velocity $`v_\beta `$), as requested. Indeed, if we now integrate Eq. (13) over the real-space coordinate $`𝐫`$, we get again the density-matrix equation in (8), but with a diagonal source term $`S_{\alpha _1\alpha _2}=v_{\alpha _1}f_{\alpha _1}^b\delta _{\alpha _1\alpha _2}`$ and a much simpler —i.e., partially diagonal— renormalization term $`\mathrm{\Delta }L_{\alpha _1\alpha _2,\alpha _1^{}\alpha _2^{}}=v_{\alpha _1}\delta _{\alpha _1\alpha _2}u_{\alpha _1\alpha _2}^{\alpha _1^{}\alpha _2^{}}(𝐫_b)`$. In the scattering-free case, the stationary solution is again described by Eq. (9). However, due to the diagonal nature of the new source term as well as of the partially diagonal structure of $`\mathrm{\Delta }L`$, Eq. (9) has now a diagonal solution: $`\rho _{\alpha _1\alpha _2}=f_{\alpha _1}\delta _{\alpha _1\alpha _2}`$. More specifically, the diagonal density-matrix elements $`f_\alpha `$ obey the following steady-state equation:

$$𝒯_{\alpha \alpha ^{}}f_\alpha ^{}=f_\alpha ^b$$
(15)
with $`𝒯_{\alpha \alpha ^{}}=u_{\alpha ^{}\alpha ^{}}^{\alpha \alpha }(𝐫_b)`$. Equation (15) is semiclassical in nature, i.e., it involves diagonal density-matrix terms only note-limit .

Figure 3 shows again results for the single-barrier potential profile previously considered. Here, the simulation based on the phenomenological injection model in Eq (1) (dashed curves) is compared to that of the new microscopic model in (15) (solid curves). As we can see, the highly non physical behaviours of Fig. 2 (solid curves) have been completely removed. Indeed, the momentum distribution in the inset is always positive-definite and the two models exhibit a very similar behaviour. We find relatively small deviations close to the device spatial boundaries, which can be ascribed to the interlevel injection coupling $`𝒯_{\alpha \alpha ^{}}`$ \[see Eq. (15)\], not present in the phenomenological injection model.

In conclusion, we have proposed a rigorous description of quantum-transport phenomena in systems with open boundaries. Our analysis has shown that the conventional Wigner-function formalism leads to unphysical results, like injection of coherent superpositions of states from the device spatial boundaries. This basic limitation has been removed by introducing a generalization of the standard Wigner-function formulation, able to properly describe the incoherent nature of carrier injection. The proposed theoretical scheme constitutes a rigorous derivation of the phenomenological injection models commonly employed in the simulation of open quantum devices.

###### Acknowledgements.

We are grateful to Aldo Di Carlo, Massimo Fischetti, Rita Iotti, Carlo Jacoboni, Paolo Lugli, and Paolo Zanardi for a number of stimulating and fruitful discussions.