text:
# 1 Introduction

## 1 Introduction

More than fifty candidates for extrasolar planets have been announced over the past six years (e.g., Schneider 2001). In each case, precise stellar radial-velocity measurements indicated the presence of a low-mass unseen companion, with a minimum mass between 1 and about 10 Jupiter masses ($`M_{Jup}`$). The identification of these unseen companions as planets relied on their masses being in the planetary range.

However, the actual masses of the planet candidates are not known. The radial-velocity data yield only $`M_2\mathrm{sin}i`$, where $`M_2`$ is the secondary mass and $`i`$ is the inclination angle of its orbital plane, which cannot be derived from the spectroscopic data. Nevertheless, the astronomical community considered the planet-candidate masses as being close to their derived minimum masses — $`M_2\mathrm{sin}i`$. This is so because at random orientation the most probable inclination is $`90^{}`$, and the expected value of $`\mathrm{sin}i`$ is close to unity.

Very recently some doubt has been cast about the validity of the random orientation assumption. Gatewood, Han, & Black (2001) and Han, Black, & Gatewood (2001) analysed the Hipparcos astrometric data of the stars hosting planet candidates together with the stellar precise radial-velocity measurements and derived in some cases very low inclination angles for the orbital planes. Han, Black, & Gatewood (2001) found eight out of 30 systems with an inclination smaller or equal to $`0.5^{}`$, four of which they categorized as highly significant. The probability of finding such small inclinations in a sample of orbits that are isotropically oriented in space is extremely small, indicating either a problematic derivation of the astrometric orbit, or, as suggested by Han, Black, & Gatewood (2001), some serious orientation bias in the inclination distribution of the sample of detected planet candidates.

However, the analysis of the Hipparcos data can be misleading. As has been shown by Halbwachs et al. (2000), one can derive a small false orbit with the size of the typical positional error of Hipparcos, about 1 milli-arc-second (=$`mas`$), caused by the scatter of the individual measurements. Therefore, one should carefully evaluate the statistical significance of any astrometric orbit of that size derived from the Hipparcos data. In Section 2 we summarize our work (Zucker & Mazeh 2001a) that evaluates the significance of the astrometric orbits by applying a permutation test to the Hipparcos data. Similarly to the results of Pourbaix (2001) and Pourbaix & Arenou (2001), we also find that the significance of all the Hipparcos astrometric orbits of the planet candidates are less than 99%, including $`\rho `$ CrB that attracted much attention after the publication of Gatewood, Han, & Black (2001) suggestion. We therefore conclude that the Hipparcos data does not prove the anisotropy of the orientations of the orbital planes of the planet candidates.

After showing that the random orientation in space is still a reasonable assumption, not confronted by any available measurement, we present in Section 3 our work (Zucker & Mazeh 2001b) that uses this assumption to derive the mass distribution of the planet candidates. This is done with MAXLIMA, a MAXimum LIkelihood MAss algorithm which we constructed to derive the mass distribution. Similar to the results of Jorrisen, Mayor & Udry (2001), we show that the mass distribution of the planet candidates is separated from the one of the secondary masses by the so-called ’brown-dwarf desert’ (e.g., Marcy & Butler 2000). This indicates that we are dealing with two different classes of objects. One is the giant planets, with masses not far from the planetary mass range, while the other is the low-mass secondaries, with stellar mass range.

One could speculate that the separation between the two different mass distributions indicates different formation processes. The commonly accepted paradigm is that planets were probably formed by coagulation of smaller, possibly rocky, bodies, whereas stars were probably formed by some kind of fragmentation of larger bodies. In other words, planets were formed by small bodies that grew larger, whereas stars, binary included, were formed by fragmentation of large bodies into smaller objects (e.g., Lissauer 1993; Black 1995). This could imply, for example, that the distribution of orbital eccentricities of giant planets and low-mass binaries would be substantially different. All the solar planets have nearly circular orbits, whereas binaries have eccentric orbits (e.g., Mazeh, Mayor, & Latham 1996). We could also expect the periods of planets to be longer than 10 years, like the giant planets in the solar system. Many studies of the newly discovered planets showed that this is not the case (e.g., Marcy, Cochran, & Mayor 2000). Moreover, following Heacox (1999) who based his analysis upon only 15 binaries and a handful of planet candidates, we show in Section 4 that within some reasonable restrictions, the eccentricity and period distributions of the two samples are surprisingly similar. Similar results have been obtained by Stepinski & Black (2001a,b,c). In Section 5 we consider a toy model that can generate the eccentricity distribution of both populations.

## 2 The Significance of the Astrometric Orbits

In this section we present our work (Zucker & Mazeh 2001a) where we evaluate for each of the extrasolar planets the statistical significance of its astrometric orbit, derived from the Hipparcos data together with its radial-velocity measurements. We first derived the best-fit orbit by assuming that the spectroscopic and astrometric solutions have in common the following elements: the period, $`P`$, the time of periastron passage, the eccentricity, $`e`$, and the longitude of the periastron. In addition, the spectroscopic elements include the radial-velocity amplitude, $`K`$, and the center-of-mass radial velocity. We have three additional astrometric elements — the angular semi-major axis of the photocenter, $`a_0`$, the inclination, $`i`$, and the longitude of the nodes. In addition, the astrometric solution includes the regular astrometric parameters — the parallax, the position and the proper motion.

In most cases the elements are not all independent. From the spectroscopic elements we can derive the projected semi-major axis of the primary orbit. This element, together with the inclination $`i`$ and the parallax, yields the angular semi-major axis of the primary, $`a_1`$. Assuming the secondary contribution to the total light of the system is negligible, this is equal to the observed $`a_0`$.

To find the statistical significance of the derived astrometric orbit in each case we applied a permutation test (e.g., Good 1994) to the Hipparcos data. For each star we generated simulated permuted astrometric data and analyzed them either together with the actual individual radial velocities of that star, or by imposing the published spectroscopic elements. Details of the analysis are given by Zucker & Mazeh (2001a).

The distribution of the falsely detected semi-major axes indicated the range of possible false detections. For example, $`a_{99}`$ — the 99-th percentile, denotes the semi-major axis size for which 99% of the simulations yielded smaller values. Consequently, an astrometric orbit is detected with a significance of 99% if and only if the actually derived semi-major axis, $`a_{derived}`$, is larger than $`a_{99}`$.

As an illustration, Figure 1A shows the histogram of the semi-major axis derived by random permutations of the Hipparcos data of HD 209458. This star’s inclination is known to be close to $`90^{}`$ through the combination of radial velocity and transit measurements (Charbonneau et al. 2000; Mazeh et al. 2000; Henry et al. 2000; Brown et al. 2001). The Hipparcos derived semi-major axis, $`a_{derived}`$, is 1.76 $`mas`$, which is marked in the figure by an arrow. One can clearly see that many random permutations led to larger semi-major axes, a fact that renders this derived value insignificant. The derived value is obviously false since the known inclination implies a value of less than a micro-arc-second.

In Figure 1B we show an opposite case, HD 164427, where the derived astrometric orbit is quite significant. Note that $`a_{derived}`$ is relatively large — 3.11 $`mas`$, which made the significant detection possible. However, this is not a planet-candidate case. The minimum mass suggests this secondary is a brown-dwarf candidate, whereas the astrometric orbit shows the secondary mass is in the stellar regime.

As of March 2001, the Encyclopedia of extrasolar planets included 49 planet candidates with minimum masses smaller than 13 $`M_{Jup}`$. We (Zucker & Mazeh 2001a) analyzed all but two of the planet candidates. One star had no Hipparcos data, and the other star is known to have two companions. Figure 2 presents our results by depicting $`a_{derived}`$ versus $`a_{99}`$. The figure indeed shows that all points fall to the right of the line $`a_{derived}=a_{99}`$. This means that all our derived astrometric motions are not significant in the 99% level. This includes the planets of $`\upsilon `$ And and HD 10697 whose derived orbits were previously published by us (Mazeh et al. 1999; Zucker & Mazeh 2000), but the new analysis renders their orbits less significant.

Note, however, that this does not mean that the orbits derived are all false. Figure 2 shows that some of the systems are close to the border line, indicating that the orbits of these systems were detected with significance close to 99%. The systems with significance higher than 90% are listed in Table 1. Here we list the Hipparcos number and the stellar name, the confidence level of the derived astrometric orbit, the derived semi-major axis, its uncertainty and the derived inclination; the derived secondary mass together with its 1$`\sigma `$ range.

To summarize, the combination of the Hipparcos data together with the radial-velocity measurements did not yield any astrometric orbit with significance higher than 99%. Apparently, the Hipparcos precision is not good enough to detect a 1 $`mas`$ orbit, even with the combination of the radial-velocity measurements. The analysis shows that the data are consistent with no astrometric detection at all, although one or two true astrometric orbits, which imply low inclinations, are still possible. However, such a finding would not prove that the orbits of the sample of planet candidates are not randomly oriented in space.

## 3 The Mass Distribution of the Extrasolar Planets

Assuming the orbits of the detected planet candidates are randomly oriented in space we can now proceed to derive their mass distribution. To do that we have to account for the unknown orbital inclination and for the fact that stars with too small radial-velocity amplitudes could not have been detected as radial-velocity variables. Therefore, planets with masses too small, orbital periods too large, or inclination angles too small were not detected.

Numerous studies accounted for the effect of the unknown inclination of spectroscopic binaries (e.g., Mazeh & Goldberg 1992; Heacox 1995; Goldberg 2000), assuming random orientation in space. Heacox (1995) calculated first the minimum-mass distribution and then used its relation to the actual mass distribution to derive the latter. This calculation amplified the noise in the observed data, and necessitated the use of quite heavy smoothing of the observed data. Mazeh & Goldberg (1992) introduced an iterative algorithm whose solution depended, in principle, on the initial guess.

Very recently Jorissen, Mayor, & Udry (2001a) studied the planet distribution by considering only the effect of the unknown inclination. Like Heacox (1995), Jorrisen, Mayor, & Udry derived first the distribution of the minimum masses and then applied two alternative algorithms to invert it to the distribution of planet masses. One algorithm was a formal solution of an Abel integral equation and the other was the Richardson-Lucy algorithm (e.g., Heacox 1995). The first algorithm necessitated some degree of data smoothing and the second one required a series of iterations. The results of the first algorithm depended on the degree of smoothing applied, and those of the second one on the number of iterations performed. In addition, Jorissen, Mayor, & Udry (2001) did not apply any correction to the observational selection effect.

We (Zucker & Mazeh 2001b) followed Tokovinin (1991, 1992) and constructed a maximum likelihood algorithm — MAXimum LIkelihood MAss, to derive an histogram of the mass distribution of the extrasolar planets. MAXLIMA derives the histogram directly by solving a set of numerically stable linear equations. It does not require any smoothing of the data, except for the bin size of the histogram, nor any iterative procedure. MAXLIMA also offers a natural way to correct for the undetected planets. This is done by considering each of the detected systems as representing more than one system with the same $`M_2\mathrm{sin}i`$, depending mainly on the period distribution. The details of the algorithm are given in Zucker & Mazeh (2001b).

To apply MAXLIMA to the current known sample of extrasolar planets we (Zucker & Mazeh 2001b) considered all known planets and brown dwarfs orbiting G- or K-star primaries as of April 2001. To acquire some degree of completeness to our sample we have decided to exclude planets with periods longer than 1500 days and with radial-velocity amplitudes smaller than 40 m/s. The values of these two parameters determine the correction of MAXLIMA for the selection effect, for which we assumed a period distribution which is flat in $`\mathrm{log}P`$. This choice of parameters also implies that our analysis applies only to planets with periods shorter than 1500 days. We further assumed that the primary mass is 1$`M_{}`$ for all systems.

The results of MAXLIMA are presented in the lower panel of Figure 3 on a logarithmic mass scale. The value of each bin is proportional to the estimated number of planets found in the corresponding range of masses in the known sample of planet candidates, after correcting for the undetected systems. To estimate the uncertainty of each bin we ran 5000 Monte Carlo simulations and found the r.m.s. of the derived values of each bin. Therefore, the errors plotted in the figure represent only the statistical noise of the sample. Obviously, any deviation from the assumptions of our model for the selection effect induces further errors into the histogram, the assumed period distribution in particular. This is specially true for the first bin, where the actual number of systems is small and the correction factor large.

To compare the mass distribution of the planet candidates with that of the stellar secondaries we plot (Zucker & Mazeh 2001b) the latter on the same scale in an adjacent panel of Figure 3. We plot here only two bins, with masses between 100 and 1000 $`M_{Jup}`$, using subsamples of binaries found by the Center for Astrophysics (=CfA) radial-velocity search for spectroscopic binaries (Latham 1985) in the Carney & Latham (1987) sample of the high-proper-motion stars (Latham et al. 2001; Goldberg et al. 2001).

Note that the upper panel does not have any estimate of the values of the bins with masses smaller than 100 $`M_{Jup}`$. This is so because the CfA search does not have the sensitivity to detect secondaries in that range. On the other hand, the lower panel does include information on the bins below 100 $`M_{Jup}`$. This panel presents the results of the high-precision radial-velocity searches, and these searches could easily detect stars with secondaries in the range of, say, 20–100 $`M_{Jup}`$. The lower panel shows that the frequency of secondaries in this range of masses is close to zero.

The relative scaling of the planets and the stellar companions is not well known (see Zucker & Mazeh 2001b for a detailed discussion). Nevertheless the comparison is illuminating. It suggests that we have here two distinct populations, separated by a ’gap’ of about one decade of masses, in the range between 10 and 100 $`M_{Jup}`$. We will assume that the two populations are the giant planets, at the low-mass side of Figure 3, and the stellar companions at the high-mass end of the figure. The present analysis is not able to tell whether the gap extends up to 60, 80 or 100 $`M_{Jup}`$.

The gap between the two populations was already noticed by many previous studies (Basri & Marcy 1997; Mayor, Queloz, & Udry 1998; Mayor, Udry, & Queloz 1998; Marcy & Butler 1998). Those papers binned the mass distribution linearly. Here we follow our previous work (Mazeh, Goldberg & Latham 1998; Mazeh 1999a,b; Mazeh & Zucker 2001) and use a logarithmic scale to study the mass distribution, because of the large range of masses, 0.5–1000 $`M_{Jup}`$, involved. The gap or the brown-dwarf desert is consistent also with the finding of Halbwachs et al. (2000), who used Hipparcos data and found that many of the known brown-dwarf candidates are actually stellar companions.

The distribution we derived in Figure 3 suggests that the planet mass distribution is almost flat in $`\mathrm{log}M`$ over five bins — from 0.3 to 10 $`M_{Jup}`$. Actually, the figure suggests a possible slight rise of the distribution toward smaller masses. At the high-mass end of the planet distribution the mass distribution dramatically drops off at 10 $`M_{Jup}`$, with a small high-end tail in the next bin. Although the results are still consistent with zero, we feel that the small value beyond 10 $`M_{Jup}`$ might be real. The dramatic drop at 10 $`M_{Jup}`$ and the small high-mass tail agree with the findings of Jorissen, Mayor, & Udry (2001), despite the differences in the algorithm used to derive the distribution, and the logarithmic scale we use for the distribution.

## 4 Eccentricity and Period Distribution of the Two Populations

Having established the difference between the mass distribution of the giant planets and that of the low-mass secondaries in spectroscopic binaries, we turn now to compare the period and eccentricity distributions of the two populations. For the latter we use (Mazeh & Zucker 2001) the results of a very large radial-velocity study of the Carney & Latham (1987) high-proper-motion sample, which yielded about 200 spectroscopic binaries (Latham et al. 2001; Goldberg et al. 2001). Goldberg (2000) separated statistically between the binaries of the Galactic halo and those coming from the disk. We consider in this section only the 59 single-lined spectroscopic binaries (=SB1s) of the Galactic disk. For the giant planet sample we use again the sample of 66 planet candidates listed in Schneider (2001) as of April 2001.

Figure 4A shows the cumulative period distribution of the two samples. The figure suggests similar general trend, except in the two ends of the distributions. We therefore plotted in Figure 4B the two distributions only in the range between 10 and 1650 days. The similarity is astounding, since the two distributions are identical. Both are consistent with a straight line, which implies a flat distribution in log P.

We speculate that at the short period range, below 10 days, some dynamical interaction changed the distribution of either one or both distributions. Such an interaction could also change the eccentricity distribution of the orbits. In order not to be distracted by this possible interaction when we consider the eccentricity distribution, we choose to consider only the eccentricities of the orbits with periods between 10 and 1650 days. The cumulative distributions are plotted in Figure 5A. We again see a similar trend in both distributions, except in both ends of the range . To illuminate the difference we plotted the density distribution in Figure 5B. We derived the distribution by convolving the actual data points with a Gaussian kernel with a width of 0.08. It is clear that both distributions peak at about 0.2–0.4. However, the distribution of the spectroscopic binaries drops sharply toward zero, whereas the planet distribution does not. The eccentricity distribution of the binaries displays a tentative ’shoulder’ at the large eccentricities, whereas that of the planets displays such a possible shoulder at the small eccentricities.

Any paradigm that assumes the two populations were formed differently has to explain why their eccentricity as well as period distributions are so much alike. Although we do not try to explain any of the two similarities, we suggest in the next section a toy model that can generate the two eccentricity distributions.

## 5 A Toy Model to Generate the Eccentricity Distributions of the Two Samples

Consider a sample of low-mass companions that orbit their parent stars in circular Keplerian orbits. For simplicity let us choose the units such that the orbital radii of all orbits are of length unity, and so are their orbital tangential velocities. Now let us introduce a Gaussian scatter to the velocities of the companions of the sample, with two independent components. One component is tangential and the other is radial. The tangential component changes the moduli of the velocities, while the radial one changes mainly their directions.

The new scatter determined the new velocity distribution. Denote the center of the distribution by $`v_0`$ and its r.m.s. by $`\sigma _v`$. Suppose that the velocity angles are distributed around $`90^{}`$, with r.m.s of $`\sigma _\theta `$. Note that the distribution has three parameters, $`v_0`$, $`\sigma _v`$ and $`\sigma _\theta `$.

We can now calculate the eccentricity distribution of the sample, and see if such a simple-minded toy model can mimic the observed distributions of the giant planets and the low-mass companions. Figure 6 compares the two. We found that we can approximate the giant planet distribution with $`v_0=0.82`$ and $`\sigma _v=0.13`$, while $`\sigma _\theta =0`$, whereas the low-mass stellar companions necessitated $`\sigma _\theta =25^{}`$, $`\sigma _v=0.05`$ and $`v_0=1.08`$.

The fact that we succeeded to mimic the two actual distributions is not surprising. As the old statistical saying goes: “You can fit an elephant with any model with two parameters, and you can make him dance with three”. However, the specific values of the parameters found are somewhat intriguing. Suppose that both populations started with Keplerian circular orbits, and two mechanisms introduced the scatter into the two populations. Suppose the nature of the mechanism that operated on the planet population was dissipative, like the dissipation generated by an interaction of a planet with a swarm of small particles in a disk. Such a mechanism could decrease the velocity without changing its direction. This would result with a null $`\sigma _\theta `$ and $`v_0`$ less than unity, the difference being of the same order of $`\sigma _v`$. On the other hand, the spectroscopic binaries could be subject to a more chaotic, eruptive disturbing mechanism, like the gravitational interaction with a few large bodies. In such a process one could expect a spread of the velocity directions and moduli, without significantly changing $`v_0`$. This simple-minded picture is consistent with our findings.

We should emphasize that the aforementioned discussion is not meant to explain how the eccentricities were formed, nor why the two distributions are similar with some definite small differences. The model might only serve as a starting point for any theoretical study to account for the observed distributions.

## 6 The Paucity of Short-Period Massive Planets

In Section 3 & 4 we have discussed the distributions of masses, periods and eccentricities of the extrasolar planets. In this section we move to examine one aspect of the inter dependence of these variables. To explore this possible dependence we performed a Principal Component Analysis (e.g., Kendall & Stuart 1958), which immediately pointed out to the significant correlation between the (minimum) masses and periods of the extrasolar planets. This is depicted in Figure 7, in which we plotted the period as a function of the (minimum) masses of the known planets, as of April 2001. We choose to plot the two axes with logarithmic scales, because the frequency of planets is flat in log M and log P, as has been shown in previous sections.

Most of the correlation between the periods and masses of the extrasolar planets could be accounted for by a selection effect, that prevents planets which are not massive enough from being discovered if their periods are too long. Such systems have radial-velocity amplitude, $`K`$, which is too small to be detected by the present planet-search projects. This is easily seen in the small-mass–long-period corner of the diagram, bounded by the $`K=25\mathrm{m}\mathrm{s}^1`$ line. There are only four planets above this line. However, a close examination of Figure 7 reveals an additional feature — a significant paucity of planets at the opposite, large-mass–short-period corner of the diagram. Only three planets appear at that corner, all marked by a circle. This is certainly not a selection effect, because planets at that part of the diagram have the largest radial-velocity amplitude, and therefore are the easiest to detect.

It is not clear yet what is the shape of the area in which we find low frequency of planets. That corner might have a rectangular shape bordered by $`\mathrm{log}(M_2\mathrm{sin}i)=0.2`$ and $`\mathrm{log}P=1.5`$, or could be of a wedge shape, bordered by the line that goes from $`(\mathrm{log}\left(M_2\mathrm{sin}i\right),\mathrm{log}P)=(0,0.5)`$ to $`(1,1.5)`$.

The three planets that we find in the small-mass–long-period corner are Gls 86, HD 195019 and $`\tau `$ Boo. Interestingly enough, all three systems are wide binaries. Els et al. (2001) discovered very recently that the star Gls 86 has a brown-dwarf companion at about 20 AU projected separation. Pourbaix & Arenou (2001) pointed out that HD 195019 is a known visual binary with a companion fainter by about 3 mag., observed at a separation of 3.5 arc-sec in 1988 (Mason et al. 2001). The angular separation of HD 195019 (=WDS 20283+1846) translates to 130 AU projected separation for a parallax of 27 $`mas`$ (ESA 1997). The third star, $`\tau `$ Boo, is also a known visual binary (WDS 13473+1727), with an M2 companion. Apparently, the period is about two thousand years (Hale 1994) and the orbit is very eccentric. The separation between the two stars has been measured in 1991 to be 3.4 arc-sec (Mason et al. 2001), which translates to about 50 AU projected distance for a parallax of 64 $`mas`$ (ESA 1997). Planets in binary systems might go through different orbital evolution, and therefore might be considered as special cases. Thus, the low frequency of planets with large masses and short periods seems to be even more real than is seen from the figure.

Statistical assessment of the significance of the low frequency found in this part of the parameter space is under way. Very simple-minded calculations that ignore both the observational selection effects and the binarity of the three stars indicate a significance at the 2–3$`\sigma `$ level. Taking into account the selection effect and the binarity of the three stars makes the significance of the low frequency even higher.

The paucity of large-mass planets with short periods and consequently small orbits might be another clue to the formation and orbital evolution of the extrasolar planets. There are now two different scenarios that account for the existence of giant planets in close-in orbits. One of them, accepted by most of the astronomical community, assumes the planets were formed out of a disc of gas and dust at a distance of 5 AU or larger, and have migrated through interaction with the disc to their present position (e.g., Lin, Bodenheimer, & Richardson 1996). The other one is that the planets were formed by some in situ disc instability (Boss 1997). In principle, our findings can be accounted for by both scenarios.

From the migration point of view, our findings might indicate that most large-mass planets halted their migration at orbital radius of the order of 0.2 AU. Obviously, the more massive the planet is, the more angular momentum and energy have to be removed from its orbital motion to enable the migration. Angular momentum and energy could be absorbed by the disc of gas and dust through generation of density waves (e.g., Goldreich & Tremaine 1980; Ward 1997) or by a planetesimal disc through gravitational interaction with the planet (e.g., Murray et al. 1998; Del Popolo, Gambera, & Ercan 2001). A too massive planet might move in until the local inner disc cannot absorb its angular momentum and energy. Such a consideration might account for a continuous dependence of the final orbital period on the planetary mass.

Interestingly enough, some studies suggested different migration scenarios for planets with small and large masses (Ward 1997). Massive planets open a gap in the disc, and subsequently go through slow, type II, migration, while small planets do not open a gap in the disc and therefore go through a relatively fast, type I, migration. The apparent paucity of short-period massive planets is consistent with such an evolutionary separation between large and small planets, if we can assume that the separation between the two types of migration occurs at a mass of about 1 $`M_{Jup}`$, and that type II migration could halt at about 0.2 AU (e.g., Lin et al. 2000).

According to the instability scenario, the mass of the formed planet depends on the available mass in the disc at the region of instability (e.g., Boss 2000). At small distances the available mass might be smaller, a fact that could result in low frequency of massive planets with short periods.

The fact that all three planets with relatively large masses and short periods are found in binary systems is intriguing. The interaction of the secondary with the protoplanetary disc could modify the structure and evolution of the disc, and therefore the formation and evolution of the planet. We obviously need more data to see whether this feature is statistically significant.

In all the aforementioned scenarios, the paucity of massive planets with short-period orbits is a natural consequence of the formation and evolutionary mechanism. However, detailed theoretical models have to be worked out so we can compare the theory with the observations. If confirmed by the discovery of more planets, the interesting input of the present analysis is the actual boundaries of the low-frequency part of the diagram. A borderline at about 1.5 $`M_{Jup}`$ and at about 30 days can help us quantitatively understand the formation and evolutionary process of extrasolar planets.

## 7 Summary

The logarithmic mass distribution derived here shows that the planet candidates are indeed a separate population, probably formed in a different way than the secondaries in spectroscopic binaries. Surprisingly the eccentricity and period distributions, with some restrictions, are very much the same.

Furthermore, the two period distributions follow strictly a straight line. This indicates flat density distributions on a logarithmic scale, inconsistent with the Duquennoy & Mayor (1991) log-normal distribution. Interestingly, flat logarithmic distribution is the only scale-free distribution, and could be argued to be the most simple distribution. Maybe the two populations were formed by two different mechanisms that still have this scale-free feature in common (Heacox 1999).

The eccentricity distribution of the sample of giant planets and that of stellar companions are similar (Stepinski & Black 2001c). In spite of the similarity, they are not identical, especially if compared to the remarkable similarity between the two period distributions. The eccentricity distributions can be attained by Keplerian orbits whose velocities are normally disturbed in the tangential and the radial directions.

We found a significant paucity of large planets with short orbital periods, and point out how, in principle, one can account for this paucity.

## Acknowledgments

We acknowledge support from the Israeli Science Foundation through grant no. 40/00. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and the Washington Double Star Catalog maintained at the U.S. Naval Observatory.

References

Basri, G., & Marcy, G. W. 1997, in AIP Conf. Proc 393, Star Formation, Near and Far, eds. S. Holt & L.G. Mundy (New York: AIP), 228

Black, D. C. 1995, ARA&A, 33, 359

Boss, A. P. 1997, Science, 276, 1836

Boss, A. P. 2000, ApJL, 536, L101

Brown, T. M., Charbonneau, D., Gilliland, R. L., Noyes, R. W., & Burrows, A. 2001, ApJ, 552, 699

Carney, B. W., & Latham, D. W. 1987, AJ, 92, 116

Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45

Del Popolo, A., Gambera, M., & Ercan, N. 2001, MNRAS, 325, 1402

Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485

Els, S. G., Sterzik, M. F., Marchis, F., Pantin, E., Endl, M., & Kürster, M. 2001, A&A, 370, L1

ESA 1997, The Hipparcos and Tycho Catalogues, ESA SP-1200

Gatewood, G., Han, I., & Black, D. 2001, ApJ, 548, L61

Goldberg, D. 2000, Ph.D. thesis, Tel Aviv University

Goldberg, D., Mazeh, T., Latham, D. W., Stefanik, R. P., Carney, B. W., & Laird, J. B. 2001, submitted to A&A

Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425

Good, P. 1994, Permutation Tests — A Practical Guide to Resampling Methods for Testing Hypotheses, (New York: Springer-Verlag)

Halbwachs, J.-L., Arenou, F., Mayor, M., Udry, S., & Queloz, D. 2000, A&A, 355, 581

Hale, A. 1994, AJ, 107, 306

Han, I., Black, D., & Gatewood, G. 2001, ApJ, 548, L57

Heacox, W. D. 1995, AJ, 109, 2670

Heacox, W. D. 1999, ApJ, 526, 928

Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000, ApJ, 529, L41

Jorissen, A., Mayor, M., & Udry, S. 2001, A&A, in press, astro-ph/0105301

Kendall, M. G., & Stuart, A. 1966, The Advanced Theory of Statistics, vol. 3, (London: Griffin)

Latham, D. W. 1985, in IAU Colloq. 88, Stellar Radial Velocities, eds. A. G. D. Philip & D. W. Latham (Schenectady, L. Davis Press) 21

Latham, D. W., Stefanik, R. P., Torres, G., Davis, R. J., Mazeh, T., Carney, B. W., Laird, J. B., & Morse, J. A. 2001, submitted to A&A

Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606

Lin, D. N. C., Papaloizou, J. C. B., Terquem, C., Bryden, G., & Ida, S. 2000, in Protostars and Planets IV eds. V. Mannings, A. P. Boss, S. S. Russell (Tucson: University of Arizona Press), 1111

Lissauer, J.J. 1993, ARA&A, 31, 129

Marcy, G. W., & Butler, R. P. 1998, ARA&A, 36, 57

Marcy, G. W., & Butler, R. P. 2000, PASP, 112, 137

Marcy, G. W., Cochran, W. D., & Mayor, M. 2000 in Protostars and Planets IV eds. V. Mannings, A. P. Boss, S. S. Russell (Tucson: University of Arizona Press), 1285

Mason, B. D., Wycoff, G. L., Hartkopg, W. I., Douglass, G. G., & Worley, C. E. 2001, Washington Double Star Catalog 2001.0, U.S. Naval Observatory, Washington

Mayor, M., Queloz, D., & Udry, S. 1998, in Brown Dwarfs and Extrasolar Planets, eds. R. Rebolo, E.L. Martin, & M.R. Zapatero-Osorio (San Francisco: ASP), 140

Mayor, T., Udry, S., & Queloz, D. 1998, in ASP Conf. Ser. 154, Tenth Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, eds. R. Donahue & J. Bookbinder (San Francisco: ASP), 77

Mazeh, T. 1999a, Physics Reports, 311, 317

Mazeh, T. 1999b,, in ASP Conf. Ser. 185, IAU Coll. 170, Precise Stellar Radial Velocities, eds. J. B. Hearnshaw & C. D. Scarfe, (San Francisco: ASP), 131

Mazeh, T. et al. 2000, ApJ, 532, L55

Mazeh, T., & Goldberg, D. 1992, ApJ, 394, 592

Mazeh, T., Goldberg, D., & Latham, D. W. 1998, ApJ, 501, L199

Mazeh, T., Mayor, M., & Latham D. W. 1996, ApJ, 478, 367

Mazeh, T., & Zucker, S. 2001, in IAU Symp. 200, Birth and Evolution of Binary Stars, eds. B. Reipurth and H. Zinnecker (San Francisco: ASP), 519

Mazeh, T., Zucker, S., Dalla Torre, A., & van Leeuwen, F. 1999, ApJ, 522, L149

Murray, N., Hansen, B., Holman, M., & Tremaine, S. 1998, Science, 279, 69

Pourbaix, D. 2001, A&A, 369, L22

Pourbaix, D., & Arenou, F. 2001, A&A, 372, 935

Schneider, J. 2001, in Extrasolar Planets Encyclopaedia http://www.obspm.fr/planets

Stepinski, T. F., & Black, D. C. 2001a, in IAU Symp. 200, Birth and Evolution of Binary Stars, ed. B. Reipurth & H. Zinnecker (San Francisco: ASP) 167

Stepinski, T. F., & Black, D. C. 2001b, A&A, 356, 903

Stepinski, T. F., & Black, D. C. 2001c, A&A, 371, 250

Tokovinin, A. A. 1991, Sov. Astron. Lett., 17, 345

Tokovinin, A. A. 1992, A&A, 256, 121

Ward, W. R. 1997, Icarus, 126, 261

Zucker, S., & Mazeh, T. 2000, ApJ, 531, L67

Zucker, S., & Mazeh, T. 2001a, ApJ, in press (astro-ph/0107124)

Zucker, S., & Mazeh, T. 2001b, ApJ, in press (astro-ph/0106042)


format: text

text:
# First Order Transition in the Spin Dynamics of Geometrically Frustrated Yb2Ti2O7.

\[

## Abstract

Using neutron diffraction, 170Yb Mössbauer and $`\mu `$SR spectroscopies, we have examined the pyrochlore Yb2Ti2O7 where the Yb3+ S = 1/2 ground state has planar anisotropy. Below $``$ 0.24 K, the temperature of the known specific heat $`\lambda `$-transition, there is no long range magnetic order. We show that the transition corresponds to a first order change in the fluctuation rate of the Yb3+ spins. Above the transition temperature, the rate, in the GHz range, follows a thermal excitation law, whereas below, the rate, in the MHz range, is temperature independent indicative of a quantum fluctuation regime.

\]

Geometrically derived magnetic frustration arises when the spatial arrangement of the spins is such that it prevents the simultaneous minimisation of all interaction energies . In the crystallographically ordered pyrochlore structure compounds R2Ti2O7, the rare earth ions (R) form a sub-lattice of corner sharing tetrahedra and a number of these compounds have been observed to exhibit frustration related behaviour . The low temperature magnetic behaviour associated with a particular rare earth depends on the properties of the crystal field ground state and on the origin (exchange/dipole), size and sign of the inter-ionic interactions. For example, the recently identified spin-ice configuration  has been linked with an Ising-like anisotropy and a net ferromagnetic interaction. Most of the published work on the pyrochlores has concerned rare earth ions with Ising-like characteristics  and there has also been some interest in the properties of weakly anisotropic Gd3+ . Less attention has been paid to the case, considered here, where the ion has planar anisotropy.

To date, in systems where geometrical frustration may be present, two different low temperature magnetic ground states have been considered. First, under the influence of the frustration the system does not experience a magnetic phase transition and remains in a collective paramagnetic state with the spin fluctuations persisting as T $``$ 0 . Second, a long range ordered state is reached through a phase transition which may be first order . Our results for Yb2Ti2O7, obtained using neutron diffraction, 170Yb Mössbauer spectroscopy and muon spin relaxation ($`\mu `$SR), evidence a novel scenario: there is a first order transition which does not correspond to a transition from a paramagnetic state to a (long or short range) magnetically ordered state. The transition chiefly concerns the time domain, and involves an abrupt slowing down of the dynamics of short range correlated spins; below the transition temperature, these spins continue to fluctuate at a temperature independent rate.

We have established the background magnetic characteristics for Yb2Ti2O7 in a separate study . The Yb3+ ion crystal field ground state is a very well isolated Kramers doublet with a planar anisotropy, $`g_{}/g_z`$ $``$ 2.5 where $`g_z`$ and $`g_{}`$ are respectively the spectroscopic factors along and perpendicular to a local  axis. The net inter-ionic interaction is ferromagnetic (the paramagnetic Curie-Weiss temperature is 0.75 K) and it is dominated by exchange (the dipole-dipole interaction is relatively small due to the modest value, 1.15 $`\mu _B`$, see below, of the Yb3+ moment). Specific heat measurements  (Fig. 1, top, inset) have evidenced a sharp peak ($`\lambda `$ transition) at $``$ 0.2 K with an associated magnetic entropy of $`0.18R\mathrm{ln}2`$, and a broad peak centered near 2 K. The total magnetic entropy was estimated to be $`0.97R\mathrm{ln}2`$  suggesting there may be a small amount of missing low temperature magnetic entropy.

Single phase, polycrystalline Yb2Ti2O7 was prepared by heating the constituent oxides up to 1400 C with four intermediate grindings. The neutron diffraction measurements were made down to 0.11 K on the G41 line (wavelength: 0.2427 nm) at the Laboratoire Léon Brillouin (Fig. 1) and to 0.065 K on the D1B line at the Institute Laue-Langevin. Fig. 1 shows that the crystal structure does not change with temperature and that below 0.2 K, there are no magnetic Bragg peaks, neither isolated nor superposed on the nuclear peaks. Had the moments of 1.15 $`\mu _B`$ (to be evidenced below) undergone long range ordering, the peaks would have easily been seen. We estimate that a magnetic Bragg peak would have remained visible if the correlation length had exceeded 2 to 3 nm, so this provides an upper bound for the correlation length. An additional low angle contribution appears in the difference spectra below $``$ 2 K (Fig. 1, bottom). Its intensity initially grows as the temperature decreases but there is no significant change on crossing the temperature of the specific heat $`\lambda `$-transition. If we interpret this additional contribution in terms of short range ferromagnetic correlations, then using a model of non-interacting spherical clusters each having a uniform magnetisation (solid line on Fig. 1, bottom)  we obtain a sphere diameter of $``$ 1.5(2) nm. It is possible that not all of the additional signal seen for 2$`\theta <15^{}`$ is of magnetic origin . We consider our estimate of the magnetic correlation length below the specific heat $`\lambda `$-transition to correspond to an upper limit. The absence of long range order contradicts previous suppositions .

Selected 170Yb Mössbauer absorption spectra are shown on the left panel of Fig. 2. At 0.036 K, a five line spectrum is observed indicating there is a “static” hyperfine field ($`H_{\mathrm{hf}}`$) which we find amounts to 115 T. In the present case, “static” means the fluctuation frequency of the field is less than the lowest measurable 170Yb value of $``$ 15 MHz. Knowing that for Yb3+ the hyperfine field is proportional to the 4f shell magnetic moment, we obtain that each of the Yb3+ carries a magnetic moment of 1.15 $`\mu _B`$. As there is no long range order, the hyperfine field is associated with the short range correlated Yb3+ moments.

In the absence of a significant quadrupole hyperfine interaction, we cannot obtain the local direction of the Yb3+ magnetic moment by directly measuring the angle it makes with the principal axis of the electric field gradient (a  direction). Instead, we make use of the property that for an anisotropic Kramers doublet, the size of the spontaneous magnetic moment is linked to the angle $`\theta `$ it makes with the local symmetry axis (a  axis) through the relation: $`M_{\mathrm{Yb}}=\frac{1}{2}g_{}\mu _B/(\mathrm{cos}\theta \sqrt{r^2+\mathrm{tan}^2\theta })`$ where $`r`$ is the anisotropy ratio $`g_{}/g_z`$ . Using the measured moment and the known $`g`$ values, we find $`\theta =44(5)^{}`$. Thus each moment does not lie perpendicular to its local  axis as would be expected if the orientation were governed only by the crystal field anisotropy. The uniform tilting of each moment away from its local  axis is probably related to the combined influence of the crystal field anisotropy and of the (presumably anisotropic) exchange interaction with the non-collinear neighbouring moments.

When the temperature is increased to 0.23 K, an additional single line sub-spectrum appears. It is linked with the fraction of the Yb3+ whose moments fluctuate “rapidly” so that the magnetic hyperfine splitting is “motionally narrowed”. The two sub-spectra (see Fig. 2, left panel at 0.24 K) are both present up to 0.26 K evidencing the coexistence of regions with “static” and “rapidly fluctuating” moments. Fig. 2, right panel, shows that as the temperature increases, there is a progressive decrease in the relative weight of the “static” hyperfine field sub-spectrum. These features clearly evidence a first order transition, taking place at slightly different temperatures in the different parts of the sample.

The single line (sub)spectrum progressively narrows as the temperature increases. As magnetic correlations are still present above 0.24 K, we attribute this change to the progressive increase in $`\nu _M`$, the fluctuation rate of $`\stackrel{}{H}_{\mathrm{hf}}`$. The relation between the dynamic line broadening, $`\mathrm{\Delta }\mathrm{\Gamma }_R`$, and $`\nu _\mathrm{M}`$, is written $`\mathrm{\Delta }\mathrm{\Gamma }_R`$ = $`(\mu _IH_{\mathrm{hf}})^2/\nu _\mathrm{M}`$ where $`\mu _I`$ is the 170Yb nuclear moment . As shown on Fig. 4, when the temperature is lowered from $``$ 1 K to just above 0.24 K, the rate decreases from $``$ 15 to $``$ 2 GHz. This decrease is linked to the slowing down of the spin fluctuations which accompanies the development of the short range correlations. Below 0.24 K, $`\nu _\mathrm{M}`$ drops to a value which is less than the lowest measurable 170Yb Mössbauer value.

The $`\mu `$SR study  was carried out at the ISIS facility (some additional measurements were also made at the PSI facility) over the range 300 to 0.04 K, mostly in a longitudinal field of 2 mT so as to decouple the small contribution of the nuclear spins to the measured depolarisation. Fig. 3 shows typical results for the time dependence of the asymmetry which is written $`aP_z^{\mathrm{exp}}(t)`$ = $`a_zP_z(t)`$ \+ $`a_{\mathrm{bg}}`$ where the first term on the right hand side originates from the sample and the second is a temperature independent constant background contribution linked with the muons stopping in the silver sample holder and in the cryostat ($`a_{\mathrm{bg}}`$ 0.065 for ISIS and close to 0 for PSI).

From 300 to 0.275 K, $`P_z(t)`$ is well represented by an exponential relaxation function (solid lines at 4.3 and 0.275 K on Fig. 3a): $`P_z(t)`$ = $`\mathrm{exp}(\lambda _zt)`$ where $`\lambda _z`$ is the muon spin-lattice relaxation rate. From 100 to $``$ 4 K, $`\lambda _z`$ is small ($``$ 0.1 MHz) and in keeping with the paramagnetic nature of the Yb3+, it shows little thermal dependence. As the temperature is lowered below $``$ 1 K, $`\lambda _z`$ increases progressively to reach 0.52(3) MHz at 0.275 K. The origin of this increase is the same as that for the dynamic broadening $`\mathrm{\Delta }\mathrm{\Gamma }_R`$ of the Mössbauer spectra, i.e. the slowing down of the electronic spin fluctuations.

As the temperature is lowered below 0.275 K, $`P_z(t)`$ first becomes moderately non-exponential and then towards 0.2 K it abruptly becomes strongly non-exponential (Fig. 3a). Below 0.2 K, $`P_z(t)`$ is independent of temperature and there are no spontaneous oscillations in accordance with the neutron result that there is no long range order. The spectrum of Fig. 3b which extends to shorter times, clearly confirms this result. The shape of $`P_z(t)`$ below 0.2 K is reminiscent of a dynamic Kubo-Toyabe (KT) decay  associated with a slowly fluctuating ensemble of moments, the decrease of $`P_z(t)`$ beyond $``$$`\mu `$s indicating that fluctuations are still present. The KT decay provides an approximate account of the shape of $`P_z(t)`$ but there is a noticeable misfit below 0.5 $`\mu `$s (dashed lines on Fig. 3). A similar misfit was also observed for Y2Mo2O7 near 0.03 $`\mu `$s at 2.5 K . A better fit is provided by the Gaussian Broadened Gaussian (GBG) model , where the single Gaussian distribution of the KT model is replaced by a collection of distributions (solid lines at 0.200 K on Fig. 3). From the fit of Fig. 3a, the mean value of the GBG distribution is $`\mathrm{\Delta }_{\mathrm{LT}}/\gamma _\mu `$ $``$ 5.7 mT, the ratio of the width of the collection of distributions to the mean value is $``$ 0.38 and the electronic fluctuation rate $`\nu _\mu `$ is $``$ 1 MHz. All are independent of temperature below 0.2 K.

To obtain an estimate of the effective value of the spin fluctuation rate, $`\nu _\mu `$, at temperatures above 0.24 K, we insert the experimental values for $`\lambda _z`$ (obtained in 2 mT) into the expression $`\lambda _z=2\mathrm{\Delta }_{\mathrm{HT}}^2`$/$`\nu _\mu `$. The value for $`\nu _\mu `$ then depends on the value of $`\mathrm{\Delta }_{\mathrm{HT}}`$, the Yb3+ \- muon spin coupling above 0.24 K. With the choice $`\mathrm{\Delta }_{\mathrm{HT}}/\gamma _\mu `$ = 31.8 mT, the $`\nu _\mu `$ obtained from the muon measurements scale remarkably well with the $`\nu _\mathrm{M}`$ obtained from the Mössbauer measurements (Fig. 4). Their common thermal dependence is well fitted (solid line on Fig. 4) by the thermal excitation law $`\nu =\nu _0\mathrm{exp}[E_{\mathrm{bh}}/(k_\mathrm{B}T)]`$, with $`\nu _0`$ = 17 GHz and mean barrier height Ebh = 0.5 K. The associated drop in $`\nu _\mu `$ at 0.24 K amounts to a factor of $`10^3`$. Choosing a different value for $`\mathrm{\Delta }_{\mathrm{HT}}`$/$`\gamma _\mu `$ would lead to different rates above 0.24 K and to a different size for the drop. For example, with $`\mathrm{\Delta }_{\mathrm{HT}}`$/$`\gamma _\mu `$ = 5.7 mT (the value we obtained for $`\mathrm{\Delta }_{\mathrm{LT}}`$/$`\gamma _\mu `$ below 0.24 K), the drop is $`10^2`$ (this is also the minimum size of the drop allowed by the Mössbauer results) and with $`\mathrm{\Delta }_{\mathrm{HT}}/\gamma _\mu `$ = 85.6 mT, obtained by scaling from Tb2Ti2O7 , it is $`10^4`$. Combining the Mössbauer and the muon results thus indicates that the spin fluctuation rate undergoes a first order change by a factor of $`10^2`$ to $`10^4`$. We mention that between 2 and 0.275 K, we find that $`\lambda _z`$ does not have the usual Lorentzian dependence on the applied longitudinal field. This anomalous behaviour will be discussed in a separate publication.

The low temperature magnetic properties of Yb2Ti2O7 are therefore unusual. In the short range correlated region from $``$ 2 K to $``$ 0.24 K, the spin fluctuation rates follow a thermal excitation law. This crystallographically ordered compound thus possesses barriers against spin reorientation. On crossing the temperature of the specific heat $`\lambda `$-peak , no long range order appears, but there is a first order drop, of two to four orders of magnitude, in the fluctuation rates of the correlated moments. As $`T0`$, the fluctuations persist at a temperature independent rate of $``$ 1 MHz, and they take place between directions which make an angle of $`44^{}`$ relative to the local anisotropy axis (a  direction). The fluctuations of the moments thus conserve this angle and hence involve spin flips or spin spirals around a  axis.

Above 0.24 K, we speculate that the spin fluctuations could link more random directions. The first order drop in the fluctuation rate would then be associated with a change in the nature of the fluctuations such that below 0.24 K, the system explores a subset of the states that are explored at higher temperatures. The described transition, first order without the appearance of any long range correlations and with a change of dynamics, evidences properties which parallel aspects of the general gas-liquid transition or of the specific water-ice transition. To our knowledge, such a transition has not previously been observed in a magnetic system.

When exotic low temperature properties have previously been observed in the rare earth pyrochlores, they have related to ions with Ising character (Tb3+, Dy3+, Ho3+). The present results show that further novel behaviour is evidenced when the rare earth has a planar anisotropy. Although Yb2Ti2O7 does not correspond to the spin-ice scenario (in addition to the planar anisotropy, there is, at most, only a small missing entropy), it shares with the spin-ice state the characteristic of short range correlated moments which continue to fluctuate between specific directions as T $``$ 0. Here however, the low temperature state is reached through a novel route: a first order dynamical transition which can be viewed as a change from a thermally excited regime to a low temperature quantum (or tunneling) regime.

We thank Isabelle Mirebeau and Collin Broholm for useful discussions and Nadine Genand-Riondet for technical assistance.


text:
# Evidence for Neutrinoless Double Beta Decay

## Abstract

The data of the HEIDELBERG-MOSCOW double beta decay experiment for the measuring period August 1990 - May 2000 (54.9813 kg y or 723.44 molyears), published recently, are analyzed using the potential of the Bayesian method for low counting rates. First evidence for neutrinoless double beta decay is observed giving first evidence for lepton number violation. The evidence for this decay mode is 97$`\%`$ (2.2$`\sigma `$) with the Bayesian method, and 99.8$`\%`$ c.l. (3.1$`\sigma `$) with the method recommended by the Particle Data Group. The half-life of the process is found with the Bayesian method to be $`\mathrm{T}_{1/2}^{0\nu }=(0.818.3)\times 10^{25}\mathrm{y}`$ (95% c.l.) with a best value of $`1.5\times 10^{25}\mathrm{y}`$. The deduced value of the effective neutrino mass is, with the nuclear matrix elements from , $`m`$ = (0.11 - 0.56) eV (95$`\%`$ c.l.), with a best value of 0.39 eV. Uncertainties in the nuclear matrix elements may widen the range given for the effective neutrino mass by at most a factor 2. Our observation which at the same time means evidence that the neutrino is a Majorana particle, will be of fundamental importance for neutrino physics. PACS. 14.69.Pq Neutrino mass and mixing - 23.40.Bw Weak-interaction and lepton (including neutrino) aspects - 23.40.-s Beta decay; double beta decay; electron and muon capture.

The neutrino oscillation interpretation of the atmospheric and solar neutrino data, deliver a strong indication for a non-vanishing neutrino mass. While such kind of experiments yields information on the difference of squared neutrino mass eigenvalues and on mixing angles, the absolute scale of the neutrino mass is still unknown. Information from double beta decay experiments is indispensable to solve these questions . Another important problem is that of the fundamental character of the neutrino, whether it is a Dirac or a Majorana particle . Neutrinoless double beta decay could answer also this question. Perhaps the main question, which can be investigated by double beta decay with high sensitivity, is that of lepton number conservation or non-conservation.

Double beta decay, the rarest known nuclear decay process, can occur in different modes:

| $`2\nu \beta \beta \mathrm{decay}:`$ | $`A(Z,N)A(Z+2,N2)+2e^{}+2\overline{\nu }_e`$ (1) |
| --- | --- |
| $`0\nu \beta \beta \mathrm{decay}:`$ | $`A(Z,N)A(Z+2,N2)+2e^{}`$ (2) |
| $`0\nu (2)\chi \beta \beta \mathrm{decay}:`$ | $`A(Z,N)A(Z+2,N2)+2e^{}+(2)\chi `$ (3) |

While the two-neutrino mode (1) is allowed by the Standard Model of particle physics, the neutrinoless mode ($`0\nu \beta \beta `$ ) (2) requires violation of lepton number ($`\mathrm{\Delta }`$L=2). This mode is possible only, if the neutrino is a Majorana particle, i.e. the neutrino is its own antiparticle (E. Majorana , G. Racah , for subsequent works we refer to , for some reviews see ). First calculations of $`0\nu \beta \beta `$ decay based on the Majorana theory have been done by W.H. Furry .

Neutrinoless double beta decay can not only probe a Majorana neutrino mass, but various new physics scenarios beyond the Standard Model, such as R-parity violating supersymmetric models, R-parity conserving SUSY models, leptoquarks, violation of Lorentz-invariance, and compositeness (for a review see ). Any theory containing lepton number violating interactions can in principle lead to this process allowing to obtain information on the specific underlying theory. The experimental signature of the neutrinoless mode is a peak at the Q-value of the decay.

The HEIDELBERG-MOSCOW double beta decay experiment in the Gran Sasso Underground Laboratory searches for double beta decay of $`{}_{}{}^{76}Ge^{76}Se`$ \+ 2 $`e^{}`$ \+ (2$`\nu `$) since 1990. It is the most sensitive double beta experiment since almost eight years now. The experiment operates five enriched (to 86$`\%`$) high-purity $`{}_{}{}^{76}Ge`$ detectors, with a total mass of 11.5 kg, the active mass of 10.96 kg being equivalent to a source strength of 125.5 mol $`{}_{}{}^{76}Ge`$ nuclei. This is the largest source strength in use.

The high energy resolution of the Ge detectors assures that there is no background for a $`0\nu \beta \beta `$ line from the two-neutrino double beta decay in this experiment.

The unique feature of neutrinoless double beta decay is that a measured half-life allows to deduce information on the effective Majorana neutrino mass $`m`$, which is a superposition of neutrino mass eigenstates. The half-life is given, when ignoring contributions from right-handed weak currents, by .

$$[T_{1/2}^{0\nu }(0_i^+0_f^+)]^1=C_{mm}\frac{m^2}{m_e^2}$$
(1)
$$m=|m_{ee}^{(1)}|+e^{i\varphi _2}|m_{ee}^{(2)}|+e^{i\varphi _3}|m_{ee}^{(3)}|,$$
(2)
Here $`|m_{ee}^{(i)}|\mathrm{exp}(i\varphi _i)`$ ($`i=1,2,3`$) are the contributions to $`m`$ from individual mass eigenstates, with $`\varphi _i`$ denoting relative Majorana phases connected with CP violation and, $`C_{mm}`$ denotes a nuclear matrix element, which can be calculated.

The effective mass is closely related to the parameters of neutrino oscillation experiments. The importance of $`m`$ for solving the problem of the structure of the neutrino mixing matrix and in particular to fix the absolute scale of the neutrino mass spectrum which cannot be fixed by $`\nu `$ \- oscillation experiments alone, has been discussed in detail in e.g. .

In this paper, we present a new, refined analysis of the data obtained in the HEIDELBERG-MOSCOW experiment during the period August 1990 - May 2000 . The analysis concentrates on the neutrinoless decay mode which is the one relevant for particle physics (see, e.g. ). A description of the HEIDELBERG-MOSCOW experiment has been given recently in , and in .

Fig. 1 shows the combined spectrum of the five enriched detectors obtained over the period August 1990 - May 2000, with a statistical significance of 54.981 kg y (723.44 molyears), around the Q-value of double beta decay. The latter has been determined to be Qββ = 2039.006(50) keV in a recent precision experiment . Fig. 2 shows the spectrum obtained with detectors Nr. 1,2,3,5 over the period August 1990 - May 2000, with a significance of 46.502 kg y. Fig. 3 shows the spectrum of single site events (SSE) obtained for detectors 2,3,5 in the period November 1995 - May 2000, under the restriction that the signal simultanously fulfills the criteria of all three methods of pulse shape analysis we have recently developed  and with which we operate all detectors except detector 1 (significance 28.053 kg y). Double beta events are single site events confined to a few mm region in the detector, corresponding to the track length of the emitted electrons, in contrast to e.g. multiple scattering $`\gamma `$ \- events. From simulation we expect that about 5$`\%`$ of the double beta single site events should be seen also as MSE. This is caused by bremsstrahlung of the emitted electrons .

All the spectra are obtained after rejecting coincidence events between different Ge detectors and events coincident with activation of the muon shield. The spectra which are taken in bins of 0.36 keV are shown in Figs. 1,2,3 in 1 keV bins, which explains the broken number in the ordinate. We do the analysis of the measured spectra with (Fig. 1) and without the data of detector 4 (Figs. 2,3) since the latter does not have a muon shield and has the weakest energy resolution. We ignore for each detector the first 200 days of operation, corresponding to about three half-lives of $`{}_{}{}^{56}Co`$ ($`T_{1/2}`$ = 77.27 days), to allow for some decay of short-lived radioactive impurities.

The background rate in the energy range 2000 - 2080 keV is found to be (0.17 $`\pm `$ 0.01) events/ kg y keV (without pulse shape analysis) considering all data as background. This is the lowest value obtained in such type of experiments. The energy resolution at the Qββ value in the sum spectra is extrapolated from the strong lines in the spectrum to be (4.00 $`\pm `$ 0.39) keV for the spectra with detector 4, and (3.74 $`\pm `$ 0.42) keV (FWHM) for those without detector 4. The energy calibration of the experiment has an uncertainty of 0.20 keV. For more experimental details see .

We analyse the spectra shown in Figs. 1,2,3 with the following methods:

1. Bayesian method, which is used widely at present (see, e.g. ). This method is particularly suited for low counting rates, where the data follow a Poisson distribution, that cannot be approximated by a Gaussian (see, also ).

2. Method recommended by the Particle Data Group .

The $`\chi ^2`$ \- method is not used since it is expected to have some limitation for application at low counting rates, as also the Maximum Likelihood Method. The limitation is due to the Gaussian approximations inherent in these methods.

To be conservative, we first concentrate on the spectra without pulse shape analysis. For the evaluation, we consider the raw data of the detectors. After that we shall turn to the analysis of the spectrum obtained with pulse shape discrimination.

Figs. 4,5,6 show the result of a Bayesian peak detection procedure. One knows that the lines in the spectrum are Gaussians with a standard deviation $`\sigma `$ = 1.70 (1.59) keV corresponding to 4.0 (3.7) keV FWHM. Given the position of a line, we used Bayes theorem to infer the contents of the line and the level of a constant background. We first describe the procedure summarily and then give some mathematical details.

Bayesian inference yields the joint probability distribution for both parameters. The background level was integrated out. This yielded the distribution of the line contents. If the distribution has its maximum at zero contents, only an upper limit for the contents can be given and the procedure does not suggest the existence of a line. If the distribution has its maximum at non-zero contents, the existence of a line is suggested and one can define the probability K that there is a line with non-zero contents. It is associated with an error interval containing with probability K the true value of the contents. The Bayesian probability K agrees with the confidence of classical statistics in the case of a Gaussian likelihood function.

For every energy E of the spectra of Figs. 1,2,3, we have determined the probability K that there is a line at E. All the remainder of the spectrum was considered to be background in this search. The result is given on the left-hand sides of Figs. 4,5,6.

We define the Bayesian procedure in some more detail. It starts from the distribution $`p(x_1\mathrm{}x_M|\rho ,\eta )`$ of the count rates $`x_1\mathrm{}x_M`$ in the bins 1..M of the spectrum - given two parameters $`\rho ,\eta `$. The distribution $`p`$ is the product

$$p(x_1\mathrm{}x_M|\rho ,\eta )=\underset{k=1}{\overset{M}{}}\frac{\lambda _{k}^{}{}_{}{}^{x_k}}{x_k!}e^{\lambda _k}$$
(3)
of Poissonians for the individual bins. The expectation value $`\lambda _k`$ is the superposition

$$\lambda _k=\rho [\eta f_1(k)+(1\eta )f_2(k)]$$
(4)
of the form factors f1 of the line and f2 of the background; i.e. f1(k) is the Gaussian centered at E with the above-mentioned standard deviation and f2(k)$``$f2 is a constant.

Since

$$\underset{k=1}{\overset{M}{}}f_\nu (k)=1for\nu =1,2,$$
(5)
one has

$$\underset{k=1}{\overset{M}{}}\lambda _k=\rho .$$
(6)
Hence, $`\rho `$ parametrizes the total intensity in the spectrum, and $`\eta `$ is the relative intensity in the Gaussian line.

Prior distributions are given by Jeffreys’ rule ($`\mathrm{\S }`$5.35 of ). The parameter $`\rho `$ is integrated out and the posterior distribution P($`\eta |x_1\mathrm{}x_M`$) of the relative contents $`\eta `$ of the spectral line is obtained.

The peak detection procedure yields lines at the positions of known weak $`\gamma `$-lines from the decay of $`{}_{}{}^{214}Bi`$ at 2010.7, 2016.7, 2021.8 and 2052.9 keV . In addition, a line centered at 2039 keV shows up. This is compatible with the Q-value  of the double beta decay process. We emphasize, that at this energy no $`\gamma `$-line is expected according to the compilations in . We do not find indications for the lines from $`{}_{}{}^{56}Co`$ at 2034.7 keV and 2042 keV discussed earlier  (Figs. 4,5).

Bayesian peak detection suggests a line at Qββ whether or not one includes detector Nr. 4 without muon shield (Figs. 4,5). The line is also suggested in Fig. 6 after removal of multiple site events (MSE), see below.

On the left-hand side of Figs. 4,5,6, the background intensity (1-$`\eta `$) identified by the Bayesian procedure is too high because the procedure averages the background over all the spectrum (including lines) except for the line it is trying to single out. Therefore on the right-hand side of Figs. 4,5,6, the peak detection procedure is carried out within an energy interval that does not contain (according to the left-hand side) lines other than the one at Qββ. Still the interval is broad enough (about $`\pm `$ 5 standard deviations of the Gaussian line) for a meaningful analysis. We find the probability K = 96.5$`\%`$ that there is a line at 2039.0 keV in the spectrum shown in Fig. 1. This is a confidence level of 2.1 $`\sigma `$ in the usual language. For the spectrum shown in Fig. 2, we find a probability for a line at 2039.0 keV of 97.4$`\%`$ (2.2 $`\sigma `$). In this case the number of events is found to be 1.2 to 29.4 with 95$`\%`$ c.l.. It is 7.3 to 22.6 events with 68.3$`\%`$ c.l..The most probable number of events (best value) is 14.8 events.

We also applied the method recommended by the Particle Data Group . This method (which does not use the information that the line is Gaussian) finds a line at 2039 keV on a confidence level of 3.1 $`\sigma `$ (99.8$`\%`$ c.l.). Such value is obtained in a wide range of the assumed width of the signal and the energy range of evaluation.

From the analysis of the single site events, we find after 28.053 kg y of measurement 9 SSE events in the region 2034.1 - 2044.9 keV ($`\pm `$ 3$`\sigma `$ around Qββ) (Table 1). Analysis with the Bayesian method of the single site event spectrum (Fig. 3), as described before, shows again evidence for a line at the energy of Qββ (Fig. 6). Analyzing the range of 2032 - 2046 keV, we find the probability of 96.8$`\%`$ that there is a line at 2039.0 keV. We thus see a signal of single site events, as expected for neutrinoless double beta decay, precisely at the Qββ value obtained in the precision experiment of . The analysis of the line at 2039.0 keV before correction for the efficiency yields 4.6 events (best value) or (0.3 - 8.0) events within 95$`\%`$ c.l. ((2.1 - 6.8) events within 68.3$`\%`$ c.l.). Corrected for the efficiency to identify an SSE signal by successive application of all three PSA methods, which is 0.55 $`\pm `$ 0.10, we obtain a $`0\nu \beta \beta `$ signal with 92.4$`\%`$ c.l.. The signal is (3.6 - 12.5) events with 68.3 $`\%`$ c.l. (best value 8.3 events).

The PDG method gives a signal at 2039.0 keV of 2.8 $`\sigma `$ (99.4$`\%`$). The analysis given in Fig. 6 (left part) shows partly drastic differences in the relative intensities of other identified lines compared to the full spectra, which should reflect their different composition of single site and multiple site events.

We emphasize that we find in all analyses of our spectra a line at the value of Qββ. Under the assumption that the signal at Qββ is not produced by a background line of presently unknown origin, we can translate the observed number of events into half-lives. We give in Table 2 conservatively the values obtained with the Bayesian method and not those obtained with the PDG method. Also given in Table 2 are the effective neutrino masses $`m`$ deduced using the matrix elements of .

We derive from the data taken with 46.502 kg y the half-life $`\mathrm{T}_{1/2}^{0\nu }=(0.818.3)\times 10^{25}`$ $`\mathrm{y}`$ (95$`\%`$ c.l.). The analysis of the other data sets, shown in Table 2, and in particular of the single site events data, which play an important role in our conclusion, confirm this result.

The result obtained is consistent with the limits given earlier by the HEIDELBERG-MOSCOW experiment , considering that the background had been determined more conservatively there. It is also consistent with all other double beta experiments - which still reach less sensitivity. A second Ge-experiment, which has stopped operation in 1999 after reaching a significance of 9 kg y yields (if one believes their method of ’visual inspection’ in their data analysis) in a conservative analysis a limit of $`\mathrm{T}_{1/2}^{0\nu }>0.55\times 10^{25}\mathrm{y}`$ (90% c.l.). The $`{}_{}{}^{128}Te`$ geochemical experiment yields $`m_\nu <1.1`$ eV (68$`\%`$ c.l.), the $`{}_{}{}^{130}Te`$ cryogenic experiment yields $`m_\nu <1.8`$ eV and the CdWO4 experiment $`m_\nu <2.6`$ eV, all derived with the matrix elements of  to make the results comparable to the present value (for references see ).

Concluding we obtain, with more than 95$`\%`$ probability, first evidence for the neutrinoless double beta decay mode.

As a consequence, on the same confidence level, lepton number is not conserved. Further the neutrino is a Majorana particle. We conclude from the various analyses given above the effective mass $`m`$ to be $`m`$ = (0.11 - 0.56) eV (95$`\%`$ c.l.), with best value of 0.39 eV. Allowing conservatively for an uncertainty of the nuclear matrix elements of $`\pm `$ 50$`\%`$ (for detailed discussions of the status of nuclear matrix elements we refer to ) this range may widen to $`m`$ = (0.05 - 0.84) eV (95$`\%`$ c.l.).

In this conclusion, it is assumed that contributions to $`0\nu \beta \beta `$ decay from processes other than the exchange of a Majorana neutrino (see, e.g. ) are negligible.

With the limit deduced for the effective neutrino mass, the HEIDELBERG-MOSCOW experiment excludes several of the neutrino mass scenarios allowed from present neutrino oscillation experiments (see Fig. 7) \- allowing mainly only for degenerate and partially degenerate mass scenarios and an inverse hierarchy 3$`\nu `$ \- scenario (the latter being, however, strongly disfavored by a recent analysis of SN1987A. In particular hierarchical mass schemes are excluded at the present level of accuracy.

Assuming the degenerate scenarios to be realized in nature we fix - according to the formulae derived in  \- the common mass eigenvalue of the degenerate neutrinos to m = (0.05 - 3.4) eV. Part of the upper range is already excluded by tritium experiments, which give a limit of m $`<`$ 2.2 eV (95$`\%`$ c.l.) . The full range can only partly (down to $``$ 0.5 eV) be checked by future tritium decay experiments, but could be checked by some future $`\beta \beta `$ experiments (see, e.g. ). The deduced best value for the mass is consistent with expectations from experimental $`\mu e\gamma `$ branching limits in models assuming the generating mechanism for the neutrino mass to be also responsible for the recent indication for as anomalous magnetic moment of the muon . It lies in a range of interest also for Z-burst models recently discussed as explanation for super-high energy cosmic ray events beyond the GKZ-cutoff .

The neutrino mass deduced allows neutrinos to still play an important role as hot dark matter in the Universe.

New approaches and considerably enlarged experiments (as discussed, e.g. in ) will be required in future to fix the neutrino mass with higher accuracy.

Acknowledgments

One of the authors - H.V. Klapdor-Kleingrothaus - would like to thank all colleagues who have contributed to the experiment over the last decade of operation. They are particularly grateful to Dr. F. Petry, Dr. B. Maier, Dr. J. Hellmig and Dr. B. Majorovits who have developed the pulse shape discrimination methods and set up the VME electronics applied in the experiment since 1995. They thank Dr. T. Kihm for his highly efficient permanent support in the field of electronics and Dr. G. Heusser for his important contribution in the field of low-level techniques. They thank Mr. H. Strecker for his invaluable technical contributions. They thank Prof. D. Schwalm for his scientific interest and his efficient support of this project and would like to thank the director of the Gran Sasso Underground Laboratory, Prof. E. Bettini and also the former directors of Gran Sasso Profs. P. Monacelli and E. Bellotti, for invaluable support. Our thanks extend also to the technical staff of the Max-Planck Institut für Kernphysik and of the Gran Sasso Underground Laboratory. We thank Perkin Elmer (former ORTEC) Company, and in particular Dr. M. Martini, and Dr. R. Collatz for the fruitful cooperation. We are grateful to Dr. G. Sawitzki and Prof. W. Beiglböck from the Institute for Applied Mathematics of the University of Heidelberg for encouraging discussions. We acknowledge the invaluable support from BMBF and DFG of this project. We are grateful to the former State Committee of Atomic Energy of the USSR for providing the enriched material used in this experiment.