text:
# Searching for WIMP Dark Matter: The case for Germanium Ionization Detectors11footnote 1Invited talk given at the 29th International Meeting on Fundamental Physics, Sitges (Barcelona), Feb.2001, to be published in the Proceedings.

## 1 Introduction

Experimental observations and robust theoretical arguments have established that our universe is essentially non-visible, the luminous matter scarcely accounting for one per cent of the critical density of a flat universe ($`\mathrm{\Omega }=1`$). The current picture describes an universe consisting of unknown species of Dark Energy ($`\mathrm{\Omega }_\mathrm{\Lambda }70\%`$) and Dark Matter ($`\mathrm{\Omega }_M2530\%`$) of which less than $`5\%`$ is of baryonic origin. Most of the Dark Matter would then be made of non-baryonic particles filling the galactic halos, at least partially, according to a variety of models. Weak Interacting Massive (and neutral) Particles (WIMPs) are favourite candidates to such non-baryonic components. The lightest stable particles of supersymmetric theories, like the neutralino, describe a particular class of WIMPs. Without entering into considerations about how about how large the baryonic dark component of the galactic halo could be, we take for granted that there is enough room for WIMPs in our halo, to try to detect them, either directly or through their by-products. Discovering this form of Dark Matter is one of the big challenges in Cosmology, Astrophysics and Particle Physics.

WIMPs can be looked for either directly or indirectly. The indirect detection of WIMPs proceeds currently through two main experimental lines: either by looking in cosmic rays experiments for positrons, antiprotons, or other antinuclei produced by the WIMPs annihilation in the halo (CAPRICE, BESS, AMS, GLAST, VERITAS, MAGIC…), or by searching in large underground detectors (SUPERKAMIOKANDE, SNO, SOUDAN, MACRO) or underwater neutrino telescopes (BAIKAL, AMANDA, ANTARES, NESTOR) for upward-going muons produced by the energetic neutrinos emerging as final products of the WIMPs annihilation in celestial bodies (Sun, Earth…). This talk will deal with the direct detection of WIMPs.

The direct detection of WIMPs relies in the measurement of the WIMP elastic scattering off the target nuclei of a suitable detector. Pervading the galactic halos, slow moving ($`300`$ km/s), and heavy ($`1010^3`$ GeV) WIMPs could make a nucleus recoil with a small energy of a few keV ($`\mathrm{T}(1100)`$ keV), at a rate which depends of the type of WIMP and interaction. Only a fraction of the recoil, QT, is visible in the detector, depending on the type of detector and target and on the mechanism of energy deposition. The so-called Quenching Factor, Q, is essentially unit in thermal detectors whereas for the nuclei used in conventional detectors it ranges from about 0.1 to 0.6. For instance for a Ge nucleus only about 1/4 of the recoil energy goes to ionization. On the other hand, the smallness of the neutralino-matter interaction cross-section makes the rates of the nuclear recoil looked for very small. The variety of models and parameters used to describe the Astrophysics, Particle Physics and Nuclear Physics aspects of the process make the neutralino-nucleus interaction rate encompass several orders of magnitude, going from 10 to $`10^5`$ c/kg day, according to the SUSY model and parameters. In fact, it is not higher than $`10^2`$ c/kg day, for the neutralino parameters which provide the most favourable relic density ($`\mathrm{\Omega }\mathrm{h}^2`$). Due to such small rate and small energy deposition produced by the WIMP elastic scattering, the direct search for particle dark matter through their scattering by nuclear targets requires ultralow background detectors of a very low energy threshold. Moreover, the (almost) exponentially decreasing shape of the predicted nuclear recoil spectrum mimics that of the low energy background registered by the detector. All these features together make the WIMP detection a formidable experimental challenge.

Customarily, one compares the predicted event rate with the observed spectrum. If the former turns out to be larger than the measured one, the particle which would produce such event rate can be ruled out as a Dark Matter candidate. That is expressed as a contour line $`\sigma `$(m) in the plane of the WIMP-nucleon elastic scattering cross section versus the WIMP mass. That excludes, for each mass m, those particles with a cross-section above the contour line $`\sigma `$(m). The level of background sets, consequently, the sensitivity of the experiment in eliminating candidates or in constraining their masses and cross sections.

However, this mere comparison of the expected signal with the experimentally observed spectrum will not be able to detect the WIMP, because such spectrum, even extremely small, could still be due to pure background sources. In other words, it would be difficult to be convinced that such signal is due to WIMPs and not to mere background sources. A convincing proof of the detection of WIMPs would need to find unique signatures in the data characteristic of them and not of the background. There exist some temporal or spatial asymmetries specific of the WIMP interaction, which cannot be faked by the background or by instrumental artifacts. They are due to the kinematics of the motion of the Earth (and of our detectors) in the galactic halo. The only distinctive signature investigated up to now is a predicted annual modulation of the WIMP signal rate due to the seasonal variation produced by the Earth’s motion with respect to the Sun. Such a seasonal affect has been found by the DAMA experiment at the $`3\sigma `$ level and has been associated to the existence of a WIMP.

The detectors used so far are: ionization detectors of Ge (IGEX, COSME, H/M, HDMS) and of Si (UCSB), scintillation crystals of NaI (ZARAGOZA, DAMA, UKDMC, SACLAY, ELEGANTS), liquid or liquid-gas Xenon detectors (DAMA, UCLA, UKDMC), calcium fluoride scintillators (MILAN, OSAKA, ROMA), thermal detectors (bolometers) with saphire absorbers (CRESST, ROSEBUD), with telurite absorbers, (MIBETA, CUORICINO) or with germanium absorbers (ROSEBUD) as well as bolometers which also measure the ionization, like that of Si (CDMS) and of Ge (CDMS, EDELWEISS). New detectors and techniques are entering the stage. Worth to be mentioned are: scintillating bolometers of calcium tungstate which measure heat and high (CRESST and ROSEBUD), and of BGO (ROSEBUD); a TPC sensitive to the direction of the nuclear recoil (DRIFT); devices which use superheated droplets (SIMPLE and PICASSO), or those which use colloids of superconducting superheated grains (ORPHEUS). There exist also projects featuring a large amount of target nuclei in segmented detectors, both with ionization Ge detectors (GENIUS, GEDEON) and cryogenic thermal devices (CUORE). Table 4 gives an overview of the experiments on direct detection of WIMPs currently in operation or in preparation. A review of neutralino dark matter can be found in Ref . WIMP direct detection are reviewed, for instance in Ref  and . Recent results are given in Ref .

## 2 Strategies for WIMP detection

The smallness and rarity of WIMP signals dictate the experimental strategies for their detection: reduce first the background, by controlling the radiopurity of the detector, its components, the shielding and the environment. The best radiopurity has been obtained in conventional Ge experiments (IGEX, H/M). In the case of the NaI scintillators (DAMA, UKDM, ANAIS), the achieved backgrounds are still one order of magnitude worse than those reached in Ge, but the use of new radiopure powder in the crystals has lead to substantial improvements (DAMA). The next step is to use discrimination mechanisms able to distinguish electron recoils (tracers of the background) from nuclear recoils (originated by WIMPs or neutrons). Two types of techniques have been applied for such purpose: a statistical pulse shape analysis (PSD) based on the different timing behaviour of both types of pulses and a background rejection method based on the identification (on an event-by-event basis) of the nuclear recoils, by measuring at the same time two different mechanisms of energy deposition, like the ionization (or scintillation) and heat, capitalizing the fact that for a given deposited energy (measured as phonons) the recoiling nucleus ionizes less than the electrons. Examples of PSD are UKDMC, Saclay, DAMA and ANAIS, whereas the event by event discrimination has been successfully applied in CDMS and EDELWEISS by measuring ionization and heat and now in CRESST and ROSEBUD by measuring light and heat.

A promising discriminating technique is that used in the two-phase liquid-gas Xenon detector with ionization plus scintillation, of the ZEPLIN series of detectors. An electric field prevents recombination, the charge being drifted to create a second pulse in addition to the primary pulse. The amplitudes of both pulses are different for nuclear recoils and gammas allowing their discrimination.

Another technique is to discriminate gamma background from neutrons (and so WIMPs) using threshold detectors—like neutron dosimeters—which are blind to most of the low Linear Energy Transfer (LET) radiation (e, $`\mu `$, $`\gamma `$). Detectors which use superheated droplets which vaporize into bubbles by the WIMP (or other high LET particles) energy deposition are SIMPLE and PICASSO. An ultimate discrimination will be the identification of the different kind of particles by the tracking they left in, say, a TPC, plus the identification of the WIMP through the directional sensitivity of the device (DRIFT).

Obviously one should try to make detectors of very low energy threshold and high efficiency to see most of the signal spectrum, not just the tail. That is the case for the bolometer experiments which seen efficiently the energy delivered by the WIMP (quenching factor is unity), (MIBETA, CRESST, ROSEBUD, CUORICINO, CDMS and EDELWEISS).

After the process of reducing and identifying the background is driven to its best, then one should search for distinctive signatures of the dark matter particles to prove that you are, indeed, seeing a WIMP. Identifying labels are derived from the Earth motion through the galactic halo, which produces two assymetries distinctive of WIMPs: First, the Earth orbital motion around the sun has a summer-winter variation which results in a small (5%) annual modulation of the WIMP interaction rates; and second, the Earth (and solar system) motion through the galactic centre produces a large ($`1`$) directional assymetry forward-backward of the recoiling nucleus. The annual modulation signature has been already explored. Pioneering searches for WIMP annual modulation signals were carried out in Canfranc (NaI-32), Kamioka (ELEGANTS) and Gran Sasso (DAMA-Xe). Recently the DAMA experiment at Gran Sasso, using a set of NaI scintillators reported in 1997 and 2000 an annual modulation effect of at the $`3\sigma `$ level interpreted as due to a WIMP of about 60 GeV of mass and scalar cross-section on protons of $`\sigma _\mathrm{p}=7\times 10^6`$ picobarns. As far as the forward-backward assymetry of the signal is concerned the DRIFT project will search for it. Information and current results of the experiments mentioned above can be found in References  and .

In particular, in the underground facility of the Canfranc Tunnel (Spain), various of the techniques mentioned above are currently employed in a search for WIMPs: Ge ionization detectors, (COSME and IGEX), Sodium iodine scintillators (NaI-32 and ANAIS), thermal detectors with Al2O3 and Ge absorbers (ROSEBUD-I) as well as and with calcium tungstate and BGO scintillating bolometers (ROSEBUD-II). To illustrate how these searches proceed we will present the case for Ge detectors, and will describe its status, achievements and results. As one of the main achievements in Ge experiments is their low radioactive background, essential ingredient in the search for dark matter, we will illustrate it with a saga of Germanium experiments performed in Canfranc and the succesive suppresion of the background obtained. Such reduction has provided the lowest raw background rate obtained so far and, consequently, the stringest exclussion plots ever derived with ionization detectors.

## 3 Germanium Experiments

The high radiopurity and low background achieved in Germanium detectors, their fair low energy threshold, their reasonable Quenching Factor (about 25%) (nuclear recoil ionization efficiency relative to that of electrons of the same kinetic energy, or ionization yield) and other nuclear merits make Germanium a good option to search for WIMPs with detectors and techniques fully mastered. The first detectors applied to WIMP direct searches (as early as in 1987) were, in fact, Ge diodes, as by-products of $`2\beta `$-decay dedicated experiments. The exclusion plots $`\sigma `$(m) for spin-independent couplings obtained by former Ge experiments \[PNNL/USC/Zaragoza (TWIN  and COSME-1), UCSB , CALT/NEU/PSI , H/M \] are still remarkable and have not been surpassed till recently by NaI experiments (DAMA) using statistical pulse shape discrimination (PSD) or by cryogenic experiments measuring both heat and ionization like CDMS and EDELWEISS. Table 5 shows the Germanium ionization detector experiments currently in operation.

A Ge detector of natural isotopic abundance (COSME) of the U. Zaragoza / U. S. Carolina / PNNL Collaboration and another one (RG-II) made of enriched 76Ge of the IGEX Collaboration are being used in Canfranc to search for WIMPs interacting coherently with the Ge nuclei of the detectors.

The COSME detector was fabricated at Princeton Gamma-Tech, Inc. in Princeton, New Jersey, using germanium of natural isotopic abundance. The refinement of newly-mined germanium ore to finished metal for this detector was expedited to minimize production of cosmogenic 68Ge. The detector is a p-type coaxial hyperpure natural germanium crystal with a mass of 254 g and an active mass of 234 g which has a long term resolution of 0.43 keV full width at half maximum (FWHM) at the 10.37 keV gallium X-ray. In its first Canfranc installation (at 675 m.w.e) the detector was placed within a shielding of 10 cm of 2000 yr. old (Roman) lead (inner layer) plus 20 cm of low activity lead (about 70 yr old). A 3 mm thick PVC box sealed with silicone closed the lead shielding to purge the radon gas. The PVC box was covered by 1 mm of cadmium and 20 cm of paraffin and borated polyethylene. All the shielding and mounting was supported by 10 cm of vibrational and acoustic insulator sandwiched within two layers of 10 cm of wood mounted on a floor of concrete (20 cm). It was operated in the former Canfranc underground facility at 675 m.w.e. in a small gallery of the Canfranc railway tunnel in the Spanish Pyrenees closed to the traffic. In that set-up, the energy threshold was $`\mathrm{E}_{\mathrm{thr}}=1.6`$ keV and the background at threshold was about 10 counts keV-1 kg-1 day-1. After about 500 days of data taking (referred as COSME-1 experiment) no signal originated by WIMP appeared and, from the background spectrum, an exclusion plot in the cross-section versus WIMP masses plane was derived, assuming WIMP-matter spin-independent couplings. In the COSME-1 data, the energy range chosen to derive the exclusion plot was from 1.6 to 8 keV where the background was, approximately, 5 counts/(keV kg day). In spite of such modest figure, the results after 130 kg day of exposure improved the exclusion plots at low masses (from 9 to 20 GeV) obtained with other Ge experiments because of the low threshold energy of COSME-1.

The COSME detector has been reinstalled, in better background conditions in the new Canfranc Underground Laboratory LSC (at 2450 m.w.e.) inside a Marinelli beaker in Roman lead (COSME-2) in a common shielding (described later) together with the three 2.1 kg enriched germanium detectors of IGEX (the International Germanium Experiment on Double Beta Decay). In its new installation, COSME-2 has an energy threshold of $`E_{thr}=2.5`$ keV, and an energy resolution of $`\mathrm{\Gamma }(\mathrm{FWHM})=0.4`$ keV at 10 keV. The average background rate, in 311 days of exposure (Mt=72.8 kg day) is 0.6 c/(keV kg day) from 2 to 15 keV and 0.3 c/(keV kg day) from 15 to 30 keV which is significantly better (more than one order of magnitude) than in COSME-1. The COSME-2 spectrum is shown in Fig. 1 together with that of COSME-1, to show the remarkable background reduction obtained in the new set-up and shielding and the disappearance of the cosmogenically induced peaks in the 8-10 keV energy region due to the about ten years elapsed between the two experiments. The numerical data of COSME 2 are given in Table 1.

The International Germanium Experiment (IGEX) is an experiment which was designed to investigate the double beta decay of Germanium. It uses three enriched detectors of 76Ge of $``$2.1 Kg each, and it has been described in detail elsewhere . One of these detectors, RG-II, has been recently prepared for use at low energies in typical WIMP searches (referred to as IGEX-DM) .

The IGEX detectors (and so RG-II) were fabricated at Oxford Instruments, Inc., in Oak Ridge, Tennessee. Russian GeO2 powder, isotopically enriched to 86% 76Ge, was purified, reduced to metal, and zone refined to $`10^{13}`$ p-type donor impurities per cubic centimeter by Eagle Picher, Inc., in Quapaw, Oklahoma. The metal was then transported to Oxford Instruments by surface in order to minimize activation by cosmic ray neutrons, where it was further zone refined, grown into crystals, and fabricated into detectors. All of the cryostat parts were electroformed using a high purity OFHC copper/CuSO4/H2SO4 plating system. The solution was continuously filtered to eliminate copper oxide, which causes porosity in the copper. A Ba(OH)2 solution was added to precipitate BaSO4, which is also collected on the filter. Radium in the bath exchanges with the barium on the filter, thus minimizing radium contamination in the cryostat parts. The CuSO4 crystals were purified of thorium by multiple recrystallization.

The RG-II detector, has a mass of $`2.2`$ kg. and an active mass of 2.0 kg, as measured with a collimated source of 152Eu in the Canfranc Laboratory, in agreement with the Oxford Instruments efficiency measurements. The full-width at half-maximum (FWHM) energy resolution of RG-II is 2.37 keV at the 1333-keV line of 60Co. Energy calibration and resolution measurements were made every 7–10 days using the lines of 22Na and 60Co. Calibration for the low energy region was extrapolated using the X-ray lines of Pb. The first-stage field-effect transistor (FET) is mounted on a Teflon block a few centimeters from the central contact of the germanium crystal. The protective cover of the FET and the glass shell of the feedback resistor have been removed to reduce radioactive background. This first-stage assembly is mounted behind a 2.5-cm-thick cylinder of archaeological lead to further reduce background. Further stages of preamplification are located at the back of the cryostat cross arm, approximately 70 cm from the crystal. The detector has preamplifier modified for the pulse-shape analysis.

The detectors shielding (shared by COSME, RG-II and two other IGEX detectors) is as follows, from inside to outside. The innermost shield consists of 2.5 tons of 2000-year-old archaeological lead forming a 60-cm cube and having $`<9`$ mBq/kg of 210Pb(210Bi), $`<0.2`$ mBq/kg of 238U, and $`<0.3`$ mBq/kg of 232Th. The detectors fit into precision-machined holes in this central core, which minimizes the empty space around the detectors available to radon. Nitrogen gas, at a rate of 140 l/hour, evaporating from liquid nitrogen, is forced into the detector chambers to create a positive pressure and further minimize radon intrusion. The archaeological lead block is centered in a 1-m cube of 70-year-old low-activity lead($`10`$ tons) having $`30`$ Bq/kg of 210Pb. A minimum of 15 cm of archaeological lead separates the detectors from the outer lead shield. A 2-mm-thick cadmium sheet surrounds the main lead shield, and two layers of plastic seal this central assembly against radon intrusion. A cosmic muon veto covers the top and sides of the central core, except where the detector Dewars are located. The veto consists of BICRON BC-408 plastic scintillators 5.08 cm $`\times `$ 50.8 cm $`\times `$ 101.6 cm with surfaces finished by diamond mill to optimize internal reflection. BC-800 (UVT) light guides on the ends taper to 5.08 cm in diameter over a length of 50.8 cm and are coupled to Hamamatsu R329 photomultiplier tubes. The anticoincidence veto signal is obtained from the logical OR of all photomultiplier tube discriminator outputs. An external polyethylene neutron moderator 20 cm thick (1.5 tons) completes the shield. The entire shield is supported by an iron structure resting on noise-isolation blocks. The experiment (referred to as IGEX-2000) is located in a room isolated from the rest of the laboratory and has an overburden of 2450 m.w.e., which reduces the muon flux down to a measured $`2\times 10^7\mathrm{cm}^2\mathrm{s}^1`$.

The data acquisition system is based on standard NIM electronics. It has been implemented by splitting the normal preamplifier output pulses of each detector and routing them through two Canberra 2020 amplifiers having different shaping times enabling noise rejection as first applied in Ref . These amplifier outputs are converted using 200 MHz Wilkinson-type Canberra analog-to-digital converters, controlled by a PC through parallel interfaces. For each event, the arrival time (with an accuracy of 100 $`\mu `$s), the elapsed time since the last veto event (with an accuracy of 20 $`\mu `$s), and the energy from each ADC are recorded. For a total counting rate smaller than 1 Hz, the dead time is negligible. The muon veto anticoincidence was done off-line with a software window of 240 $`\mu `$s.

For the data taken in a new set-up, described below, (and referred to as IGEX-2001), the adquision system previously used has been complemented by recording also the pulse shapes of each event before and after amplification by mean of two 800 MHz LeCroy 9362 digital scopes. These are analysed one by one by means of a method based on wavelet techniques which allows us to access the probability of this pulse to have been produced by a random fluctuation of the baseline. This probability is used as a criterium to reject events coming from electronic noise or microphonics. According to the calibration of the method, it works very efficiently for noise events above 4 keV.

## 4 Results obtained from the IGEX set-ups

The detector IGEX-RG-II in both set-ups features an energy threshold of 4 keV and an energy resolution of 0.8 keV at the 75 keV Pb x-ray line. We will show first the results obtained with the set-up and shielding described above. They correspond to 30 days of analyzed data (Mt=60 kg-days) and they will be labelled IGEX-2000 (Table 2 and Fig. 2). The background rate recorded was $`0.4`$ c/(keV-kg-day) between 4–10 keV, $`0.12`$ c/(keV-kg-day) between 10–20 keV, and $`0.05`$ c/(keV-kg-day) between 20–40 keV. Fig. 2 shows the RG-II 30-day spectrum in this shielding. The numerical data are given in Table 2. The corresponding $`\sigma `$(m) exclusion plot is given in Fig. 3 following the procedure explained below. Earlier results from the COSME detector (COSME-1) are included in Fig. 3 (dot-dashed line). The recent results obtained in its current set-up (COSME-2), corresponding to the spectrum of 311 days depicted in Fig. 1 and listed in Table 1 are also shown in its corresponding $`\sigma `$(m) exclusion contour of Fig. 3 (thick-dashed line).

The exclusion plots resulting from the IGEX-2000 data are derived from the recorded spectrum Fig. 2 in one-keV bins from 4 keV to 50 keV. The method followed in deriving the plot has been the same for all the detectors. As recommended by the Particle Data Group, the predicted signal in an energy bin is required to be less than or equal to the (90% C.L.) upper limit of the (Poisson) recorded counts. The derivation of the interaction rate signal supposes that the WIMPs form an isotropic, isothermal, non-rotating halo of density $`\rho =0.3`$ GeV/cm3, have a Maxwellian velocity distribution with $`\mathrm{v}_{\mathrm{rms}}=270`$ km/s (with an upper cut corresponding to an escape velocity of 650 km/s), and have a relative Earth-halo velocity of $`\mathrm{v}_\mathrm{r}=230`$ km/s. The cross sections are normalized to the nucleon, assuming a dominant scalar interaction. The Helm parameterization is used for the scalar nucleon form factor, and the recoil energy dependent ionization yield used is the same that in Ref  $`\mathrm{E}_{\mathrm{vis}}=0.14(\mathrm{E}_{\mathrm{REC}})^{1.19}`$. The exclusion plots obtained from the COSME-2 spectrum (Fig. 1 and Table 1) and from COSME-1 (Fig. 1) have been calculated, as previously noted, in the same way as for IGEX. As shown in Fig. 3, the IGEX-2000 results (thick-solid line) exclude WIMP-nucleon cross-sections above 1.3x10-8 nb for masses corresponding to the 50 GeV DAMA region. Also shown is the combined germanium contour (thin-dashed line), including the last Heidelberg-Moscow data (recalculated from the original energy spectra with the same set of hypotheses and parameters), the DAMA experiment contour plot derived from Pulse Shape Discriminated spectra, and the DAMA region corresponding to their reported annual modulation effect. The IGEX-2000 exclusion contour improves that of other Germanium experiments for masses between 20 GeV and 200 GeV, which includes the mass region corresponding to the neutralino tentatively assigned to the DAMA modulation effect and results from using only raw data without background subtraction. The COSME-2 exclusion contour also slightly improves the Ge-combined plot for masses between 20 and 40 GeV.

Then, a new set-up was constructed, with the purpose of further improving the background. In particular all the detectors except IGEX RG-II were removed from the shielding, leaving RG-II alone. A new, thorough cleaning of the inner layers of the bricks of lead was performed, and the free space left was stretched. In this new set-up, the detector RG-II is surrounded by not less than 40-45 cm of lead of which 25 cm are archaeological. Also the muon veto covers now more completely the ensemble, because three dewars have been removed. A substantied improvement is that of the neutron shielding which has been enlarged to up to 40 cm of thickners , consisting of polyethylene blocks and borated water tanks.

The results obtained whit this new shielding are from a recent run and correspond to an exposure of Mt=80 kg days. The spectrum obtained in this new set-up, and labelled IGEX-2001, is shown in Figure 4 compared with the previous IGEX-2000 spectrum. The numerical data are given in Table 3.

The energy threshold of the detector is, as in previous runnings, 4 keV and the FWHM energy resolution at the 75 keV Pb X-ray line is 800 eV. The background rate recorded was $`0.21`$ c/keV/kg/day between 4–10 keV, $`0.10`$ c/keV/kg/day between 10–20 keV, and $`0.04`$ c/keV/kg/day between 25–40 keV. As it can be seen, the background below 10 keV has been substantially reduced (about a factor 50%) with respect to that obtained in the previous set-up of IGEX-2000, essentially due to the improved shielding (both in lead and in polyethylene-water).This suggests that the neutrons could be an important component of the low energy background in IGEX.

The new exclusion plot derived in this improved conditions IGEX-2001 is shown in Fig. 5 (thick solid line). It improves the IGEX-2000 exclussion contourn (thick dashed line) as well as that of the other previous germanium ionization experiments (and in particular that of the last result of Heidelberg-Moscow experiment now specifically depicted by the thick dotted line) for a mass range from 20 GeV zone to 200 GeV, which encompass, as already said, that of the DAMA mass region. In particular, this new IGEX result excludes WIMP-nucleon cross-sections above 7 $`\times 10^6`$ pb for masses of 40-60 GeV and enters the DAMA region. IGEX excludes the upper left part of this region. That is the first time that a direct search experiment without background discrimination, but with very low (raw) background, enters such region. Also shown for comparison are the contour lines of the other experiments, CDMS and EDELWEISS (thin dashed line), which have entered that region by using bolometers which also measure ionization. The DAMA region (closed line) corresponding to the $`3\sigma `$ annual modulation effect reported by that experiment and the exclusion plot obtained by DAMA NaI-0 (thin solid line) by using statistical pulse shape discrimination are also shown. A remark is in order: for CDMS two contour lines have been depicted according to a recent recommendation, the exclusion plot published in Ref.  (thin dotted line) and the CDMS expected sensitivity contour (thin dot-dashed line).

Data collection is currently in progress and some strategies are being considered to further reduce the low energy background. Another 50 % reduction from 4 keV to 10 keV (which could be reasonably expected) would allow to explore practically all the DAMA region in 1 kg y of exposure. In the case of reducing the background down to the flat level of 0.04 c/kg/keV/day (currently achieved by IGEX for energies beyond 20 keV), the DAMA region would be widely surpassed. In Figure 6 we plot the exclusions obtained with a flat background of 0.1 c/kg/keV/day (dot-dashed line) and of 0.04 c/kg/keV/day (solid line) down to the current 4 keV threshold, for an exposure of 1 kg year. As can be seen, the complete DAMA region (m=$`52_8^{+10}`$ GeV, $`\sigma ^p`$=($`7.2_{0.9}^{+0.4}`$)x10-9 nb) could be tested with a moderate improvement of the IGEX performances.

Table5

A new experimental project on WIMP detection using larger masses of Germanium of natural isotopic abundance (GEDEON, GErmanium DEtectors in ONe cryostat) is planned. It will use the technology developed for the IGEX experiment and it would consist of a set of $``$1 kg germanium crystals, of a total mass of about 28 kg, placed together in a compact structure inside one only cryostat. This approach could benefit from anticoincidences between crystals and a lower components/detector mass ratio to further reduce the background with respect to IGEX.

The GEDEON single cell is a cylindrical cryostat in electroformed copper (dimensions 20 cm diameter $`\times `$ 32 cm long) hosting 28 germanium crystals which share the same common copper cryostat (0.5 mm thick). The Ge crystals, are arranged in four plates suspended from copper rods. The cell is embedded into a precision-machined hole made in a Roman lead block providing a shield of 20 cm, and surrounded by another lead shielding 20 cm thick. A cosmic veto and a large neutron shield complete the shielding.

The preliminary MC estimated background in the $`150`$ keV region ranges from $`2\times 10^2`$ to $`2\times 10^3`$ c/keV kg day, according to the level of radioimpurities included as input. The radiopurity assays have been carried out in the Canfranc Laboratory for the lead and copper components of the shielding. The background final goal of GEDEON, below 100 keV, would be hopefully in the region of $`10^3`$ c/keV kg day and this value has been used to calculate anticipated $`\sigma `$(m) exclusion plots in the most favourable case. The expected threshold assumed has been $`\mathrm{E}_{\mathrm{thr}}=2`$ keV and the energy resolution in the low energy region has been taken $`\mathrm{\Gamma }1`$ keV. The MC estimated background of the GEDEON unit cell (28 crystals) is given in Fig. 7. A detailed study is in progress to assess the physics potential of this device. The exclusion plot which could be expected with such proviso for 24 kg y of exposure is shown in the Figure 6. Moreover, following the calculations presented in, GEDEON would be massive enough to search for the WIMP annual modulation effect and explore positively an important part of the WIMP parameter space including the DAMA region. A second phase of GEDEON with four cryostats (112 detectors and a total mass of 92 kg of Ge) is also being considered.

## Acknowledgments

The data presented here regarding Ge ionization detectors resulted from collaborative research with the members of the COSME and the IGEX collaborations, formed by C.E. Aalseth, F.T. Avignone, R.L. Brodzinski, S. Cebrián, E. García, W.K. Hensley, I.G. Irastorza, I.V. Kirpichnikov, A.A. Klimenko, H.S. Miley, J. Morales, A. Ortiz de Solórzano, S.B. Osetrov, V.S. Pogosov, J. Puimedón, J.H. Reeves, M.L. Sarsa, A.A. Smolnikov, A.G. Tamanyan, A.A. Vasenko, S.I. Vasiliev, J.A. Villar members to whom I am deeply indebted. The results contained in this paper have already been published in the open literature. I wish to thank especially S. Cebrián and I.G. Irastorza for their invaluable collaboration in the making of the exclusion plots and to J. Morales and I.G. Irastorza for many useful discussions on IGEX-DM. I thank S.Cebrian for allowing me the use of the MC results for GEDEON and J.Puimedon for the making of the radiopurity measurement. The present work was partially supported by the CICYT and MCyT (Spain) under grant number AEN99-1033.

Table4


format: text

text:
# The formation of Uranus and Neptune among Jupiter and Saturn

## 1 Introduction

The growth of Uranus and Neptune in the outer Solar System is not readily accounted for by conventional models of planet formation. A low primordial density of planetesimals and weak solar gravity would have made the process of accretion slow and inefficient. In direct N-body simulations of accretion among (approximately) Earth-mass bodies beyond 10 AU, performed with three different computer codes, little accretion is found to take place over timescales of $`10^8`$ years, and by extrapolation, over the age of the Solar System (Levison & Stewart 2001). Earlier simulations by Brunini & Fernandez (1999) showed accretion of the ice giants in several $`\times 10^7`$ years with the same initial conditions, but later simulations, performed with an improved integrator, require that bodies be enhanced in radius by at least a factor of ten, relative to bodies having the density of Uranus and Neptune, in order to recover the previous result (Brunini 2000). Therefore, Uranus and Neptune are unlikely to have formed from a late stage of mergers among large protoplanets, analogous to the putative final phase of planet formation in the terrestrial zone (eg. Wetherill 1996, Chambers & Wetherill 1998). The oligarchic growth model, in which the principal growth mode is accretion of small planetesimals by a protoplanet, also produces timescales which are too long (Kokubo and Ida 2000, Thommes 2000), unless the feeding zones of Uranus and Neptune can be replenished quickly enough with low random velocity planetesimals from elsewhere in the nebula (Bryden, Lin and Ida 2000).

Thommes, Duncan & Levison (1999), hereafter TDL99, develop an alternative model to in-situ formation for the origin of Uranus and Neptune. Beginning with four or more planetary embryos of 10-15 M in the Jupiter-Saturn region, they explore through N-body simulation the evolution of the system after one of these bodies (and in one case, two at the same time) accretes a massive gas envelope to become a gas giant. They find that the remaining giant protoplanets are predominantly scattered outward. Dynamical friction with the planetesimal disk subsequently recircularizes their orbits, which leads, in about half the simulations performed, to a configuration quite similar to the present outer Solar System, with the scattered giant protoplanets taking the roles of Uranus, Neptune and Saturn. These results suggest that Uranus and Neptune are actually potential gas giant cores which formed in the same region as Jupiter and Saturn, but lost the race to reach runaway gas accretion.

Here, we explore this model in more detail, and perform further simulations. Sections 2 and 3 motivate our choice of initial conditions for the simulations. Section 4 discusses the mechanism for transporting a proto-Uranus and -Neptune to the outer Solar System. N-body simulation results are presented in Section 5. The effect on the asteroid and Kuiper belts is discussed in Section 6. We summarize and discuss our findings in Section 7.

## 2 Available material in the Jupiter-Saturn region

Hayashi (1981) estimated the minimum primordial surface density of solids in the outer Solar System to be

$$\sigma _{min}(a)=2.7(a/5AU)^{3/2}g/cm^2.$$
(1)
The requirement that Jupiter and Saturn’s cores be massive enough to have initiated runaway gas accretion suggests that they are $``$ 10 M in mass (Mizuno et al 1978, Pollack et al 1996). Interior models are consistent with such a core mass, but also allow a coreless Jupiter (Guillot 1999). We assume in this work that both Jupiter and Saturn began as $``$ 10 M bodies; putting gas giant cores and ice giant cores on the same footing is necessary in the “strong” version of our model, though it is not essential to the basic mechanism; we discuss variations on the model in Section 7.

A surface density of 2.7 g/cm2 was likely too low to form a $``$ 10 M body at 5 AU before the gas was removed from the protoplanetary disk. Lissauer (1987) finds that a surface density of 15-30 g/cm2 is needed to allow formation of Jupiter’s core on a timescale of $`5\times 10^510^6`$ years, while the model of Pollack et al (1996), which includes concurrent accretion of solids and gas, produces Jupiter in less than $`10^7`$ years with 10 g/cm2. The formation of giant planet cores may have been triggered at least in part by the enhancement in the solids surface density beyond the “snow line”, where water goes from being a gaseous to a solid constituent of the protoplanetary disk. In fact, outward diffusion and subsequent freezing of water vapor from the inner Solar System may have resulted in a large local density enhancement around 5 AU, perhaps yielding a surface density even higher than 30 g/cm2 (Stevenson & Lunine 1988). Here we assume a power law surface density,

$$\sigma (a)=\sigma _0(a/5AU)^\alpha .$$
(2)
The above discussion suggests $`\sigma `$ 10 to 30 g/cm2 as a plausible surface density at 5 AU. Allowing, also, the exponent $`\alpha `$ to vary between 1 and 2, one obtains a total mass in the Jupiter-Saturn (J-S) region in excess of 40 M, and as high as 180 M. Therefore, it is likely that the region originally contained significantly more solids than ended up in the cores of Jupiter and Saturn.

## 3 Implications of oligarchic growth

Runaway growth (eg. Wetherill & Stewart 1989, Kokubo & Ida 1996) transitions to a slower, self-limiting mode called “oligarchic growth” (Ida & Makino 1993, Kokubo & Ida 1998, 2000) when the largest protoplanets are still orders of magnitude less than an Earth mass everywhere in the nebula. Oligarchic growth has previously only been demonstrated to take place interior to about 3 AU (Weidenschilling & Davis 2000). Though Kokubo and Ida make estimates of protoplanet mass and growth timescales in the giant planet region, they point out that their simulations are restricted to annuli which are narrow compared to their radii, and thus cannot make strong predictions about how oligarchic growth works over a wide range in semimajor axis (Kokubo and Ida 2000). However, their most recent simulations (Kokubo & Ida 2000b) span a range of 0.5 to 1.5 AU, and show oligarchic growth proceeding in an outward-expanding “wave” over time. Also, Thommes (2000) and Thommes, Duncan & Levison (2001) perform numerical simulations which suggest that oligarchic growth does take place in the outer Solar System, and proceeds on approximately the timescales predicted using the approach of Kokubo & Ida (2000).

Assuming that protoplanets remain approximately evenly spaced in Hill radii while growing—as is characteristic of oligarchic growth—their final mass is given by

$$M=(2\pi )^{3/2}\left(\frac{2}{3M_{}}\right)^{1/2}p^{3/2}n^{3/2}\sigma ^{3/2}a^3$$
(3)
where $`a`$ is the protoplanet semimajor axis, $`n`$ is the spacing between adjacent protoplanets in Hill radii $`r_H`$, and $`p`$ is the fraction of the total mass in the zone $`[anr_H/2,a+nr_H/2]`$ incorporated into the protoplanet (Kokubo and Ida 2000). Using numerical simulations, Kokubo & Ida (1998) show that $`n510`$. Models of giant planet formation by concurrent planetesimal and gas accretion suggest a surface density profile $`\sigma a^2`$ (Pollack et et 1996). Using this in Eq. 3, one obtains a protoplanet mass independent of semimajor axis. Adopting $`\sigma _0=10g/cm^2`$ and a spacing of n=7.5 Hill radii, one must then set the accreted mass fraction to about 0.8 in order to obtain a protoplanet mass of 10 M . With this spacing, between three and five such bodies fit between 5 and 10 AU. It is likely, therefore, that this region originally contained more than just the future solid cores of Jupiter and Saturn. At the same time, the more recent simulations of Thommes, Duncan and Levison (2001, preprint) indicate that accretion is quite inefficient in the outer Solar System; even in the Jupiter-Saturn region, the above value of $`p`$ is probably overly optimistic, and a correspondingly higher density of solids in this region may have been necessary to form gas giant cores. But as we find (Section 5.5), a higher-density disk actually tends to increase the “success rate” of the model.

## 4 Scattering and circularization of the protoplanets

In the nucleated instability picture of gas giant formation, there are three distinct phases of growth (Pollack et al 1996): In Phase 1 (which itself encompasses the sub-phases of runaway and oligarchic growth; see Section 3 above), a solid core grows until it depletes most of the material in its feeding zone. Phase 2 is characterized by much slower growth, with gas accretion gradually coming to dominate over solids accretion. After a time of order several Myrs, the protoplanet contains comparable masses of gas and solids; around this time, the third phase of runaway gas accretion sets in and proceeds on a timescale of $`10^5`$ years. One of the protoplanets in the Jupiter-Saturn region must have been the first to reach this point. Shorter formation timescales at smaller heliocentric radii argue for Jupiter having formed before Saturn. The relatively long time spent by a protoplanet on the “plateau” of Phase 2 means that several protoplanets could plausibly find themselves there by the time the first of them makes it to Phase 3. Furthermore, if the winner has a broad enough margin of victory, its rivals will have only had time to accumulate a few M of gas, thus resembling the present-day ice giants in both solids and gas mass.

A body increasing its mass from 10 M to Jupiter’s mass, MJ = 314 M, expands its Hill radius by a factor of about three. The adjacent protoplanets will therefore have a high likelihood of being gravitationally scattered. In a three-body system—Sun, Jupiter and a single protoplanet—a Jupiter-crossing protoplanet would continue to have close encounters with its much more massive scatterer, and so would remain coupled to it, unless it were scattered onto an unbound orbit and left the system altogether. In reality, however, one expects that as it crosses beyond the Jupiter-Saturn region, a body will encounter a less accretionally evolved part of the protoplanetary disk, consisting predominantly of much smaller bodies. As a result, the protoplanet will experience dynamical friction. This will tend to reduce the eccentricity of its orbit. If the eccentricity decays enough, the protoplanet’s perihelion will be lifted away from Jupiter, thus decoupling it from its scatterer. Insofar as one can neglect interactions with other scattered protoplanets, its eccentricity will then monotonically decrease until it reaches equilibrium with the planetesimals. The eccentricity rate of change due to dynamical friction on a body of mass $`M`$ in a swarm of mass $`m`$ bodies can be expressed as (Weidenschilling et al 1997)

$$\frac{de_M^2}{dt}=C(me_m^2Me_M^2)K_e$$
(4)
and a similar expression gives the time evolution of the inclination. $`e_m^2^{1/2}`$ is the RMS eccentricity of the mass $`m`$ bodies, and $`K_e`$ is a definite integral depending on the inclinations and eccentricities, defined in Stewart & Wetherill (1988). The coefficient C is given by

$$C=\frac{16G^2\rho _{sw}^{}{}_{1}{}^{}ln\mathrm{\Lambda }}{v_K^3(e_1^2+e_2^2)^{3/2}},$$
(5)
where $`G`$ is the gravitational constant, $`v_K`$ the local Keplerian velocity, $`\rho _{sw}^{}{}_{1}{}^{}`$ the spatial mass density of the swarm, and $`\mathrm{\Lambda }`$ is approximately the ratio of the maximum encounter distance between the body $`M`$ and a member of the population of the swarm, to the maximum separation that results in a physical encounter (Stewart & Wetherill 1988).

An equilibrium is therefore reached when

$$e_M^2Me_m^2m$$
(6)
This also means that as long as

$$e_M^2Me_m^2m,$$
(7)
the eccentricity decay rate of the body $`M`$ will only depend on the spatial mass density of the planetesimal disk, and will be essentially independent of the individual masses of the planetesimals. This feature will be exploited in the numerical simulations below.

## 5 N-body simulations

Putting the pieces together, one can now envision a scenario in which a) a rapidly growing Jupiter scatters its smaller neigbours outward, and b) these “failed cores” are decoupled from Jupiter and ultimately evolve onto circular, low-inclination orbits in the outer Solar System. TDL99 performed three series of numerical simulations to test this model; about half of the simulations produced a stable final configuration qualitatively similar to the present-day Solar System. Here, we expand on the earlier work, performing simulations to assess the effect of other planetesimal disk density profiles and of proto-Jupiter being initially not the innermost giant protoplanet. We also investigate how sensitive this scenario is to the timing of Saturn’s final stage of growth, relative to that of Jupiter. Simulations are once again performed using the proven SyMBA symplectic integrator (Duncan, Levison & Lee 1998). This integrator is able to handle close encounters among massive bodies while preserving the symplectic properties of the method of Wisdom & Holman (1991).

### 5.1 Initial conditions: Set 1

For the baseline set of simulations, a planetesimal disk of surface density

$$\sigma _1=10(a/5AU)^2g/cm^2$$
(8)
is used, with the disk extending from 5 to 60 AU. The above density profile is steeper than those of TDL99 ($`\sigma a^1`$ and $`a^{3/2}`$). The total disk mass between 10 and 60 AU (106 M) is about half that of the first two sets of runs in our earlier work (216 M), and similar to that in the third run (119 M). The truncation at 60 AU is to keep the number of bodies in the simulations tractably small; in reality the planetesimal disk may have extended for hundreds of AU, far beyond the presently visible Kuiper belt. Observational evidence does point to the possibility of a truncation of the belt at $``$ 50 AU (eg. Chiang & Brown 1999). But if a “Kuiper cliff” exists, it is not necessarily primordial.

Since we truncate the disk at a location where the planetesimal surface density is already very low (0.07 g/cm$`2`$ for the density profile of Eq. 8), one expects this not to have a strong effect on the evolution of any scattered protoplanets which cross the region $`r>60`$ AU. As a test, we perform four pairs of runs comparing the evolution of a 10 M body with initial $`a=100`$ AU, $`e=0.90.95`$ in planetesimal disks with a surface density as prescribed by Eq. 8, in one case extending to 60 AU, in the other to 200 AU. No systematic difference in the 10 M body’s semimajor axis evolution over the first few Myrs is apparent between the two different disk sizes. The eccentricity, however, tends to be damped more rapidly in a 200 AU disk. We expect, therefore, that the subsequent runs somewhat underestimate the effectiveness with which those scattered protoplanets which cross beyond 60 AU are circularized, if the disk was in fact larger. In particular, a planetesimal disk with a radius of hundreds of AU may allow the retention of protoplanets which in our runs are scattered strongly enough to become unbound; we will explore this possibility in future work.

The simulated disk is made up of equal-mass “planetesimals”, each having a mass of 0.2 M. At twice the mass of Mars, these bodies far exceed the actual characteristic mass of planetesimals in the early outer Solar System. In reality planetesimals are thought to have been on the order of 1 to 100 km in size, thus with a mass of $`10^{12}10^6`$ M (eg. Lissauer 1987). The unrealistically large masses are chosen to keep the number of bodies manageably low, at slightly over 500. As mentioned in Section 4, when the large bodies have eccentricities high enough that Eq. 7 is satisfied, the eccentricity decay timescale will be effectively independent of the planetesimal masses. Thus despite the large planetesimal masses, the initial strength of eccentricity damping will be realistic. The equilibrium eccentricity condition Eq. 6, however, does depend on the planetesimal mass, so the equilibrium eccentricity reached by a large body among the planetesimals in the simulation will tend to be unrealistically high. Of course, a true equilibrium between protoplanets and planetesimals will not exist anyway, since mutual perturbations among the protoplanets will also have an effect.

The initial planetesimal eccentricities and inclinations are given a Rayleigh distribution in $`e`$ and $`i`$, (Kokubo and Ida 1992) with $`e^2^{1/2}=0.05`$, $`i^2^{1/2}=0.025=1.4^o`$. In the numerical integrations, the planetesimals are treated as a “second-class”, non-self-interacting population. Thus they are perturbed in their Keplerian orbits only by forces from the protoplanets, not each other. The protoplanets, on the other hand, are subject to forces from each other as well as from the planetesimals. This serves two purposes: It makes the simulations run much faster, since for $`N`$ second-class bodies, the computation time scales as $`N`$ instead of $`N^2`$. Also, it prevents unrealistically strong self-stirring of the disk. Of course, not modeling planetesimal interactions means that collective planetesimal effects are not accounted for. Wave phenomena could have had an important effect on the evolution of the planetesimal disk velocity distribution (eg. Ward and Hahn 1998), provided the disk was sufficiently massive and dynamically cold. However, it is unlikely that significant wave phenomena could persist once the initial scattering has taken place and the planetesimal disk has been stirred by eccentric 10 M bodies.

Also unmodeled is nebular gas, either as a source of aerodynamic drag (relevant for small planetesimals; eg. Adachi, Hayashi and Nakazawa 1976), or as a source of tidal forces (relevant for bodies $``$ 0.1 M; eg. Ward 1997). The former mechanism will keep the planetesimal disk more dynamically cold, but this does not affect the strength of dynamical friction on larger bodies until $`e_M^2Me_m^2m`$. The latter effect is thought to cause the inward migration of $``$ 1 - 10 M objects on timescales short compared to their formation times. This of course constitutes a potentially severe problem not only for our model, but for any model of giant planet formation which requires the accumulation of large solid cores. Addressing this problem is beyond the scope of our present work, but we summarize some reasons why it may in reality have been less severe, in Section 7.

The 10 M bodies are given a density $`\rho =0.25\rho _{}=1.5`$ g/cm3, roughly equal to that of Uranus and Neptune. This gives them radii of $`2.18\times 10^4`$ km. Four such bodies are put in the simulation, initially on nearly circular and uninclined orbits. The orbits are spaced by 7.5 mutual Hill radii, starting from an innermost distance of 6 AU to allow for later inward migration by Jupiter. Thus the bodies’ initial semimajor axes are 6.0 AU, 7.4 AU, 9.0 AU and 11.1 AU. Between $`5`$ and 12 AU, the disk is depleted in planetesimals so that the surface density is still given by Eq. 8. Since the large bodies are spaced proportionally to their semimajor axes, their distribution is consistent with a surface density $`a^2`$.

Integrating the surface density from $`5`$ to 12 AU, the total mass is 51.5 M. With 40 M of this in the large bodies, this leaves 10.5 M in planetesimals in the Jupiter-Saturn region, in addition to 94.7 M between 12 AU, and the outer disk edge at 60 AU. In summary, then, the initial conditions amount to a state where Phase 2 of giant planet formation has been reached between 5 and 12 AU, with $`80\%`$ of the planetesimals having accreted through oligarchic growth into four bodies of 10 M each, while no large bodies have yet formed beyond this region (Fig. 2). The adoption of equal-mass protoplanets is a simplification; apart from stochastic variations among the oligarchic growth endproducts, one can expect some intermediate bodies, with masses perhaps up to a few M. Such bodies are likely to ultimately be cleared, along with the planetesimals, from the giant planet region. However, they may end up playing a role in the dynamics of the trans-Neptunian region; see Petit, Morbidelli and Valsecchi (1999).

For computational reasons, the inner simulation radius is chosen as 1 AU; any body whose orbit penetrates this boundary is eliminated from the system. This is done because a limitation of the SyMBA integrator used here is its inability to handle close perihelion passages. Although a new version of SyMBA has since been developed which removes this restriction (Levison and Duncan 2000), the runs presented here predate this development. In any case, this limitation has little relevance for the runs presented here; typically, less than ten planetesimals are lost at the inner boundary over the course of a 5 Myr simulation. The base timestep is chosen as 0.05 years, giving 20 steps per orbital period for an orbit with its semimajor axis at the inner radius, and over 200 steps per orbit for Jupiter. Experimentation shows that this timestep is small enough that the energy of the system is well-conserved. Runs initially go to 5 Myrs; in cases where the system still appears to be undergoing rapid evolution at this point, the runs are extended by another 5 Myrs.

### 5.2 Set 1 results

To model gas accretion, SyMBA was modified to allow a subset of bodies to have artificially time-varying masses (in addition to any changes in mass resulting from the accretion of other bodies). For the runs in Subset 1, it is assumed that the innermost protoplanet undergoes runaway gas accretion first, and grows into Jupiter. The simulations start at the time when this happens, which should be a few million years into the life of the solar system (see above). Runaway gas accretion is simulated by increasing the body’s mass over the first $`10^5`$ years of simulation time, from its original mass of 10 M to 314 M, approximately the present Jupiter mass. A linear growth in mass is used; this is deemed appropriate since the actual time evolution of mass during runaway growth is highly uncertain. Also, as the simulations will show, $`10^5`$ years is roughly the response time of the system, and the system’s subsequent evolution is therefore unlikely to be affected in a systematic way by the exact form of the time evolution of the runaway-phase mass growth.

Set 1 consists of eight alternate realizations of a run, differing only in the initial phases of the four protoplanets; for each, the angles $`\mathrm{\Omega }`$, $`\omega `$ and $`M`$ are randomly generated. As will be seen, the stochasticity of the system ensures that this difference in phases is sufficient to bring about a very different evolution in each of the versions of the run.

Fig. 3 shows the evolution of one of the eight runs, denoted as Run 1F, which after 5 Myrs produces final protoplanet orbits that bear a particularly close resemblance to those of the giant planets in the present-day Solar System. Semimajor axis versus time is plotted for each of the four protoplanets; the innermost one has grown into Jupiter after the first $`10^5`$ years. By this time, the protoplanet orbits begin to mutually cross, and strong scattering occurs. The protoplanet plotted in red briefly has its semimajor axis increased to greater than 100 AU. However, dynamical friction acts to reduce eccentricities universally, decoupling the protoplanets from Jupiter and from each other. By about 1.2 Myrs, none of the protoplanets are on crossing orbits anymore. After about 3 Myrs, the bodies no longer undergo any changes in semimajor axis greater than a few AU on a million-year timescale. At this point the orbits are well spaced and all eccentricities are $`0.05`$, with no large fluctuations. Due to the large stochastic variations among runs, the exact timescales differ, but we find that orbits typically become noncrossing after less than 5 Myrs in these and subsequent runs.

Subsequent semimajor axis evolution proceeds by scattering of planetesimals by protoplanets, rather than scattering of the protoplanets off each other and Jupiter. As planetesimals are scattered among Jupiter and the protoplanets, the former experiences a net loss of angular momentum while the latter experience a gain. Thus Jupiter’s orbit shrinks, while those of the protoplanets expand (Fernandez and Ip 1996, Hahn and Malhotra 1999). This phase takes place over a timescale of several tens of Myrs and ends, at the very latest, when the planetesimals have been cleared from among the planets. Because the length of migration during this phase is only a few AU, and to save time, most of the runs are stopped after 5 Myrs. As an example, Fig. 4 shows Run 1F continued to 50 Myrs; the net migration of the outer two protoplanets subsequent to 5 Myrs is only $``$ 1 - 4 AU outward.

The semimajor axes at which the scattered bodies end up are very noteworthy, if one compares them to the present orbits of the giant planets (Fig. 1). At 5 Myrs, the outer two protoplanets are at 16 and 26 AU. The innermost one is at 11 AU, while Jupiter is at 5.5 AU. This configuration of orbits is very similar to that of the present Solar System, where Uranus and Neptune are at 19 and 30 AU respectively, Saturn is at 9.6 AU, and Jupiter is at 5.2 AU. And at $`5\times 10^7`$ years, after some more net outward migration, the outer two protoplanets’ semimajor axes are even closer to those of Uranus and Neptune (though “proto-Saturn”, having also moved outward, is further away from its present orbit). Eccentricities and inclinations are likewise very close to their present values.

The end states of all the runs in Set 1 are summarized in Fig. 5. Depicted are snapshots of eccentricity versus semimajor axis at 5 Myrs, except for Runs 1A and 1H, which were continued to 10 Myrs. These were extended because one or more of the protoplanets still had a high but decreasing eccentricity at 5 Myrs.

Runs 1A, 1D, 1F, 1G and 1H result in a final ordering of orbits that is at the very least broadly consistent with the present Solar System: Jupiter is the innermost body, with the other three bodies interior to the region of the Kuiper belt, and eccentricities low enough that no protoplanets cross each other or Jupiter. Out of these five runs, 1D and 1F in particular resemble the present Solar System. Of course, to actually reproduce the Solar System, another important event has to take place: The next protoplanet beyond Jupiter must also undergo a runaway gas accretion phase to acquire an envelope of mass $``$ 80 M and become Saturn. With the innermost core on a stable orbit in the vicinity of Saturn’s present location, though—as it is in runs 1D and 1F—the time delay between Jupiter and Saturn’s runaway phases is not strongly constrained. A more detailed investigation of the role of Saturn’s growth will follow in Section 5.7.

In contrast, Uranus and Neptune must somehow have been prevented from later reaching runaway gas accretion. One possibility is that they simply ran out of time, still caught in the long plateau of the second giant planet growth phase (Pollack et al 1996; see Section 4) when the nebular gas was removed. Alternatively, it may be that the gas disk was truncated early on by photoevaporation between the orbit of Saturn, and the eventual orbit of Uranus (Shu, Johnstone and Hollenbach 1993). This possibility is discussed further in Section 7.

The final state of the planetesimal disk differs substantially among the runs. In 1B, 1C and 1E, the planetesimal disk is largely unperturbed over most of its radial extent, while in 1F and 1H, eccentricities have been greatly increased throughout the entire disk. 1A, 1D and 1G are intermediate cases. The extent of the perturbation simply depends on how much of the disk is crossed by the protoplanets over the course of the run. In Run 1F, for example, Fig. 3 shows that one body’s aphelion spends some time beyond 100 AU. On the other hand, in 1C, no protoplanet’s aphelion ever goes further out than $`Q_{max}`$ 18 AU. In those runs where $`Q_{max}<60`$ AU (the disk radius), $`Q_{max}`$ corresponds closely to where the planetesimal disk makes a transition from perturbed to unperturbed. The outer limit of the planetesimal disk’s excitation by scattered protoplanets has been referred to as the “fossilized scattered disk” (TDL99). The contemporary scattered disk, by contrast, consists of objects which were scattered by Neptune after the latter had attained its current orbit (Duncan and Levison 1997). A more detailed discussion of the fossilized scattered disk follows in Section 6.2.

Three of the runs produce systems at 5 Myrs that are irreconcilably different from our own. In 1B, one of the protoplanets has merged with Jupiter. In 1C, a protoplanet has been scattered onto an orbit interior to Jupiter, in the region of the present-day asteroid belt, with its semimajor axis at 3.4 AU, its perihelion at 2.8 AU and its aphelion at 3.9 AU. It attains a stable orbit not crossing Jupiter even though the only planetesimals available for damping interior to Jupiter are those few which are also scattered there. This is possible because the protoplanet is scattered not just by Jupiter, but also by the other protoplanets. Even assuming the protoplanet could be subsequently removed from this region, such an event would most likely have cleared much of the asteroid belt. Finally, in 1E, two of the protoplanets have merged. It should be noted that mergers—or, indeed, ejections—which reduce the number of protoplanets are not intrinsically a problem, since extra ones may have existed. Scenarios with five protoplanets (Jupiter + Saturn + 3) will be explored in Section 5.4.

### 5.3 Set 2: Dependence on initial ordering

How strongly does the final configuration of the system depend on the initial ordering of the protoplanets? The next set of runs, Set 2, uses the same initial conditions, within a random variation in the protoplanet phases, as Set 1. However, for these runs it is the second-innermost protoplanet, rather than the innermost one, which undergoes simulated runaway gas accretion. One can reasonably expect that this will favour an outcome like 1C (Fig. 5), where a protoplanet is scattered inward into the region of the asteroid belt.

The end states of the runs are shown in Fig. 6. In six of the eight runs, a protoplanet has indeed ended up interior to Jupiter. However, in two cases (2C and 2D), Jupiter is the innermost body. Thus it appears that if Jupiter does not grow from the innermost protoplanet, the likelihood of ending up with a final configuration similar to our Solar System declines, though such an outcome continues to be quite possible.

### 5.4 Set 3: Dependence on number of cores

How sensitively does the end state depend on the initial number of core-sized bodies? In the next set of simulations, an extra 10 M protoplanet is added. All protoplanets are more tightly spaced, by 6.5 instead of 7.5 mutual Hill radii. Starting, again, from 6.0 AU, the outermost protoplanet is therefore initially at 12.2 AU. The surface density of planetesimals in the region of the protoplanets is reduced to keep the average surface density unchanged at $`\sigma _1=10(a/5\text{AU})^2`$ g/cm2. The innermost protoplanet is, again, increased in mass to that of Jupiter over the first $`10^5`$ years.

Fig. 7 shows the endstates of the runs. This time all runs initially have a length of $`10^7`$ years, since the larger number of bodies take longer to decouple from each other. Run 3A is continued to $`1.5\times 10^7`$ years, because after $`10^7`$ years some of the protoplanets are still on crossing orbits. Eccentricities are uniformly low and the protoplanet orbits are well-spaced for the most part, thus the systems have a good chance of being stable indefinitely. However, once all the planetesimals have been scattered from among the protoplanets, so that the latter are no longer subject to dissipative forces, some of these systems may still become unstable. This caveat applies to all the runs presented in this work, but generally speaking, larger numbers of bodies increase the potential for instability (Levison, Lissauer and Duncan 1998).

All of the protoplanets remain in six of the eight runs, thus resulting in systems with one too many planets relative to the present Solar System. However, in Runs 3B and 3F, one of the protoplanets is ejected from the Solar System, leaving the right number of bodies behind. In both of these runs, a protoplanet ends up with a semimajor axis within $`10\%`$ of present-day Uranus and Neptune, respectively, though both “Saturns” are too far out. One may conclude that with one extra initial protoplanet in the Jupiter-Saturn region, scattered protoplanets continue to be readily circularized, and the resulting systems tend to look like ours with one extra outer planet. However, a system with four giant planets remains a possible outcome. Also, the subsequent formation of Saturn from one of the protoplanets may trigger more ejections, especially in those cases where the inner protoplanets are more closely spaced, such as 3D and 3H.

### 5.5 Set 4: A more massive planetesimal disk

In this set of runs, a number of parameters are changed to simulate a system with a more massive planetesimals disk. The protoplanets are now 15 M bodies. The planetesimal disk surface density profile is still $`a^2`$, but it is scaled up to be 15 g/cm2 at 5 AU; that is,

$$\sigma _2=15(a/5\text{ AU})^2g/cm^2.$$
(9)
Other minor differences are a legacy of chronologically earlier runs. The innermost protoplanet is initially at 5.3 AU, and successive protoplanets are spaced by only 5.8 mutual Hill radii. Thus the outermost protoplanet is initially at 9.0 AU. The individual planetesimals have a mass of 0.24 M.

The innermost protoplanet, as before, has its mass increased to 314 M over the first $`10^5`$ simulation years. The end states of the runs are shown in Fig. 8. All except C yield the correct number and ordering of bodies, and eccentricities are uniformly low. This “success rate” is higher than that of Set 1, in which only five out of eight runs yield qualitatively the correct orbital configuration. This is accounted for by the more massive planetesimal disk; it provides stronger dynamical friction, so that scatterings of protoplanets tend to be less violent, and subsequently, orbits tend to be circularized and mutually decoupled more quickly. Runs A, D and H end up with protoplanet orbits that are particularly close to those of Saturn, Uranus and Neptune. Jupiter systematically ends up at too small a heliocentric distance, indicating that the initial distance of 5.3 AU is too small. Also, the larger disk mass gives the protoplanets and Jupiter more planetesimals to scatter, and thus increases the distance they travel due to angular momentum exchange (Section 5.2).

### 5.6 Set 5: A shallower disk density profile

In this set of runs, a shallower planetesimal disk surface density, $`a^{3/2}`$, is used. The disk now begins at 10 AU, with the interior region being initially occupied solely by the protoplanets. In other words, it is assumed that in the epoch from which the runs start, all but a negligible mass of the planetesimals among the protoplanets has been swept up or scattered from the region. The surface density profile is given by

$$\sigma _3=1.8(a/10\text{ AU})^{3/2}g/cm^2$$
(10)
and thus a total mass of 123 M is contained in the disk between 10 and 60 AU. Extending this planetesimal disk inward to 5 AU would yield a surface density there of only 5 g/cm2, and a total mass between 5 and 10 AU of only 25 M. Despite this, we still put four 15 M in this region; thus, it is assumed that the original planetesimal surface density profile in this region was steeper, perhaps due in part to redistribution of water vapor from other parts of the disk to the vicinity of the snow line (Stevenson & Lunine 1988).

This set consists of twelve runs, each to 5 Myrs. As before, the inner body’s mass is increased to 314 M over the first $`10^5`$ years of simulation time. The endstates are shown in Fig. 9. Eight out of twelve runs possess the right number and ordering of bodies. Eccentricities of the protoplanets and Jupiter are $`0.1`$ in five of these. This “success rate” is close to that of Set 1 (which has a similar total disk mass); this model thus does not appear to be highly sensitive to changes in density profile alone.

Various degrees of disruption of the planetesimal disk can again be seen. Most show a sharp transition between a disrupted region crossed by the scattered planetesimals, and a largely undisturbed outer region. In Run 5A, for example, this transition occurs at slightly below 40AU, while in Run 5F, it is located between 45 and 50 AU. Runs 5B and 5L show strong disruption throughout the entire disk, indicating that all of it was crossed by one or more protoplanets. Another state can be seen in Run 5G, where all eccentricities in the outer part of the belt are uniformly raised. This occurs when a protoplanet crosses the outer disk with an inclination high enough that it spends most of its orbit above or below, rather than inside, the disk. The disk planetesimals are then excited primarily by long-range secular effects rather than by short-range scattering encounters (TDL99). Another example of this effect can be seen in Run 1G above.

### 5.7 Set 6: The role of Saturn

Thus far, the only gas giant in the simulations has been Jupiter. The innermost of the scattered protoplanets does tend to end up near the present location of Saturn. However, to reproduce the Solar System, it must at some point accrete $``$ 80 M of nebular gas. We therefore investigate what effect the subsequent growth of a Saturn-mass object has on our model.

A Set 1 run which produced a good Solar System analogue (1F; Fig. 3) is used as a starting point. At about 1 Myr, the protoplanets and Jupiter are no longer on crossing orbits. At this time, the initially outermost protoplanet (plotted in blue) has become the innermost one, closest to Jupiter at $``$ 10 AU. We perform a set of five runs which branch off from this point. In these runs, the innermost protoplanet has its mass increased to that of Saturn over a $`10^5`$ year interval, starting at 1, 1.2, 1.4, 1.5 and 1.6 Myrs, respectively. We choose to start only after the protoplanets are on noncrossing orbits in order to avoid ending up with an eccentric “Saturn”; the eccentricity evolution of an initially eccentric giant planet in a gas disk is uncertain (Lin et al 2000).

The endstates are shown in Fig. 10.

No protoplanets have been lost from the system by the end of the runs. A protoplanet still has a high eccentricity in the 1 Myr case; this is because this protoplanet’s perihelion is still very close to the innermost protoplanet’s orbit at 1 Myr, and it suffers strong perturbations as the latter grows to Saturn’s mass. In the other cases, however, Saturn’s formation does not cause large eccentricities in the protoplanets. This is as one would expect; at its final mass, and at 10 AU, Saturn’s Hill radius is 0.45 AU, and by $`1.2\times 10^6`$ years, the closest protoplanet’s perihelion is at $``$ 14 AU, almost $`9r_H`$ away, thus unlikely to be in reach of strong scattering. We leave for future work the effect of Saturn’s gas accretion on systems with more protoplanets, such as those in Set 3 in Section 5.4 above. Such systems tend to have weaker stability, and thus may be more susceptible to disruption by Saturn’s final growth spurt.

One effect visible in Fig. 10 is that the semimajor axis of Saturn at 5 Myrs tends to be smaller than that of the innermost scattered protoplanet—the putative proto-Saturn—in those runs where the gas accretion of Saturn is not modeled (see for example Fig. 5). When a protoplanet grows to Saturn’s mass, its subsequent migration speed is much slower, since the rate of migration depends on how much mass in planetesimals it scatters relative to its own mass. This counteracts the tendency of the innermost scattered protoplanet to end up at a semimajor axis larger than that of present-day Saturn in the other runs, where only Jupiter grows.

## 6 Scattered protoplanets and the small body belts

The simulations presented above show that scattering of giant planet core-sized protoplanets is a violent event, which leaves a strong dynamical signature on the surrounding planetesimal disk. The asteroid belt and the Kuiper belt are therefore the natural places to look for evidence of large scattering events in the Solar System’s early history.

### 6.1 The asteroid belt

It is those members of the asteroid belt larger than about 50 km in diameter which are of interest in inferring properties of the early Solar System; smaller bodies cannot be primordial because they could not have survived intact for the age of the Solar System (eg. Petit, Morbidelli and Valsecchi 1999). This population displays a degree of dynamical excitation that is not readily explained. The other puzzling feature of the asteroid belt is its severe mass depletion—by at the very least a factor of $`10^3`$—relative to the amount of mass the region originally contained, extrapolated from the terrestrial zone (Weidenschilling 1977).

In the context of the model presented here for the early evolution of Uranus and Neptune, the initial violent scattering of protoplanets is perhaps the most obvious candidate to look to for perturbation of the asteroid belt. Indeed, in numerous runs, protoplanets spend some time interior to the orbit of Jupiter, crossing part of the asteroid belt region for up to $`10^4`$ years at a time. One might expect that such an occurence would wreak havoc in the the asteroid belt. To explore this possibility, we perform six additional runs with planetesimals added in the asteroid belt region. The individual bodies have a mass of $`0.024M_{}`$, and are distributed with a surface density of

$$\sigma _{belt}=8.0(a/1\text{ AU})^1g/cm^2$$
(11)
between 2.5 and $`4.5AU`$. This shallow density profile is the same as the one used by Chambers and Wetherill (1998) in the terrestrial region, which in turn was chosen to be more consistent with the large densities required at larger heliocentric distances to form Jupiter and Saturn. The protoplanet masses in this case are $`15M_{}`$. We find that eccentricities can get excited to their present values in this way, though only down to the crossing protoplanet’s minimum perihelion distance, which in none of the runs performed reaches the inner edge of the belt. Inclinations fare more poorly; protoplanets seldom raise them much above $`10\mathrm{°}`$, which is less than the present median inclination of asteroids beyond 2.5 AU. Also, very little mass is scattered out of the belt while the protoplanets are crossing it. Fig. 11 shows “before and after” snapshots of the run which, out of the six, displays the strongest disruption of the region interior to Jupiter by scattered protoplanets. As can be seen, few bodies, apart from those that are nearly Jupiter-crossing, attain inclinations above $`15\mathrm{°}`$.

Scattered protoplanets may have had a more indirect role in the excitation and mass depletion of the asteroid belt. If the majority of the solids in the asteroid belt accreted to form a system of planet-sized bodies (raising eccentricities and inclinations of the remaining planetesimals in the process), scattered protoplanets crossing into the region may have contributed to making this system unstable, akin to the model of Chambers (1999). Also, if proto-Saturn was still sufficiently close to Jupiter when it accreted its massive gas envelope, the subsequent migration of the gas giants may have been enough to sweep the inclination-exciting $`\nu _{16}`$ secular resonance through most of the asteroid belt (cf. Gomes 1997, Levison et al 2001) As an example, in one of the runs from Set 6, Saturn commences growing at $`4\times 10^5`$ years and reaches its final mass at $`5\times 10^5`$ years, at which point Jupiter and Saturn are at 5.7 and 8.4 AU, respectively. This places the $`\nu _{16}`$ resonance at $`3.5`$ AU, and it moves inward as Jupiter and Saturn move apart. In contrast, Gomes (1997), using a more moderate range of migration for Jupiter and Saturn, finds that only the region inward of 2.7 AU is crossed by the $`\nu _{16}`$ resonance.

### 6.2 The Kuiper belt

In the present Solar System, a new class of Kuiper belt object (KBO) has recently been identified (Duncan and Levison 1997, Luu et al 1997). These objects have semimajor axes and eccentricities such that they lie near the locus of Neptune-encountering objects shown in Fig. 1. They are thought to be part of a population referred to collectively as the scattered disk—formerly low-eccentricity KBOs which have had their orbits changed by close encounters with Neptune. Many of the simulations in Section 5 show an analogous class of planetesimals in their “Kuiper belt” regions. However, these fall on the locus of orbits crossing not the final semimajor axis of the outermost protoplanet, but the furthest aphelion distance of any of the protoplanets during their initial high-eccentricity phase. Since these orbits are no longer being crossed by a protoplanet, they will be stable over long times. One can refer to these structures as “fossilized” scattered disks (TDL99), because they preserve part of the dynamical history of the planetesimal disk. Such structures only appear in runs where the initial scattering was strong enough that one or more protoplanets had their aphelia increased to well beyond the final semimajor axis of the (ultimately) outermost protoplanet.

Observations of our Solar System’s Kuiper belt do indeed reveal an anomalously high degree of excitation (eg. Petit, Morbidelli & Valsecchi 1999, Malhotra, Duncan & Levison 2000). The eccentricities and, to a lesser degree, inclinations of bodies in mean-motion resonances with Neptune, particularly the 2:3 resonance at 39.5 AU, can be explained by resonance sweeping during Neptune’s migration, as can the paucity of objects on nonresonant orbits interior to 39 AU (Malhotra 1995). However, the high inclinations found beyond $``$ 41 AU in what is commonly called the “classical” Kuiper belt, cannot be explained in this way. Petit, Morbidelli and Valsecchi (1999) propose large (up to 1 M) Neptune-scattered planetesimals as the mechanism which stirred and cleared the belt. But even when such bodies remain in the belt for 100 Myrs, the inclinations they raise are almost always less than $`20\mathrm{°}`$.

Can the excitation of the Kuiper belt be better accounted for if we are thus far only seeing the part of it that is interior to the locus of a fossilized scattered disk? The inclinations raised by a protoplanet can be directly obtained from the simulations of Section 5. The runs of Set 1, Set 3 and Set 5 will be used for comparison. Inclinations for Sets 1, 3 and 5 are shown in Figs. 12, 13 and 14, respectively. Inclinations beyond the outermost protoplanet are excited up to a maximum of about $`30\mathrm{°}`$, similar to those observed in the classical Kuiper belt today (see Fig. 1) However, observations of the Kuiper belt are biased against high-inclination objects. Brown (2001) derives a de-biased inclination distribution function for the classical belt:

$$f_t(i)=\mathrm{sin}(i)\left[a\mathrm{exp}\left(\frac{i^2}{2\sigma _1^2}\right)+(1a)\mathrm{exp}\left(\frac{i^2}{2\sigma _2^2}\right)\right]$$
(12)
with $`a=0.93\pm 0.02`$, $`\sigma _1=2.2\pm _{.6}^{.2}`$, and $`\sigma _2=17\pm 3`$. One can define a parameter $`i^{}\mathrm{cos}^1(\overline{\mathrm{cos}(i)})`$ to give a measure of the characteristic inclination of a population. For the de-biased distribution above, $`i^{}=21\mathrm{°}`$. However, in the runs presented here, the largest $`i^{}`$ in the region corresponding to the classical belt (between the outermost planet’s 2:3 and 1:2 mean-motion resonance) is only $`15\mathrm{°}`$. Thus, although higher inclinations are produced here than in the large Neptune-scattered planetesimals model of Petit, Morbidelli and Valsecchi (1999), the inferred full velocity distribution of the classical Kuiper belt still cannot be accounted for.

It is possible to estimate with a simple numerical experiment if a planet as large as Uranus or Neptune can in principal excite the Kuiper belt to observed values. This experiment consists of a single Uranus-mass planet on an orbit with $`a=45AU`$, $`e=0.2`$, and $`i=25\mathrm{°}`$. The planet is embedded in a swarm of 400 massless test particles informally spread from 35AU to 55AU, with initial $`e=0.01`$ and $`i=1\mathrm{°}`$. Fig. 15 plots the $`i^{}`$ of the particles as a result of scattering off of the planet. It shows that a planet the mass of the ice giants can indeed excite the Kuiper belt to $`i^{}=20^{}`$ in a million years.

However, in none of the runs we performed was such a high inclination imparted on a protoplanet. Alternatively, a less inclined protoplanet may be able to reproduce the Kuiper belt inclinations if it remains on a belt-crossing orbit significantly longer than 1 Myr (and thus longer than in our runs). This may require a less massive trans-Saturnian planetesimal disk, in order to increase the dynamical friction timescale, though that in turn would increase the chances of protoplanets being scattered out of the Solar System altogether. A longer circularization time will also result if a protoplanet can be decoupled from Jupiter and the other protoplanets at a larger semimajor axis, so that the planetesimal surface density averaged over its orbit is lower. For the latter situation to have a better chance of occurring, the planetesimal disk needs to be extended to larger heliocentric distance. We will investigate this issue further in future work.

## 7 Discussion

The conventional picture of Uranus and Neptune’s formation, whereby the ice giants accrete near their present heliocentric distances, has grave problems. Numerical simulations have not been able to produce $``$ 10 M objects in the trans-Saturnian region in the lifetime of the Solar System without significantly increasing protoplanet radii, or invoking dissipational forces and planetesimal disk densities too large to be consistent with a physically plausible protostellar disk.

Building on our previous work (TDL99), we have performed additional N-body simulations of the evolution of the outer Solar System starting at the time when the first gas giant forms. At this point, we assume that a number of $``$ 10 M objects have formed at a heliocentric distance of roughly 5 to 10 AU, as is suggested by the oligarchic growth model (Kokubo & Ida 1998, 2000). Using a variety of different initial conditions, we find as before that the accretion of Jupiter’s gas envelope causes the remaining protoplanets to become violently unstable. In most cases they are scattered onto high-eccentricity orbits in the trans-Saturnian region. With most of its orbit now crossing the largely pristine trans-Saturnian planetesimal disk, a scattered protoplanet experiences dynamical friction and has its eccentricity rapidly damped. As a result, the protoplanets tend to end up on nearly circular, well-spaced orbits on a Myr timescale, with semimajor axes comparable to those of Saturn, Uranus and Neptune. Of the simulations which initially contain a total of four giant protoplanets and form Jupiter from the innermost, the majority produce final orbital configurations similar to that of our outer Solar System. Such systems are produced—though with lower probability—even if ones adds an extra protoplanet, or lets the second-innermost protoplanet grow into Jupiter. These results strengthen our earlier conclusion that if Uranus and Neptune shared the same birthplace as the gas giants, they could then readily have been delivered to their present orbits.

The role of migration in the formation of Uranus and Neptune was previously investigated numerically by Ipatov (1991), based on an idea by Zharkov and Kozenko (1990). Ipatov concludes that planetary embryos of a few M may have originated just outside the orbit of Saturn, to migrate outward and later grow into Uranus and Neptune, provided that they did not acquire high eccentricities during this process. In contrast, we find that Uranus and Neptune could have originated from anywhere in the Jupiter-Saturn region, and that initially high eccentricities—which nearby bodies naturally tend to acquire during Jupiter’s final growth phase—are in fact a powerful mechanism for rapidly transporting them outward. Also, the long growth timescales in the outer Solar System suggest that Uranus and Neptune likely already completed most of their growth in the Jupiter-Saturn region; even with a “head start” of a few Earth masses, the formation of Uranus- and Neptune-mass objects much beyond 10 AU within the age of the Solar System is unlikely (Levison and Stewart 2001, Thommes, Duncan and Levison 2001)

Can one find any observational support for this model in the present-day Solar System? The high inclinations in the classical Kuiper belt point to strong dynamical excitation in the past, and the simulations performed here do produce high inclinations in this region as a natural side effect. However, the simulations all fall short of reproducing the mean debiased inclination of the classical Kuiper belt. Strong observational support would be provided by the discovery of a fossilized scattered disk in the Kuiper belt, and a dynamically colder population beyond. It is tempting to link the trans-Neptunian object 2000 CR105 with a fossilized scattered disk; its high eccentricity (0.8) is characteristic of a scattered disk object, but recent observations (Gladman et al 2001) have established that its perihelion is at 44 AU, far beyond the reach of Neptune. However, none of the fossilized disk objects in our simulations acquire semimajor axes as high as that of 2000 CR105 (216 AU). Fig. 16 shows a plot of perihelion distance versus semimajor axis for Set 5, revealing only one case (B) where one or more objects simultaneously acquire a semimajor axis of $``$ 100 AU and a perihelion distance significantly further out than the (circularized) outermost protoplanet. All other sets of runs fare even more poorly. 5B is a run in which a protoplanet spends a long time at high eccentricity, and excites particularly high planetesimal inclinations in the disk. A long circularization time may thus be an important ingredient in reproducing both objects like 2000 CR105 and the high inclinations of the Kuiper belt; we will address this possibility in future work.

Findings regarding the deuterium to hydrogen (D/H) ratios of Uranus, Neptune and comets are particularly interesting in the context of this model. From infrared observations, Feuchtgruber et al (1999) find that the D/H ratios of the ice giants are lower than those of the Oort cloud comets Halley (Eberhardt et al 1995), Hyakutake (Bockelée-Morvan et al. 1998) and Hale-Bopp (Meier et al 1998) by a factor of approximately three, a difference large compared to the uncertainties of the measurements. Oort cloud comets are thought to have originated primarily in what is today the Uranus-Neptune region (Duncan, Quinn & Tremaine 1987, Fernandez 1997), and the notion of a common birthplace is supported by the comets’ similar D/H ratios. Though a sample size of three is clearly very small, this discrepancy between the comets and the ice giants would seem to present a further problem for any scenario in which Uranus and Neptune form in the trans-Saturnian region, since they should then share the chemical composition of the comets. However, this is exactly what one would expect if the ice giants originally formed at a smaller heliocentric distance, where higher temperatures would have made for a lower D/H ratio.

An aspect not modeled in any of the simulations is the gravitational interaction of the bodies with a gaseous disk. It has been shown that gas disk tidal forces can cause rapid inward migration of protoplanets (eg. Ward 1997). In fact, the speed of migration may be peaked for objects of $``$ 10 M, taking place on timescales of $`10^5`$ years or less. This peak corresponds to the transition between so-called Type I migration, where a body’s resonant interaction with the gas disk gives rise to a torque imbalance, to Type II migration, where the object opens a gap in the gas disk and is subequently locked to the disk’s viscous evolution. Tidal migration therefore poses a problem for any model of giant planet formation: how do they avoid spiralling into the central star as they form? However, tidal torques may not in fact operate throughout the whole disk. Gammie (1996) develops an accretion disk model which, beyond $``$ 0.1 AU, only transports angular momentum in a relatively thin surface layer. Thus the bulk of the disk would have zero viscosity, and objects embedded in it would be subject to neither Type I nor Type II orbital decay. Also, the disk may be truncated early on by photoevaporation from the central star (Shu, Johnstone and Hollenbach 1993), or from surrounding stars, as is seen in the Orion Nebula proplyds (Johnstone, Hollenbach and Bally 1998). Bodies scattered beyond the truncation distance would then be safe from nebula tides, and those near the edge would experience a net positive torque and migrate outward. Of course, such a scenario has the added advantage that, as discussed in the above works, it accounts for Uranus and Neptune having no massive gas envelopes.

In constructing the simulations presented here, we have appealed to the oligarchic growth model as a plausible guide for our choice of initial protoplanet masses and orbits. One can of course envision a number of variations. For instance, two gas giants may form on initially widely separated orbits. The waves they raise in the disk will tend to clear the gas between them, and Type II migration will cause their orbital separation to decrease (Kley 2000). If ice giant sized bodies are trapped in between, they will be prevented from accreting more gas, and will be scattered as the stable region between the gas giants shrinks to nothing. Also, if future measurements of Jupiter constrain its core to be much smaller, or even absent, this will invalidate the assumption that Jupiter grew from an ice giant sized nucleus. Though the scattering of ice giants could still take place in principal, one would then need to explain how a small body won the gas accretion race against much larger bodies (if the core is small), or how ice giant sized bodies managed to form before the birth of the first gas giant from an unnucleated disk instability (if there is no core; eg. Boss 2000). Alternatively, one could search for a way in which accretion could continue to take place in close proximity to a gas giant, with scattering delayed until ice giant sized bodies form. For such a scenario Type I tidal effects, which cause the eccentricity to decay on an even shorter timescale than the semimajor axis (eg. Papaloizou & Larwood 2000), may actually be helpful. It may then be the dispersal of the gas which triggers scattering. Doing away with any reference to a specific formation process, the most general statement of our results is: Ice giant sized bodies can be scattered from the Jupiter-Saturn region by gas giant sized bodies, to ultimately end up on Uranus- and Neptune-like orbits.

EWT is grateful for support from the Center for Integrative Planetary Science, as well as from the National Sciences and Engineering Research Council (NSERC) during the earlier part of this work. MJD is grateful for support from NSERC. HFL is grateful for support from NASA’s Planetary Geology & Geophysics, Origins of Solar Systems, and Exobiology programs.