text:
# Orbital migration and the frequency of giant planet formation

## 1 Introduction

Extrasolar giant planets (EGPs) have been detected by the radial velocity method at orbital distances from several AUs (e.g., 47 UMa c, Gl614b, $`ϵ`$ Eri b, 55 Cnc d) to several hundredths of an AU (51 Peg b et al.) from their central stars (see Mayor & Queloz mq (1995); Butler & Marcy bm (1996); Butler et al. butler97 (1997); Noyes et al. noyes (1997); Cochran et al. cochran (1997); Santos et al. santos (2000); Vogt et al. vogt (2000); Hatzes et al. hatzes (2000); Butler et al. butler2001 (2001); Santos et al. santos2001 (2001); Fischer et al. fischer2002 (2002); Marcy et al. marcy2002 (2002); and many others; see also the review by Marcy et al. 2000a ). Unless giant planets form in place within 1 AU of low mass stars – unlikely in the context of published formation models (see, for example, Guillot et al. guillot96 (1996)) – the observed range of orbital semi-major axes implies that dramatic orbital changes occur after formation. Previous work showed that in gaseous disks and even subsequent particulate disks, giant planets can move from formation distances of around 5 AU to a wide range of final distances, a process known as orbital migration (Lin & Papaloizou lp86 (1986); Lin et al. linetal (1996); Takeuchi et al. takeuchi (1996); Ward 1997a ; 1997b ; Trilling et al. trillingetal (1998); Murray et al. murray (1998); Bryden et al. bryden (1999); Kley kley99 (1999, 2000); Del Popolo et al. delp (2001); Tanaka et al. tanaka (2002)). The first direct determination of the radius of an extrasolar giant planet via transit measurement (Charbonneau et al. c00 (2000), Henry et al. h00 (2000)) supports rapid inward migration of giant planets, since the large planetary radius requires close proximity to the parent star when the planet’s internal entropy was much larger than at present (Guillot et al. guillot96 (1996), Burrows et al. burrad (2000)), and formation in place is generally considered implausible (see, e.g., Guillot et al. guillot96 (1996)). Early inward migration is commensurate with the idea of migration caused by disk-planet interactions, as disk lifetimes are not longer than $`10^7`$ years (e.g., Zuckerman et al. zuck (1995)).

In this work, we address the following questions: (1) How efficient is orbital migration? (2) What population of planets survives the migration process? and (3) How does this produced population compare to the observed EGP population? We answer these questions by allowing a large population of giant planets to evolve and migrate in circumstellar disks with various initial conditions (one planet per disk) and determining the final semi-major axis and mass distributions. We do not employ any stopping mechanisms, but instead allow only those planets to survive whose migration timescales are longer than their disks’ lifetimes. We compare the final semi-major axis and mass distributions of the surviving model planets to those of the observed EGPs. We reproduce the fraction of planets in small orbits, and predict the frequency of and orbital distributions for giant planets that are as yet unobservable. Of those stars which do form giant planets, a few percent should have giant planets ultimately residing in small orbits, and around one quarter of stars which form planets retain planets in orbits of semi-major axis several AU and beyond. Given the current discovery statistics and their uncertainties, our results imply that giant planet formation is very efficient around low-mass stars: we find that 10% to 80% of young stars form planets, with the uncertainty dominated by the initial planetary mass distribution (and with some poorly-known uncertainties associated with discovery statistics). The presently detectable portion of the giant planet population represents only about 25% to 33% of the total extant population of giant planets.

The model we present here is a simple one with a minimum set of physical assumptions. By neglecting specific stopping mechanisms, we provide results that represent a baseline for planet formation statistics, bereft of specific additional assumptions that any given stopping mechanism would require. It is hoped that the calculations and results presented here are transparent enough that others can use them to explore the particular effects that specific assumed stopping mechanisms would have. In short, we use the simplest possible model in order to explore the consequences of and for planet formation.

## 2 Assumptions in the model

In order to approach these problems in a simple and clear way, we must make several assumptions about the problem of planetary migration. These assumptions are listed here and explained below: (1) we use the simple impulse approximation for the Type II migration; (2) we employ no ad hoc mechanisms to stop planetary migration; (3) we assume that gas accretion onto planets migrating in a gap is small, and neglect this accretion; and (4) we neglect Type I migration.

### 2.1 Impulse approximation

In the migration model presented here, we use the impulse approximation described in Lin & Papaloizou (lp86 (1986)) to calculate the torques between the planet and disk and determine the planet’s orbital motion. We have previously shown (Trilling et al. trillingetal (1998)) that the impulse approximation yields quantitatively similar results to much more computationally intensive methods of calculating orbital migration (for example, the WKB approximation of Takeuchi et al. takeuchi (1996)) in the dissipation regime relevant to this work. Given our current knowledge of the migration process, we prefer to explore the consequences of the simplest possible model rather than to add complexity which does not guarantee more accurate results.

### 2.2 Stopping mechanisms

Halting a giant planet’s inward migration at small semi-major axes is even less well understood than initiating and maintaining migration. Various stopping mechanisms have been proposed, including tides (Lin et al. linetal (1996), Trilling et al. trillingetal (1998)); magnetic cavities (Lin et al. linetal (1996)); scattering by multiple planetesimals (Weidenschilling & Marzari weid (1996)); resonant interactions by two (or more) protoplanets (Kley kley2000 (2000)); and mass transfer from the planet onto the star (Trilling et al. trillingetal (1998)). All stopping mechanisms have problems: Tides are insufficient because the magnitude of the tidal torque is hardly ever enough to counterbalance the inward torque; magnetic cavities have not been demonstrated to work in a quantitative way; the presence of more than one planetesimal or planet requires a great deal of serendipity; and mass transfer only works in a small percentage of cases. Additionally, stopping mechanisms must not only produce 51 Peg b-like planets, at small orbital distances, but must also produce planets like HD108147b at 0.1 AU, HD28185b at 1.0 AU, and 47 UMa c at 3.8 AU (Pepe et al. pepe (2002); Santos et al. santos2001 (2001); Fischer et al. fischer2002 (2002)). None of the above mechanisms adequately halts planets at the wide range of observed distances from the parent stars.

In this work, we presume that none of the above mechanisms are effective in stopping more than a small fraction of migrating planets, and instead consider planets stopping for what is essentially a statistical reason: the planet stops migrating when the disk dissipates. We show that by considering a physically reasonable distribution of initial conditions for planets and disks, our model produces a range of final orbital semi-major axes for extrasolar planets – including small, 51 Peg b-like separations – without resorting to any exotic stopping mechanisms. It is not necessary to produce many planets very close to their stars: the current detection rate suggests that approximately 1% of low-mass main sequence stars (F, G, K dwarfs) have extrasolar giant planets in small orbits (Vogt et al. vogt (2000)), and extrasolar planets in all orbits exist around perhaps a few percent of all stars in the radial velocity surveys (see, e.g., Marcy et al. 2000a ). We note in passing that one M-dwarf exhibits two companion objects of at least 0.5 and 2 Jupiter masses (Marcy et al. marcy98 (1998); Delfosse et al. delfosse (1998); Marcy et al. 2001a ); the statistics for M dwarfs are thus very poor and indeed the environment for planet formation may be very different from that of the more massive dwarfs. In what follows we confine ourselves to F through K stars and refer to these conveniently as low-mass stars.

### 2.3 Planetary mass accretion during migration

Giant planets form gaps in their gaseous circumstellar disks; when this occurs, $`\dot{M_p}`$, the planetary mass accretion rate, drops substantially from its pre-gap value as the local surface density is substantially reduced. It is an outstanding question how much, if any, material may flow across the gap from the disk and be accreted onto the planet (see Bryden et al. bryden (1999); Kley kley99 (1999); Bryden et al. bryden2000 (2000); Kley kley2000 (2000); Lin et al. linppiv (2000); Kley et al. kley2001 (2001)). The results shown in those works suggest that accretion timescales ($`\tau _{acc}M_p/\dot{M_p}`$) are between $`3\times 10^5`$ and $`3\times 10^6`$ years for the viscosities, mass ratios, and other parameters we use in our models. We have adopted a characteristic accretion timescale of $`10^6`$ years and performed several tests to determine the effect of including mass accretion in the migration models (Figure 1). Note that the accretion timescale is comparable to the migration timescale; we allow both linear and exponential growth with this characteristic timescale. The results from these test cases show that when all disk parameters are held constant, the range of variation between the migration timescales without mass accretion and timescales with mass accretion (all cases, linear and exponential) is around 25% from the no-growth migration timescale of $`10^6`$ years (Figure 1).

Including mass accretion for a population of migrating planets has a small effect on the overall results, at best. For more massive planets, which migrate quite slowly anyway (see Sections 3.2 and 4), the effect of considering mass accretion is very much smaller, or completely negligible. Additionally, the primary parameters which drive migration (diskmass, viscosity, disk lifetime) are not known to better than 25%, and even the mass accretion rate itself is a product of some underlying assumptions about the flow of material near, around, and onto a migrating planet. Attention to the detail of mass accretion during migration is not warranted in the face of these larger unknowns; instead, we choose to use a simpler model (ignoring post-gap formation accretion) which may neglect some details for the sake of clarity and understanding. Lastly, because migration timescales are only affected by mass accretion in a minor way, the overall semi-major axis and final mass distributions of the statistical study (see below) do not change significantly in the presence or absence of mass accretion. We therefore neglect mass accretion for the remainder of the present work; all subsequent discussions concerning the statistical study implicitly refer to model runs without mass accretion.

### 2.4 Type I migration

Giant planets, because of their large ratio of hydrogen and helium relative to heavier elements, must form in gaseous protoplanetary disks. Regardless of whether a giant planet forms by core accretion (e.g., Wuchterl et al. wuchterl (2000)) or by direct collapse (e.g., Boss boss97 (1997, 1998, 2000, 2001)), when a planet grows to a certain mass ($``$10-30 $`M_{Earth}`$, although there is some debate about what this critical mass actually is), a gap is formed in the disk and spiral density waves travel away from the planet (see, for example, Goldreich & Tremaine gt (1980); Lin & Papaloizou lp86 (1986); Ward & Hourigan hourigan (1989); Lin & Papaloizou lp93 (1993)). Interactions of these waves at the planet’s Lindblad resonances transfer angular momentum between the planet and the disk; this angular momentum exchange in turn produces a change in orbital distance of the planet. Usually (but not always), these interactions lead to a net loss of angular momentum and a decrease in the planet’s orbital semi-major axis. Such planet-disk interactions have been described in detail (see, for example, Goldreich & Tremaine gt (1980); Ward & Hourigan hourigan (1989); Lin & Papaloizou lp86 (1986, 1993); Lin et al. linppiv (2000); Ward & Hahn wardppiv (2000)). The type of migration executed with a planet in a gap is referred to as Type II (see Ward & Hahn wardppiv (2000)).

So-called Type I migration is executed by bodies which are not massive enough to form gaps in their disks (i.e., less than $``$10-30 $`M_{Earth}`$) (Ward 1997a ; 1997b ). Ward has shown that Type I migration may be very fast (less than $`10^6`$ or even $`10^5`$ years) which is obviously detrimental to planet growth and survival: if migration timescales are so short for Earth mass objects, it is unclear how any planet could survive long enough to accrete enough gas mass to open a gap and initiate the slower Type II migration (which itself is perilous for planet survival). Clearly, there is a problem with planet survival in the face of Type I migration. We do not intend to solve this problem in this paper. Instead, we neglect Type I migration for the following reasons.

If Type I migration were really as dangerously efficient as has been proposed, then formation and survival of planets (particular at distances greater than 0.1 AU) would be nearly impossible, since the Type I migration timescales are shorter than any possible disk lifetime. But, since we know that planets survived both at small and intermediate semi-major axes (many EGPs) and also relatively large semi-major axes (Jupiter), Type I migration cannot be as universally destructive as proposed. For lack of a better quantification of the Type I effect, in the present work we neglect the effect of Type I migration on the initial conditions used in our models. Furthermore, all planets we consider in this paper are too massive to undergo Type I migration. Thus, we can safely neglect Type I effects in our calculations. We reiterate that this does statement does not refute Type I migration nor the importance of understanding and characterizing its timescales; however, it is clear that for planets which have reached $``$1 Jupiter mass via any formation mechanism, Type I migration is no longer relevant.

## 3 Planet migration and the statistical study of planet survival

### 3.1 The migration model and disk properties

We have applied a model of disk-induced migration (Trilling et al. trillingetal (1998)) to a set of planetary systems whose initial characteristics are distributed over a representative range of planet formation environments. We use the same method and model as used in Trilling et al. (trillingetal (1998)), which includes torques on the planet due to disk-planet interactions, due to star-planet tidal torques, and due to planetary mass loss. Mass loss significantly prolongs planetary lifetimes for only a very small population of model planets, but is included here for completeness. Star-planet tidal torques affect a larger population, but we shall see that, by far, the largest effect is the magnitude of the disk-planet torque. A planet’s survival, which is enabled essentially entirely by dissipation of the gaseous disk before infall onto the star, depends on its mass, and on the disk mass, disk lifetime, and disk viscosity.

We have performed a statistical study of planet fates in systems in which each of disk mass, disk lifetime, and disk viscosity is varied in turn. Initial planet masses are taken to be integral values between 1 and 5 Jupiter masses (1 Jupiter mass = 1 $`\mathrm{M}_\mathrm{J}`$ = $`2\times 10^{30}`$ g), inclusive, evenly distributed, with one planet per system. Planetary accretion after gap formation is neglected (see Section 2.3). For the three varied parameters, we adopt nominal values and vary each parameter over a range spanning a factor of 30 larger and smaller than the the nominal value, with 200 model runs evenly spaced per decade in log space, for a total of 600 model runs per planet per varied parameter. For example, in varying disk mass, we adopt a nominal disk mass of 0.02 $`M_{\mathrm{}}`$ (after Beckwith et al. beckppiv (2000)), and calculate model cases between 0.02/30 and $`0.02\times 30M_{}`$, or 0.0007 and 0.6 $`M_{\mathrm{}}`$. This range spans disk masses from less than the minimum mass disk for a system containing a Jupiter to a disk so massive that direct binary star production is more likely than giant planet formation (see, for example, reviews by Beckwith & Sargent bs (1993) and Beckwith et al. beckppiv (2000)). Over this entire range we test 600 model runs in evenly spaced intervals of 0.005 log units. Similarly, we use a nominal disk lifetime and alpha viscosity of $`3\times 10^6`$ years and $`3\times 10^3`$, respectively (after Zuckerman et al. zuck (1995), Shakura & Sunyaev shakura (1973), and Trilling et al. trillingetal (1998)) and test parameter values from disk lifetimes of $`10^5`$ to $`9\times 10^7`$ years and disk viscosities of $`10^4`$ to $`9\times 10^2`$. In all, we carried out 1800 model runs per planet mass and 9000 model runs total; each model run has one parameter varied while the other two are given their nominal values. This very fine sampling of disk parameters allows us to simulate a continuum of physical properties of planet-forming systems. All planets’ initial orbital semi-major axes are 5.2 AU. Giant planets forming farther out would spend more time migrating inward, all else being equal, and hence would have a larger chance of survival (see below).

The result of our computations is a three-dimensional grid of models, based on the nominal values and spanning the range of parameter space, for each initial planet mass. To produce a gaussian statistical population in which the nominal value is most likely and the extreme values are least likely, the significance of each model run was weighted by the gaussian probability of the value of the varied parameter according to the standard gaussian formalism:

$$P(u)=\frac{1}{\sigma \sqrt{2\pi }}e^{\left\{\frac{1}{2}\left(\frac{uu_0}{\sigma }\right)^2\right\}}$$
(1)
where $`u=\mathrm{log}x`$ and $`u_0=\mathrm{log}x_0`$, $`x`$ is equal to the value of the varied parameter, $`x_0`$ is the nominal value for the parameter being varied, and where $`\sigma `$, which describes the width in log space of the gaussian parameter distribution we use, is equal to 0.5. In other words, the half-width of the gaussian physical parameter distribution is a factor of $`\sqrt{10}`$, or around 3.16; thus, a parameter value which is 3.16 times greater or less than the nominal value has a probability one half of the probability that of the nominal value: $`P(u_h)=\frac{1}{2}P(u_0)`$ for $`u_h=u_0\pm 0.5`$ which is the same as $`x_h=x_0(3.16)^{\pm 1}`$.

We use this gaussian weighting in all results discussed in this paper. Weighting is achieved by assigning a significance to each surviving model planet as $`\delta \times P(u)`$ where $`\delta `$ is equal to 1 for planets which survive and 0 for planets which do not. Thus, in production of histograms and discussing the ensemble results for the 9000 model runs, we explicitly incorporate the probability of the planet’s formation conditions and circumstellar disk. Thus, planets from extreme disks which survive but whose initial conditions are extremely unlikely are not overrepresented in the statistical results. Note that because probabilities are assigned after all model runs are completed and are simple functions of the nominal values and a description of a gaussian distribution, other probability weightings can be installed to represent other experiments. Future work includes deriving both $`x_0`$ and $`\sigma `$ more directly from observations (see, e.g., Gullbring et al. gull (1998); Hartmann et al. hartmann (1998)) and testing non-gaussian probability distributions.

Our models calculate orbital evolution for $`10^{10}`$ model years, and planets are defined to have survived when they have non-zero semi-major axes after this time. Very few planets (much less than 1%) undergo significant orbital evolution (due to star-planet tides) after the disk dissipates.

### 3.2 Results of the statistical study

The results from our statistical study are shown in Figures 2 and 3. Of the 9000 initial planets, 30% survive. This survival rate assumes a flat mass function (see below), and is normalized to the total probability of all surviving planets; in this and subsequent discussions, all results are post-weighting. We find that around 0.8% of the initial planets have final semi-major axes less than 0.1 AU. Approximately 6.5% have final semi-major axes between 0.1 and 1 AU. The rest of the surviving planets (23% of the original planets) have final distances from their stars greater than 1 AU. We find that 2.5% of all surviving planets are found at less than 0.1 AU; 21% are found between 0.1 and 1 AU; and 76% of surviving planets have semi-major axes greater than 1 AU. The majority (70%) of the population of initial planets migrates too fast; these planets lose their mass onto the central star (by Roche lobe overflow), and do not survive. As noted in Section 2.2, other stopping mechanisms at the inner edge of the disk may operate and would increase somewhat the fraction of planets retained in orbits with very small semi-major axes. These mechanisms will have a small or negligible effect on the overall statistics, except for the possible existence of a magnetic cavity, which effect on the migration has yet to be quantified.

The most likely fate of a planet which forms near 5 AU is that it migrates inward too fast and ultimately is lost onto the central star. The vast majority of planets that do survive, though, tend to reside at greater than 1 AU from their central stars. Our statistical study is robust because we have sampled planet migration over a broad distribution of disk masses, viscosities, lifetimes and initial planet masses: the variables that are key to migration and its termination. We have also carried out an initial study of the effects of initial semi-major axis on planet migration and survival. While it is not easy or obvious to parameterize these results, it is clear that for identical initial disk conditions, planets which form farther out migrate less rapidly. This is because the surface density is lower at larger semi-major axes, so that the torque between the planet and disk is less. The effect on planet survival is that more planets, of lower masses, survive (see below). Therefore, if 5 AU is the minimum semi-major axis for giant planet formation, then the results presented here represent a minimum for planet survival; conversely, if 5 AU corresponds to the largest semi-major axis for planet formation, then the planet survival found here is the upper limit. We choose 5.2 AU as a nominal middle ground. The fraction of planets which survive, the ratio of planets at small semi-major axes compared to at large semi-major axes, and the overall shape of the surviving planet mass distribution are quite similar even when somewhat different initial semi-major axes are considered (see below).

## 4 Final orbital semi-major axis distribution

We have analyzed our statistical results in terms of the final distribution of orbital semi-major axes of the planets, in comparison to the observed EGP distribution. The results are shown in Figure 4. The observed EGP distribution is relatively flat, that is, equal numbers of planets throughout the range of 0.05 to around 0.7 AU; there is the suggestion of a rise from 0.7 AU through just over 1 AU; and a decrease at greater semi-major axes than just over 1 AU (this decrease is due at least in part to observational bias). (The exception to this apparently flat distribution is the slight enhancement in planets at orbital distances around 0.05 AU. In our previous work (Trilling et al. trillingetal (1998)), we described how tidal and mass loss torques can cause a small “piling-up” of planets near the tidal limit, which is typically around 0.05 AU.) In contrast, the surviving planets in our statistical study are predominantly at large ($`>`$ 1 AU) semi-major axes, including a rise at just under 1 AU to a peak at 2-4 AU. The hint of an upturn in the observed distribution at less than 1 AU echoes the pattern in the model distribution and may suggest the larger population that we predict at greater than 1 AU. The physical explanation for the peak in model planets at separations greater than 1 AU is that the time during which a planet is migrating ($`10^5`$ to $`10^6`$ years) can be relatively short compared to both the lifetime of the disk ($`10^7`$ years) and the timescale required to initiate migration (variable, but at least $`10^6`$ years, or more) (Trilling et al. trillingetal (1998)). Therefore, the number of planets which survive at small distances ($`<`$ 1 AU) is small. The majority of planets that survive do so by creating very large gaps in their disks and by requiring long timescales to excavate the farthest regions of their large gaps, as follows.

More massive planets create larger gaps than smaller planets, everything else being equal (Takeuchi et al. takeuchi (1996); Trilling et al. trillingetal (1998)). For the most massive planets, the length of time it takes to fully form the gap can be quite long, approaching or even surpassing the lifetime of the disk, as such gaps are many AUs wide. The gap is completely formed when it reaches its equilibrium size (see Takeuchi et al. takeuchi (1996) and Trilling et al. trillingetal (1998)). During the stage in which a gap is growing but has not yet reached its full equilibrium size, there is no Type I migration since the planet has already opened a gap; there is also little Type II migration because most Lindblad resonances already fall within the growing gap. Thus, the planet is hardly migrating during this time in which the gap is still growing to its equilibrium size. If the timescale for gap growth is longer than the disk lifetime then the planet will not migrate far. Therefore we find that, overall, close-in surviving model planets have smaller masses and distant surviving models planets have larger masses; this rough trend is seen in the observed data as well (see below). Note that relatively large mass planets (compared to the disk mass) are required, in general, for these arguments about large gap formation to be applicable.

Since the vast majority of planets either migrate too fast and are lost onto the star or remain out near their formation location, disk and planet parameters must be “just right” to halt a planet at small distances ($`<`$ 1 AU). The conclusion of this analysis is that it is very difficult to make planets survive at small distances from the parent star, and that to do so requires a large population of giant planets which still reside at large distances from their parent star as well as a population which has migrated too fast and been lost onto the central star. From our model results, we find that for every planet found at small semi-major axes, around 3 planets must exist at much larger distances, depending on the initial mass distribution (see below). The radial velocity technique is most sensitive to planets at small semi-major axes, because the magnitude of the stellar wobble is greatest for close-in planets. Additionally, since the highest precision radial velocity surveys have relatively short baselines of data, planets with longer periods are only beginning to be identified (e.g., 55 Cnc d, with $``$15 year orbit (Marcy et al. marcy2002 (2002))). Therefore, the majority of the extant population of giant planets we predict would, at present, be beyond the detection capabilities of the radial velocity searches. In other words, a large population of Jupiter-like planets is as yet undetected.

## 5 Implications

### 5.1 Giant planet mass distribution

We have also analyzed our surviving statistical population for its distribution of final planet mass, in comparison to the observed EGP mass distribution. This comparison is complicated by the fact that for most EGPs, the known quantity is not planetary mass but $`M_p\mathrm{sin}i`$. We assume $`\mathrm{sin}i`$ = 1, for simplicity. We shall see in any case that it is not the absolute mass that matters but the slope and shape of the mass distribution curve. The slope and shape will be roughly unchanged between $`M_p\mathrm{sin}i`$ and $`M_p`$ assuming a random distribution of system inclinations relative to the Earth. For a random distribution of inclinations, the average of $`\mathrm{sin}i`$ is $`\pi `$/4. With increasing numbers of extrasolar planets known, the average $`\mathrm{sin}i`$ for the population gets closer to the value for a random distribution. However, because of the detection bias of the radial velocity technique (biased towards detection of systems with $`i`$ near 90 degrees), the average of $`\mathrm{sin}i`$ for the detected systems is actually between $`\pi `$/4 and unity (Marcy et al. 2000a ).

The comparison among mass functions used in this work is shown in Figure 5. By assumption, in our model we have started with equal numbers of planets with masses 1, 2, 3, 4, and 5 $`\mathrm{M}_\mathrm{J}`$ (i.e., after accretion has taken place), or a flat initial mass distribution. As shown in Figures 23, and 5, most of the surviving planets in the statistical study are the more massive planets; planets with the largest initial mass (5 $`\mathrm{M}_\mathrm{J}`$) survive preferentially. The observed EGP mass distribution, however, is biased toward smaller mass planets (Figure 5). Since the radial velocity technique is most sensitive to larger masses, the relative dearth of large mass planets discovered is significant. The present observed EGP population can instead be obtained with an initial mass distribution biased strongly toward smaller mass giant planets. This “preferred” initial distribution is shown in Figure 5 as the dashed line. The preferred mass function is derived from the observed mass distribution and the ratio of initial model planets to surviving model planets:

$$\frac{\mathrm{\#}\mathrm{of}\mathrm{initial}\mathrm{preferred}\mathrm{mass}\mathrm{distribution}\mathrm{planets}\mathrm{with}\mathrm{M}_\mathrm{i}}{\mathrm{\#}\mathrm{of}\mathrm{observed}\mathrm{EGPs}\mathrm{with}\mathrm{M}_\mathrm{i}}=\frac{\mathrm{\#}\mathrm{of}\mathrm{initial}\mathrm{flat}\mathrm{mass}\mathrm{distribution}\mathrm{planets}\mathrm{with}\mathrm{M}_\mathrm{i}}{\mathrm{\#}\mathrm{of}\mathrm{surviving}\mathrm{flat}\mathrm{mass}\mathrm{distribution}\mathrm{planets}\mathrm{with}\mathrm{M}_\mathrm{i}}$$

where initial planet masses are $`\mathrm{M}_\mathrm{i}`$ = 1, 2, 3, 4, and 5 $`\mathrm{M}_\mathrm{J}`$ and mass bins of $`\pm `$0.5 $`\mathrm{M}_\mathrm{J}`$ are used. The preferred initial mass distribution by definition produces a final mass distribution identical to that which is presently observed.

Model planet masses and semi-major axes produced using the preferred the initial mass distribution are shown in Figures 2 and 3. The final distributions of semi-major axis and planet mass produced using the preferred initial mass distribution are shown in Figures 4 and 5. As with a flat initial mass distribution, in order to produce the observed extrasolar planet characteristics from the preferred initial mass distribution, many planets with large semi-major axes must be produced. The preferred initial mass distribution results in 0.2% of all initial planets surviving with separations less than 0.1 AU; 1.6% of all initial planets survive and have 0.1 AU $`<a<`$ 1 AU; and 3.8% of all initial planets survive and have $`a>`$ 1 AU. These percentages correspond to 4.4%, 28%, and 68% of all surviving planets to be found in those three semi-major axis ranges, respectively. Thus, the ratio of as-yet undetected planets to detected planets should be around 2:1, compared to a ratio of 3:1 obtained with the flat mass distribution. Thus, the total percentage of stars with planets must be 3-4 times greater than the percentage given by the current observational results. With the preferred initial mass distribution, 5.6% of all initial planets survive at any semi-major axis; with the flat initial mass distribution, this total survival rate is 31%. With the preferred initial mass distribution, fewer planets survive than in the case of the flat initial mass distribution because the smaller mass planets that dominate the preferred initial distribution preferentially are destroyed<sup>1</sup><sup>1</sup>1 Ward (1997a ; 1997b ) concludes that Type II migration timescales are independent of planet mass. This conclusion is valid in the limit that the planet’s mass is small compared to the disk mass, resulting in a planet which is fully entrained in the viscous evolution of the disk. In the case of a finite mass ratio between planet and disk, however, there is a feedback torque from the planet onto the disk, and the planetary migration timescale is seen to be a function of planetary masses (Lin & Papaloizou lp86 (1986); Trilling et al. trillingetal (1998); Bryden et al. bryden2000 (2000); Nelson et al. nelson (2000); Kley kley2000 (2000); D’Angelo et al. dang (2002)). In the regime in which we are working, planetary migration timescales are functions of planet masses in the sense that, in general, smaller planets migrate faster and more massive planets migrate more slowly. This is partially a function of gap size, as described in Section 4.. Our results imply that there must still be a substantial population of giant planets at large semi-major axis that has not yet been detected. The ratio of extant planets at large semi-major axes to extant planets at small semi-major axes is less when the preferred initial mass distribution is used (around 2 to 1) compared to when the nominal (flat) initial mass distribution is used (closer to 3 to 1).

### 5.2 Planet-forming efficiency

Planet forming efficiency, that is, the number of stars that form giant planets, must be relatively high. This is because many planets which form are subsequently destroyed, and many planets which exist have not yet been detected, as follows:

$$\frac{1\mathrm{planet}\mathrm{detected}}{100\mathrm{stars}}\times \frac{3\mathrm{to}4\mathrm{extant}\mathrm{planets}}{1\mathrm{planet}\mathrm{detected}}\times \frac{20\mathrm{planets}\mathrm{formed}}{1\mathrm{to}6\mathrm{extant}\mathrm{planets}}\frac{1080\mathrm{planets}\mathrm{formed}}{100\mathrm{stars}}.$$

The first term, the number of observed stars which have close-in planets, is from Vogt et al. (vogt (2000)) and Marcy et al. (2000a ); the second term, that 2 to 3 planets exist for each planet which has been discovered so far, is from Sections 4 and 5.1; the third term, the fraction of planets formed which survive, is from Sections 3.2 and 5.1. Note that our model only incorporates 1 planet per star, so this means 10% to 80% of F, G, and K dwarf stars form (at least) one giant planet. This is a very large percentage, which bodes well for observing large numbers of giant planets around young stars (see below). The primary source of uncertainty is the third term, based on the unknown initial mass distribution.

Our model suggests a high efficiency of planet formation. Since planets more easily survive migration in lower mass disks, planet formation in such disks leads to greater planet survival. Additionally, the shape of the preferred initial mass distribution suggests a natural cutoff towards larger mass planets: larger mass planets may be more difficult to form (in agreement with Lineweaver & Grether (linew (2002))). In other words, giant planet formation must be a relatively efficient use of the disk gas: it cannot be a process requiring primarily disks at the upper end of the mass range we have considered, but must also occur in disks with masses closer to the minimum we have selected. Our minimum mass disk is slightly less than the mass of Jupiter itself, and represents an extremum. Exclusion of extreme diskmasses from the cohort changes the overall required efficiency very little because these extreme cases have very low probabilities and therefore do not contribute much to the total survival probability. Overall we require that the small- or moderate-disk masses participate in the formation of giant planets, and these giant planets are then overrepresented in the final distribution because of slower migration rates.

Finally, it is important to emphasize that high planet forming efficiency does not mean that most low mass stars currently have Jupiter-like planets. Planetary migration is very unforgiving: our statistical study shows that around 75% to almost 95% of formed planets migrate too fast and meet their demise at their central stars. A large percentage (10% to 80%) of stars and disks must produce planets so that a few percent may survive close to their central stars and be detected by radial velocity searches. But there should then be a population of giant planets in intermediate (2-4 AU) and larger orbits, around perhaps 2%-3% of F, G, and K dwarf stars. This fraction agrees with the estimates of Zucker & Mazeh (2001b ). Giant planet formation may be occurring around many young stars today, and searches for indications of planet formation around these stars should have a high success rate.

## 6 Discussion and some uncertainties

Our preferred initial mass distribution and mechanism for producing close-in giant planets requires a high frequency of planet formation in disks, and a bias toward making giant planets with smaller masses (in agreement with Lineweaver & Grether linew (2002)). The latter conclusion is robust so long as the distribution of observed system inclinations as seen from Earth is roughly random; that is, as long as the measured $`M_p\mathrm{sin}i`$ yields actual masses on average less than a factor of two larger. To date, five systems with known inclinations (HD209458, 55 Cnc, HD210277, $`\rho `$ CrB, and $`ϵ`$ Eri) all have $`i30`$ degrees and thus have mass factors less than 2 (Charbonneau et al. c00 (2000); Henry et al. h00 (2000); Trilling & Brown 55cnc (1998); Trilling et al. omnibus (2000); Greaves et al. greaves (1998)). Additionally, by tidal arguments, $`\upsilon `$ And, HD75289, HD187123, 51 Peg, and HD217107 are all constrained to have $`i30`$ degrees as well, although there are no direct observations as yet to confirm this (Trilling egptides (2000)); the Hipparcos data implies an inclination $``$25 degrees for $`\upsilon `$ And (Mazeh et al. mazeh (1999)). Additionally, from dynamical and astrometric arguments, the multiple-planet systems Gl876, HD168443, and $`\upsilon `$ And likely have inclinations that are far from face-on (Laughlin & Chambers laugh2001 (2001); Marcy et al. 2001b ; Chiang et al. chiang (2001)). Therefore, most systems for which inclinations are known or suspected appear to have mass factors less than around 2. Astrometric studies of EGP systems show a range of results but few certain inclinations, to date (Zucker & Mazeh zm2000 (2000); Halbwachs et al. halb (2000); Gatewood et al. gatewood (2001); Han et al. han (2001); Pourbaix pour (2001); Pourbaix & Arenou poura (2001); McGrath et al. mcgrath (2002)).

Giant planet formation is as yet not well enough understood to permit prediction of a preferred initial mass distribution as a function of disk mass. Models in which core accretion precedes accumulation of gas require a substantial surface density of solids to permit accretion to occur within the disk lifetime (Lissauer lissauer (1987)), with the resulting gas accretion occurring fairly rapidly once a threshold mass of solids is accumulated (Wuchterl et al. wuchterl (2000)). There is thus no simple and direct link between the disk gas mass and the final masses of the giant planets formed. The presence of Uranus and Neptune in our own system is suggestive of a process that is sensitive to depletion of the gas either due to small available amounts of gas beyond a certain distance from the primary, or to long core accretion times that stretch beyond the disk lifetime. However, alternative models for the formation of Uranus and Neptune exist that do not tie them directly to conditions in the 20-30 AU region of the gas disk (Thommes et al. thommes (1999)).

In sum, we do not as yet have a direct link between the mass distribution of giant planets and the mass of the disk. Thus our preferred initial mass distribution is a general prediction: whatever mechanism forms giant planets tends to prefer the smaller mass objects. This is at least crudely consistent with the inference that the combined giant planet/brown dwarf mass function appears to show a minimum around approximately 10 $`\mathrm{M}_\mathrm{J}`$ (Marcy & Butler marbutl2000 (2000); Zucker & Mazeh 2001a ), suggestive of a planet-forming process biased to lower masses and a brown dwarf-forming process (direct collapse or disk instability) that is most efficient at much higher masses.

In our study we did not consider minimum masses below 1 $`\mathrm{M}_\mathrm{J}`$. The lack of data for objects significantly below $`M_p\mathrm{sin}i`$ of 0.5 $`\mathrm{M}_\mathrm{J}`$ hampers our ability to extend the calculations, and we must await further and more sensitive data (e.g., SIM, see below) to establish the frequency of lower mass planets. However, the radial velocity surveys have begun to detect planets with masses as small as 0.12 $`\mathrm{M}_\mathrm{J}`$ (HD49674b; Butler et al. butler2002 (2002)). This implies that orbital evolution takes place for Saturn-mass and smaller planets. Further observations will show the effects of planet mass on orbital evolution.

Likewise, we have arbitrarily imposed planet formation at 5.2 AU. Formation of giant planets out to double or even quadruple that distance is possible, based on our own solar system. We have performed a small set of model runs for planets with larger initial semi-major axes, out to 10 AU. The overall behavior and survival of these migrating planets is not qualitatively different from that of the entire suite of 9000 models. The general effect of larger initial semi-major axes is an increase in the survival rate (because migration times are longer – see Papaloizou & Larwood paplar (2000)). This would produce a lower required planet-forming efficiency and an initial mass distribution which is biased toward slightly higher masses relative to the preferred initial mass distribution, since more small mass planets survive. Eventually, however, a distance must be reached in all disks at which accretion times become too long and giant planet formation yields first to ice giants (Uranus and Neptune) and then debris. The current cohort of close-in giant planets does not provide a constraint on the outer orbital radius at which giant planet formation could occur. A more complete sample of observed EGPs at larger semi-major axes is required. It is possible that production of giant planets is actually preferred at the position in the disk where water ice first condenses out, because the surface density of ice cold-trapped there (and hence of total solids) is very high (Stevenson & Lunine stevlun (1988)). However, testing this hypothesis must await a more complete sample of planetary systems.

Our model says nothing about the distribution of orbital eccentricities of the known radial velocity companions. Many of the known planets with semi-major axes beyond the tidal circularization radius have significant eccentricities. This suggests a process of orbital evolution commonly associated with the formation and earliest evolution of giant planets, and yet one which failed to act on our own four giant planets. The discussion of whether radial migration in gaseous or particulate disks produces eccentric orbits is an ongoing one (e.g., Ward 1997a ; 1997b ; Bryden et al. bryden (1999); Kley kley2000 (2000); Papaloizou et al. pap2001 (2001); Murray et al. murray2002 (2002)). The presence of high eccentricities suggests that multiple giant planet formation may be the norm, and that interactions among giant planets after disk migration ends pump up the orbital eccentricities. The stochastic nature of such interactions precludes their inclusion in an analysis such as the one we present. Further, until the orbital inclinations can be determined, an important clue to the genesis of the eccentricities is missing. Regardless of the mechanism(s) by which eccentricity is created, except for extreme cases, our results regarding the semi-major axis distribution are hardly affected. Our statistical results regarding planet formation frequency hinge on the frequency of the close-in giant planets versus those farther out. Additional modest semi-major axis evolution associated with eccentricity pumping does not alter our results significantly. Dramatic three-planet interactions that would leave giant planets in vastly altered orbits and greatly affect our statistics are arguably rare. In summary, while we do not calculate eccentricities induced by gas drag in our one-dimensional model, the results based on the statistics of semi-major axis distribution remain robust.

The migration behavior and orbital evolution of planets in a multiple-giant planet system could differ from the results presented here. In these cases, both the gas-planet interactions and the planet-planet interactions could cause orbital evolution. Since we know that systems with multiple giant planets exist (e.g., $`\upsilon `$ And (Butler et al. butler99 (1999)), 47 UMa (Fischer et al. fischer2002 (2002)), HD168443 (Marcy et al. 2001b ), 55 Cnc (Marcy et al. marcy2002 (2002)), and, to a lesser extent, our own), quantitative studies of these systems are needed. For some cases of multi-planet systems, planet survival statistics could be different from our results because of the additional planet-planet interactions. Additionally, if Type II migration should operate significantly more or less efficiently than we have calculated here (either through the presence of another planet or some other disk-planet interaction), the overall survival and statistics we have derived here could change substantially.

## 7 Conclusions

We have presented a simple, baseline model of planet migration. Our simple model allows us to study the statistical behavior of planets migrating in disks without complications introduced from relatively unconstrained processes, like stopping mechanisms. Our overall results can be useful as a baseline statistical result of planet formation and survival. We have shown that most giant planets, under nominal initial conditions, migrate rapidly relative to the disk lifetime and are destroyed before the circumstellar disk dissipates. We have found that the observed extrasolar planets must represent only the tip of the iceberg – perhaps 25% to 33% of the total extant population has been detected – based on the required distribution of initial parameters. Most planets which have formed and survived reside outside the current detection limit of the radial velocity searches (a similar conclusion was reached by Armitage et al. arm (2002)); this prediction is supported with new detections reported in Vogt et al. (vogt2002 (2002)). We have shown that the planetary initial mass distribution must be biased toward smaller masses in order to produce the shape of the observed EGP mass distribution, yielding our preferred initial mass distribution (see also Armitage et al. arm (2002)); again, recent detections by Vogt et al. (vogt2002 (2002)) support the conclusion that the lowest mass planets are the most common. Finally, because migration in high mass disks destroys almost all giant planets formed therein, low mass disks – within a factor of 10 or so of the mass of Jupiter – must be capable of forming giant planets (essentially in agreement with recent observational results by Carpenter (carp (2002))). Hence we find that giant planet formation must be a relatively efficient process in disks. Because high mass planets migrate more slowly than low mass planets in a given disk, overall there should be more massive planets at intermediate and large semi-major axes than found close-in. The population of close-in planets, in contrast, should be dominated by smaller mass planets. Zucker & Mazeh (zm2002 (2002)) have shown that this effect is, in fact, observed and statistically significant.

Migrating giant planets may be detrimental to terrestrial planet survival, if terrestrial planets form coevally with giant planets. Planets interior to a migrating giant planet would be disrupted and lost from the system. This of course assumes that smaller planets do not migrate, although they too likely migrate, potentially on even shorter timescales than giant planets (Ward 1997a ; 1997b ). If terrestrial planets form after the dissipation of gas in the protoplanetary disk, then disruption by a migrating giant planet may be less of a risk (excepting giant planet migration caused by planet-planetesimal interaction (Murray et al. murray (1998))). Kortenkamp & Wetherill (kort (2000)) have considered the case of terrestrial planet formation when Jupiter has both its current and a larger heliocentric distance of 6.2 AU. They have found that accumulation of rocky bodies may be easier with Jupiter at a larger heliocentric distance (that is, pre-migration). It is possible that formation of the terrestrial planets in our Solar System may reveal clues about giant planet migration in our planetary system; certainly, these studies are also relevant to the formation of small, rocky planets in other planetary systems. In the context of large scale migrations, if terrestrial planet formation requires that giant planets have not migrated through the terrestrial zone around 1 AU, then only around 20% (flat initial mass distribution) or 3% (preferred) of all planet-forming systems qualify (i.e., a few percent of all late-type stars; see sections 3 and 5.2). If, however, the constraint is merely that there be no giant planet in the immediate vicinity of the terrestrial planet zone at the onset of terrestrial planet formation (perhaps $`10^7`$ years), then the vast majority (93% for flat initial mass distribution, 99% for preferred) of planet-forming systems qualify. This number corresponds to nearly all planet-forming systems which is around 10% to 80% of late-type stars.

Although radial velocity detections of more distant giant planets will become possible as the time baseline of observations increases, astrometric techniques are more sensitive to giant planets in large orbits. SIM, the Space Interferometry Mission (Danner & Unwin sim (1999)), will do a thorough job of detecting giant planets, Uranus-mass objects, and even smaller bodies from small to large semi-major axes ($``$10 AU), with maximum sensitivity achieved for planets approximately 0.3 $`\mathrm{M}_{\mathrm{Earth}}`$ around 3-5 AU. SIM’s ability to test our predictions of the preferred initial giant planet mass distribution will be limited largely by mission lifetime. However, SIM can also analyze disks around young stars with high precision (target resolution of 1 microarcsecond), perhaps mapping out the signature of gaps created by migrating (or non-migrating) giant planets and giving us a rough snapshot of the time-dependent mass distribution of planets during the migration phase itself. (A 1 AU gap at 100 parsecs is 10 milliarcseconds.) For giant planets in the largest orbits, i.e., 20-30 AU from their parent star, direct imaging techniques may be the only practical method for detection since astrometric techniques would require baselines of decades or more.

Our results suggest that techniques to study planet formation around young stars – radial velocity; high resolution imaging of young stellar systems; searches for gaps such as with SIM and potentially also the Space InfraRed Telescope Facility (SIRTF); searches for planetary outflows (Quillen & Trilling quillen (1998)); or searches for other indirect evidence, like cometesimals scattered onto stars (Quillen & Holman qh (1999)) – should ultimately have a very high success rate. We anticipate that observational data will show that planet formation is taking place around 10% to 80% of low mass pre-main sequence stars, and that planet searches around main sequence stars will have a much lower success rate. Planet formation is an “easy come, easy go” business, with many planets created and many planets destroyed, and with an important minority – including our own Jupiter – surviving.

###### Acknowledgements.

The authors thank Doug Lin, Peter Bodenheimer, Geoff Bryden, Hubert Klahr, Geoff Marcy, Fred Adams, and Norm Murray for many useful discussions. We thank an anonymous referee for useful suggestions. D.E.T. acknowledges a NASA GSRP grant, and NASA grants to Doug Lin, Peter Bodenheimer, and Robert H. Brown. J.I.L. acknowledges support from the NASA Origins Program. W.B. acknowledges support from the Swiss National Science Foundation. D.E.T. thanks E. Stein for hosting him for useful work sessions, and NASA’s IRTF which accommodated him for some of the writing of this paper.


format: text

text:
# Electronic Thermal Conductivity of Multi-Gap Superconductors with Application to MgB2

## Abstract

The remarkable field dependence of the electronic thermal conductivity observed in MgB2 can be explained as a consequence of multi-gap superconductivity. A key point is that for moderately clean samples, the mean free path becomes comparable to coherence length of the smaller gap over its entire Fermi surface. In this case, quasiparticle excitations bound in vortex cores can easily be delocalized causing a rapid rise in the thermal conductivity at low magnetic fields. This feature is in marked contrast to that for anisotropic or nodal gaps, where delocalization occurs only on part of the Fermi surface.

The unexpected discovery of superconductivity in MgB2 with a relatively high $`T_c=38`$ K Nagamatsu01  aroused great interest and was soon followed by experiments which established phonon mediated $`s`$-wave superconductivity, e.g., a B-isotope effect Budko01 , a coherence peak in 11B nuclear relaxation rate Kotegawa01  and an exponential dependence for temperatures $`T10`$ K Yang01 ; Manzano02 . Theoretical studies concluded that the coupling of the holes in the $`2p_\sigma `$-bands of the B-planes to bond stretching modes was strong and primarily responsible for superconductivity. The electron-phonon coupling on the parts of the Fermi surface associated to $`2p_\pi `$-bands is much weaker Kortus01 ; An01 ; Kong01 ; Bohnen01 .

Despite its standard origin, superconductivity in MgB2 has several unusual properties pointing towards a more complex nature. One aspect is the presence of two gaps of different magnitude. Their ratio is estimated as $`r=\mathrm{\Delta }_0^\mathrm{S}/\mathrm{\Delta }_0^\mathrm{L}0.3`$–0.4 based on various experiments Yang01 ; Tsuda01 ; Szabo01 ; Chen01 ; Giubileo01 ; Bouquet01 . Evidence for two gaps is also provided by the rapid rise of the specific-heat coefficient, $`\gamma _s(H)`$, at very low magnetic fields Yang01 ; Bouquet01 . Orbital dependent superconductivity has been proposed theoretically by several authors with the primary (secondary) gap associated with the $`\sigma `$\- ($`\pi `$)-bands Liu01 ; Choi01 .

Recent studies of the inplane thermal conductivity in a magnetic field show an unusual field dependence Sologubenko02 . For fields both parallel and perpendicular to the $`c`$-axis, the electronic thermal conductivity $`\kappa _s(H)`$ exhibits a steep increase in the low-field region, suggesting a large release of mobile quasiparticles in the mixed state. This contrasts strongly with the behavior of conventional $`s`$-wave superconductors, where quasiparticles bound in the vortex cores give very little contribution to $`\kappa _s`$ except very close to $`H_{c2}`$ Lowell70 ; Vinen71 . At first glance, a small secondary gap $`\mathrm{\Delta }_0^\mathrm{S}`$ in multiband models would provide enough carriers for transport at low fields. However, they would be non-mobile carriers inherent in their $`s`$-wave character. It is our aim here to reexamine thermal transport for multi-gap superconductors and show the drastic influence of sample purity on the characteristic behavior of $`\kappa _s(H)`$.

The measured MgB2 samples are regarded as being in the moderately clean regime: experimental estimates of the mean free path give $`\mathrm{}500`$$`800`$ Å, compared with the inplane coherence length $`\xi _{ab}120`$ Å derived from $`H_{c2}`$, which is determined by the primary superconducting $`\sigma `$-band Sologubenko02a . The relevant length scale here to be compared with $`\mathrm{}`$ is, however, $`\xi _0^\mathrm{S}`$ of the passive $`\pi `$-band. Thus, we may consider the quasiparticles in the $`\sigma `$-band in the moderately clean regime, while those in the $`\pi `$-band with $`\xi _0^\mathrm{S}\xi _0^\mathrm{L}/r`$ can be in marginally clean regime. The numerical calculation based on the Bogoliubov-de Gennes framework shows that low-energy states in the smaller gap are loosely bound in vortex cores Nakai02 . Moreover, recent scanning tunneling spectroscopy measurements confirmed a large $`\xi _0^\mathrm{S}500`$ Å in MgB2 Eskildsen02 . Since the magnitude of the secondary gap is small all over the Fermi surface, we expect a distinctively different behavior of $`\kappa _s(H)`$ compared to single-band superconductors with an anisotropic gap or even gap nodes.

In view of these circumstances, we analyze the field dependence of $`\kappa _s`$ and the density of states (DOS), $`N_s`$. For this purpose we introduce Pesch’s solution Brandt67 ; Pesch75 ; Klimesch78  for the quasiclassical formalism Serene83 ; Rammer86 , which is known to be valid for the range of purity in question. We show that the quasiparticle excitations in the small gap bound in the vortex cores can easily be delocalized in the marginally clean regime, causing a rapid rise at low magnetic field. This field dependence is definitely stronger than that obtained for any of the single-band models. On the contrary, superclean samples should exhibit a behavior very similar to that of conventional $`s`$-wave superconductors.

We restrict our considerations to the case of the inplane thermal current with $`Hz`$. Thus inplane impurity scattering is the most important for the thermal transport. The scattering matrix between $`\sigma `$\- and $`\pi `$-bands is assumed to be small Mazin02  because of the different parity of the two orbitals under the reflection $`zz`$. Thus, we neglect interband impurity scattering completely and discuss contributions from each bands independently. In order to calculate $`\kappa _s(H)`$ and $`N_s(H)`$, we introduce the quasiclassical propagators,

$$\widehat{g}(\omega _n,\widehat{𝒌},𝑹)=\left(\begin{array}{cc}g& f\\ f^{}& g\end{array}\right)\frac{i}{\pi }𝑑\xi \widehat{\tau _3}\widehat{G}(\omega _n,𝒌,𝑹),$$
(1)
where $`\widehat{G}`$ is the Nambu-Gorkov Green’s function matrix with the fermionic Matsubara frequency, $`\omega _n`$, the center of mass coordinate, $`𝑹`$, and the relative momentum $`𝒌`$. $`\tau _3`$ is the $`z`$-component of the Pauli matrices acting on the particle-hole space and $`\widehat{𝒌}𝒌_\mathrm{F}/|𝒌_\mathrm{F}|`$ is the unit wave vector at the Fermi surface. They satisfy the normalization condition $`\widehat{g}^2=\widehat{1}`$ and obey the Eilenberger equations ($`\mathrm{}=c=k_\mathrm{B}=1`$ hereafter),

$$[\left(i\stackrel{~}{\omega }_n+e𝒗_\mathrm{F}𝑨(𝑹)\right)\widehat{\tau }_3\widehat{\mathrm{\Delta }}(\widehat{𝒌},𝑹),\widehat{g}]+i𝒗_\mathrm{F}\mathbf{}_𝑹\widehat{g}=0,$$
(2)
supplemented by the gap and Maxwell equations. We introduce the gap matrix

$$\widehat{\mathrm{\Delta }}=\left(\begin{array}{cc}0& \mathrm{\Delta }(\widehat{𝒌},𝑹)\\ \mathrm{\Delta }^{}(\widehat{𝒌},𝑹)& 0\end{array}\right),$$
(3)
and the renormalized frequency, $`\stackrel{~}{\omega }_n=\omega _n+\sigma (\omega _n)`$, where $`\sigma `$ is the diagonal element of the impurity self-energy determined by the Born or the $`T`$-matrix approximation in this study. We neglect vertex corrections.

Instead of solving these transport-like equations self-consistently, we adopt the Brandt-Pesch-Tewordt (BPT) approximation Brandt67 . In this approximation, an Abrikosov solution is used for vortex lattice structures and the spatial dependence of the magnetic field is replaced by the external uniform field $`H`$. Only the uniform component $`\overline{g}`$ is kept, since the higher Fourier $`𝑲`$-components of $`g(𝑹)`$ decrease rapidly as $`\mathrm{exp}(\mathrm{\Lambda }^2K^2)`$, $`\mathrm{\Lambda }=1/\sqrt{2eH}`$ being the magnetic length. On the other hand, the exact spatial dependence of the anomalous propagators is taken into account including the phase variation due to the vortices. Although this theory was designed to work well for $`HH_{c2}`$, especially in strongly type-II superconductors like MgB2, a detailed comparison to numerical solutions yields good agreement both for $`s`$\- and $`d`$-wave superconductors over almost the whole field range Dahm02 . This numerical study also shows that the frequently applied Volovik-theory Volovik93 , yielding $`\gamma _sH`$ for an $`s`$-wave gap and $`\gamma _s\sqrt{H}`$ for gaps with lines of zeros, is restricted to the very low-field region. This indicates the importance of quasiparticle transfer between vortices even in the relatively low-field region Ichioka99 .

By means of BPT, the solutions in eq. (2) can be obtained Pesch75  formally (after analytic continuation) as

$$\overline{g}_{\widehat{𝒌}}(\omega )=\left[1i\sqrt{\pi }[2\mathrm{\Lambda }\overline{\mathrm{\Delta }}(\widehat{𝒌})/v_\mathrm{F}(\widehat{𝒌})]^2W^{}(u)\right]^{1/2},$$
(4)
where $`u=\stackrel{~}{\omega }[2\mathrm{\Lambda }/v_\mathrm{F}(\widehat{𝒌})]`$, $`W(u)=e^{u^2}\mathrm{erfc}(iu)`$ and $`\stackrel{~}{\omega }=\omega +i\sigma `$. Here $`\overline{\mathrm{\Delta }}(\widehat{𝒌})`$ denotes the spatial average of the gap and $`v_\mathrm{F}(\widehat{𝒌})`$ is the component of $`𝒗_\mathrm{F}`$ perpendicular to the field. The real part of $`\overline{g}_{\widehat{𝒌}}(\omega )`$ is nothing but the angle-dependent DOS normalized by the normal-state DOS, $`N_\mathrm{n}`$. In order to get the closed-form solution, we use the Born approximation for the $`s`$-wave scattering self-energy, i.e., $`\sigma (\omega )=\overline{g}_{\widehat{𝒌}}(\omega )/2\tau _n`$, where $`\tau _n`$ is the lifetime in the normal state and $`\mathrm{}`$ represents angular average over the Fermi surface. Then, we can determine the self-consistent $`\sigma (\omega )`$ numerically. From the linear response of the thermal current $`j_{\mathrm{h}i}`$ to the temperature gradient $`_jT`$, we obtain the thermal conductivity tensorKlimesch78 ; Serene83 ; Rammer86  as

$`\kappa _s^{ij}=v_\mathrm{F}^2N_n{\displaystyle _0^{\mathrm{}}}𝑑\omega \left({\displaystyle \frac{\omega }{T}}\right)^2\mathrm{sech}^2\left({\displaystyle \frac{\omega }{2T}}\right)\times `$
$`\text{ }\times \widehat{k}_i\widehat{k}_j\mathrm{Re}[\overline{g}_{\widehat{𝒌}}(\omega )]\mathrm{Re}[\tau _{\widehat{𝒌}}(\omega )].`$ (5)
The comparison with the simple kinetic theory defines the transport lifetime, $`\mathrm{Re}[\tau _{\widehat{𝒌}}(\omega )]`$:

$$\frac{1}{2\tau _{\widehat{𝒌}}(\omega )}=\sigma (\omega )+\sqrt{\pi }\frac{2\mathrm{\Lambda }\overline{\mathrm{\Delta }}^2(\widehat{𝒌})}{v_\mathrm{F}(\widehat{𝒌})}\frac{\mathrm{Re}[W(u)\overline{g}_{\widehat{𝒌}}(\omega )]}{\mathrm{Re}[\overline{g}_{\widehat{𝒌}}(\omega )]}.$$
(6)
Here scattering by the vortices appears in addition to quasiparticle broadening due to impurities. Note that eqs. (4)–(6) can be reduced to the conventional expressions in the $`H=0`$ limit Ambegaokar65 . Moreover one finds $`\overline{g}_{\widehat{𝒌}}=1`$ and $`\tau _{\widehat{𝒌}}=\tau _n`$ in the normal state.

We concentrate on the $`T0`$ limit in this paper. The gap function is factorized as

$$\overline{\mathrm{\Delta }}(\widehat{𝒌})=r\mathrm{\Delta }_0f_{\widehat{𝒌}}\sqrt{1H/H_{c2}},$$
(7)
where $`r`$ represents the smaller gap, ($`0<r<1`$), while $`r=1`$ is used for the larger gap or the single-band case. The shape of the averaged gap function $`\overline{\mathrm{\Delta }}(\widehat{𝒌})`$ is given by $`f_{\widehat{𝒌}}`$, e.g., $`f_{\widehat{𝒌}}=1`$ for an isotropic $`s`$-wave, $`f_{\widehat{𝒌}}=\widehat{k}_x^2\widehat{k}_y^2`$ for $`d_{x^2y^2}`$-wave and $`f_{\widehat{𝒌}}=1/\sqrt{1+a\widehat{k}_z^2}`$ for anisotropic $`s`$-wave Posazhennikova02 . We use the square-root field dependence inferred from the Ginzburg-Landau theory.

We discuss now $`\kappa _s^{xx}(H)`$ for MgB2 and the other cases based on this theory. For MgB2 we use for simplicity a spherical (cylindrical) Fermi surface for the $`\pi `$\- ($`\sigma `$-) band and the parameters $`nN_n^\pi /N_n^\sigma =1.5`$, $`qv_\mathrm{F}^\pi /v_\mathrm{F}^\sigma =1.5`$ and $`r\mathrm{\Delta }_0^\mathrm{S}/\mathrm{\Delta }_0^\mathrm{L}=\mathrm{\Delta }_0^\pi /\mathrm{\Delta }_0^\sigma =0.35`$. The impurity scattering rate for the $`\sigma `$ band is moderate, $`\eta 1/2\tau _\mathrm{n}\mathrm{\Delta }_0^\sigma =0.3`$. These parameters are within the range of current estimates Liu01 ; Sologubenko02 ; Sologubenko02a ; Eskildsen02 . In Fig. 1, the contribution from the $`\pi `$-band shows a rapid rise for very low fields, while that from the $`\sigma `$-band displays rather conventional behavior. This rapid rise is caused by the drastic enhancement of the quasiparticle lifetime of the smaller gap over the entire Fermi surface as vortices are introduced. In contrast, as we demonstrate for anisotropic $`s`$-wave (ani. $`s`$) and $`d_{x^2y^2}`$ ($`d`$) in Fig. 1, the delocalization of quasiparticles occurs only on parts of the Fermi surface. Here, the anisotropy parameter $`a=15`$ was used. We adopted the unitarity limit, $`\delta =\pi /2`$ in the $`T`$-matrix self-energy, i.e., $`\sigma =\overline{g}_{\widehat{𝒌}}/2\tau _\mathrm{n}(\mathrm{cos}^2\delta +\overline{g}_{\widehat{𝒌}}^2\mathrm{sin}^2\delta )`$ for $`d_{x^2y^2}`$-wave Vekhter99 .

The sum of both bands gives $`\kappa _s^{xx}(H)`$ for MgB2 in Fig. 2. The overall features reproduce the experimental data (squares) Sologubenko02  well with the two-band model (2$`s`$) for $`\eta =0.3`$. Similarly, the single-band isotropic $`s`$-wave model with $`\eta =0.08`$ ($`s`$) gives a reasonable fit for Nb (triangles) Lowell70 . The transport properties depend sensitively on the purity of samples, as we can see by considering $`\kappa _s^{xx}(H)`$ of the two-band model for the superclean regime, $`\eta =0.01`$ (dashed line). This shows a behavior similar to that of a conventional $`s`$-wave superconductor (dotted line). In the limit $`\eta 0`$, putting $`\stackrel{~}{\omega }=0`$ in eqs (4)–(6), we obtain the low-field expression for the $`\pi `$-band as $`\kappa _s^{xx}/\kappa _n=(\pi ^{3/2}/5\sqrt{2})(q^2\eta /r^3)H/H_{c2}`$. Thus, even in the case of small $`r`$ the low-field dependence of $`\kappa _s`$ remains small due to the factor $`\eta `$ in the numerator. In other words, in the low-field region the excited quasiparticles are almost localized in the vortex cores even in the case of the smaller gap. On the other hand, the slope of $`N_s(0)`$ is considerably enhanced for small $`r`$ as $`N_s(0)/N_n=(\pi ^2/8\sqrt{2})(q/r)\sqrt{H/H_{c2}}`$. It would be interesting to test this predicted change of behavior for $`\kappa _s(H)`$ in high-quality samples.

The thermal conductivity $`\kappa _s`$ is governed by two characteristic quantities, the DOS and the transport lifetime. We analyze both here in order to elucidate the origin of the above behavior. The field dependence of the DOS is shown in Fig. 3, where all parameters are the same as those used in Fig. 1. The sharp rise of the DOS is consistent with experimental observations of $`\gamma _s(H)`$ in the polycrystalline samples Yang01 ; Bouquet01 . Even though there is a big difference between the two-band model and the single-gap models in $`\kappa _s^{xx}(H)`$, $`N_s(H)`$ shows no drastic differences apart from the presence of a residual DOS in the $`d`$-wave case as $`H0`$.

The appropriate measure of quasiparticle delocalization is the transport lifetime in the plane. We discuss the lifetime of quasiparticles in the passive $`\pi `$ band, where $`r=0.35`$, $`v_\mathrm{F}^\mathrm{S}/v_\mathrm{F}^\mathrm{L}=1.5`$ and $`\widehat{𝒌}z`$. The field dependence of $`\mathrm{Re}[\tau _s/\tau _n]`$ is shown in Fig. 4 for $`\eta =0.01`$, $`0.07`$ and $`0.3`$. As expected, the transport lifetime changes drastically, if the marginally clean regime ($`\eta =0.3`$) is approached, showing a rapid rise in the low-field region. The enhancement of the quasiparticle lifetime occurs over the entire Fermi surface. In addition, the slope of the DOS is much enhanced as shown in Fig 3. These effects yield cooperatively the steep rise in the thermal conductivity shown in Fig. 2. The superclean regime ($`\eta =0.01`$), in contrast, gives only a weak field dependence for low fields due to the quasiparticle localization in this case.

Finally, we comment on the thermal conductivity for $`HH_{c1}`$ which is not covered by our theory. Experimentally, a sudden drop of $`\kappa _s`$ is observed as the magnetic field barely exceeds $`H_{c1}`$Sologubenko02 . The reduction of $`\kappa _s`$ is usually attributed to the decrease of the phonon contribution, since phonons are scattered by the quasiparticles in the vortex cores. Usually this kind of mechanism leads to a more gentle reduction of $`\kappa _s`$ Lowell70 ; Vinen71 . For MgB2 the conditions are more complex. For $`T\mathrm{\Delta }_0^\mathrm{S}`$ all quasiparticle contributions are frozen out in the zero-field limit and they remain localized in the vortex cores for $`HH_{c1}`$. There is a stronger scattering of phonons from core states in multi-gap models. Since the core DOS is considerably larger for the $`\pi `$-band (DOS $`E_\mathrm{F}/(r\mathrm{\Delta }_0^\mathrm{L})^2`$) than for the $`\sigma `$-band Caroli64 . For $`T\mathrm{\Delta }_0^\mathrm{S}`$ the quasiparticles in $`\pi `$-band are sufficiently excited to contribute to the zero-field thermal conductivity. When vortices appear, this quasiparticle contribution is also reduced by scattering at the vortices with localized quasiparticles in the $`\sigma `$-band. This effect in combination with the phonon effect leads to an even stronger drop of $`\kappa _s`$. These simple considerations of the multi-gap effect are in good qualitative agreement with the experiment Sologubenko02 .

In summary, we have discussed the inplane thermal conductivity and the DOS in a magnetic field along the $`z`$-axis in the multi-gap superconductor MgB2. The rapid rise of $`\kappa _s(H)`$ in the low field region is a special feature of a multi-gap superconductor in moderately clean samples. Even in the presence of a small gap, we predict conventional behavior for superclean samples. This sensitivity to sample quality has to be carefully taken into account in the interpretation of thermal transport data in a multi-gap superconductor.

We would like to thank H.R. Ott, I. Vekhter, A.V. Sologubenko, K. Izawa, M. Matsumoto, K. Voelker for fruitful discussions. H.K. is supported by a Grant-in-Aid for encouragement of Young Scientists from Monkasho of Japan. This work was also supported by the NEDO of Japan and Swiss National Fund.


text:
# 1 Introduction

## 1 Introduction

According to Promislow and Pletcher  during the last 20 years there have been incredible advances in our undestanding of the evolution of development, in molecular evolution, and so on. However, theories of ageing have remained relatively restricted to the mutation accumulation theory of Medawar  or to the antagonistic pleiotropy theory of Williams . In their opinion, there are major evolutionary questions that could play a very important role in the evolution of senescence. One of these questions concerns the conflict between sexes.

In flies and many other species females mate several times in a short period, and store the sperm of different males. Thus, males have evolved mechanisms to maximize the ability of their sperm to compete with those of other males, at the expenses of the female’s fitness. For instance, one protein found in the fruit fly male ejaculate has been shown to be similar to a spider neurotoxin .

Here we introduce into the Penna model for biological ageing , which is based on the mutation accumulation theory, such a kind of male strategy, in order to study the consequences on population survival.

## 2 The Penna model

In the asexual version of the Penna model (see  for several applications) the genome of each individual is represented by a computer word of 32 bits. Each bit corresponds to one “year” in the individual lifetime, and consequently each individual can live at most for 32 “years”. A bit set to one means that the individual will suffer from the effects of a deleterious inherited mutation (genetic disease) in that and all following years. One step of the simulation corresponds to reading one bit of all genomes. Whenever a new bit of a given genome is read, we increase by one the individual’s age. The rules for the individual to stay alive are: 1) The number of inherited diseases (bits set to 1) already accumulated until its current age must be lower than a threshold $`T`$, the same for the whole population. 2) There is a competition for space and food given by the logistic Verhulst factor $`V=1N(t)/N_{max}`$, where $`N_{max}`$ is the maximum population size the environment can support and $`N(t)`$ is the current population size. We will consider $`N_{max}`$ five times larger than the initial population $`N(0)`$. At each time step and for each individual a random number between zero and one is generated and compared with $`V`$: if it is greater than $`V`$, the individual dies independently of its age or genome.

If the individual survives until a minimum reproduction age $`R`$, it generates $`b`$ offspring in that and all following years. The offspring genome is a copy of the parent’s one, except for some randomly chosen deleterious mutations introduced at birth. That is, if a bit 1 is randomly tossed in the parent’s genome, it remains 1 in the offspring genome; however, if a bit zero is randomly tossed, it is set to 1 in the mutated offspring genome. In this way, for the asexual reproduction the offspring is always as good as or worse than the parent.

The sexual version of the Penna model was first introduced by Bernardes , followed by Stauffer et al. , who adopted a slightly different strategy. We are going to describe and use the second one. Now individuals are diploids, with their genomes represented by two bit-strings that are read in parallel. In order to count the accumulated number of mutations and compare it with the threshold $`T`$, it is necessary to distinguish between recessive and dominant mutations. A mutation is counted if two bits set to 1 appear at the same position in both bit-strings (inherited from both parents) or if it appears in only one of the bit-strings but at a dominant position (locus). The dominant positions are randomly chosen at the beginning of the simulation and are the same for all individuals.

The population is now divided into males and females. After reaching the minimum reproduction age $`R`$, a female randomly chooses a male with age also equal to or greater than $`R`$ to mate. To construct one offspring genome first the two bit-strings of the mother are cut in a random position (crossing), producing four bit-string pieces. Two complementary pieces are chosen to form the female gamete (recombination). Finally, $`m_f`$ deleterious mutations are randomly introduced. The same process occurs with the male’s genome, producing the male gamete with $`m_m`$ deleterious mutations. These two resulting bit-strings form the offspring genome. The sex of the baby is randomly chosen, each with probability $`50\%`$. This whole strategy is repeated $`b`$ times to produce the $`b`$ offspring. The Verhulst killing factor already mentioned works in the same way as in the asexual reproduction.

## 3 Male selfish strategy and results

To simulate the male ability to increase its own reproductive success we attribute to half of the male population a label, to play the role of a selfish gene. The offspring of a selfish male may inherit this gene with a given probability $`p_s`$. Non-selfish males produce only non-selfish offspring. The selfish male has the advantage that whenever it mates, the female generates $`b+1`$ offspring (instead of $`b`$ as the non-selfish males). On the other hand, the female pays the bill for this male advantage (as usual). We tested two different prices to be paid:

case 1) If a female mates with a selfish male she loses forever part of her energy to dispute for food and space; That is, her probability to die due to the Verhulst factor increases $`20\%`$ (or equivalently, for her the carrying capacity $`N_{max}`$ becomes $`20\%`$ smaller). Such a punishement occurs only once; if she mates again with another selfish male, nothing new happens.

case 2) Whenever a female mates with a selfish male, her genetic death age (programmed since birth as the age where the $`T`$-th genetic disease would appear) is decreased by one, i.e. the female dies one year before the expected date. Note that in this case there is always a price to be paid, that does not depend on the value of any random number, contrary to the previous one.

For both cases we observe that if the selfish gene inheritance probability $`p_s`$ is below a given value $`p_c`$, after many generations no selfish male remains in the population. When $`p_s`$ increases above $`p_c`$, the percentage of selfish males also increases, reaching $`100\%`$ for $`p_s=1`$. In figure 1 we show how this percentage increases with increasing $`p_s`$, for case 1). Squares correspond to a larger population. The curve for case 2) is analogous. There is a dynamical phase transition between an absorbing state in which selfish males are extinct, and a state where they do survive. The order parameter of this transition is the fraction shown in figure 1.

In both cases the total (already stable) population fluctuates around half of the value of the initial population, independently of the $`p_s`$ value. In the first case the male and female population sizes are always the same, as well as their survival probabilities as a function of age. (The survival probability $`S(a)`$ is defined as the ratio between the number of individuals with age $`a+1`$ and the number of individuals with age $`a`$.)

However, in the second case the female survival probability decreases with increasing values of $`p_s`$, in comparison to the male survival probability (which does not change at all). Results are shown in figure 2. Also the female population sizes are slightly smaller than the male ones, when $`p_s>p_c`$. For an initial population of 200,000 individuals (half of each sex), when $`p_s<p_c`$ both the male and female populations stabilize around 87,000. When $`p_s=1`$, the male population stabilizes around 88,000 individuals and the female one, around 84,000.

## 4 Conclusions

We introduce into the Penna model for biological ageing a male strategy that improves its sperm ability to compete with sperm from other males. We obtain that depending on how this strategy works, it may decrease the female survival probability as a function of age. Considering as an order parameter the percentage of males carrying such an ability, we observe a phase transition between populations where no such males survive (after many generations) and populations where they succeed, depending on the probability of their offspring to inherit this “selfish gene”.

Acknowledgements: To PMMH at ESPCI for the warm hospitality, to Sorin Tănase-Nicola for helping us with the computer facilities and to the Brazilian agency FAPERJ for financial support.