text:
# Statistical properties of exoplanets.

## 1 Introduction

The most remarkable feature of the sample of known extra-solar planets is undoubtedly the variety of their orbital characteristics, which challenges the conventional views of planetary formation. Amongst the most peculiar candidates are the giant planets orbiting very close to their parent stars, in contrast to the prediction of the standard model (e.g. Pollack et al. 1996) that they formed first from ice grains111Such grain growth provides the supposed requisite solid core around which gas could rapidly accrete, over the lifetime of the protoplanetary disk ($`10^7`$ y) in the outer region of the system where the temperature of the stellar nebula is not too high. The accommodation of this scenario to the present observations requires that the planets undergo a subsequent migration process bringing them close to the central star (see e.g. Lin et al. 1996; Ward 1997). Alternative points of view invoke in-situ formation (Wuchterl 2000a, b), possibly triggered through disk instabilities (Boss 2001, 2002). Note however that, even in such cases, subsequent disk-planet interactions leading to migration are expected to take place as soon as the planet is formed.

The number of known extra-solar planets exceeding 100, a statistically significant sample is now available from which we can determine meaningful distributions of planetary characteristics (Udry et al. 2003a; Marcy et al. 2003), and so try to point out useful constraints for the planet-building models and then possibly discriminate between the different proposed scenarios.

In a new series of papers, we will try to emphasize the emerging properties of planet-host stars and characteristics of the different orbital-element distributions of exoplanetary systems and discuss their implications for our understanding of planetary formation and evolution. This paper (Paper I) will be dedicated to the period/separation distribution of exoplanets. Sharp features in the period versus mass diagram start to emerge and provide strong observational constraints for the migration scenario. More specifically, we re-discuss the lack of massive planets on short-period orbits recently pointed out (Udry et al. 2002; Zucker & Mazeh 2002; Pätzold & Rauer 2002). We also emphasize the clear emergence of a planet shortage in the 10–100 d period range (a period valley) and we point out a sharp lower mass limit for the detected planets on “longer”-period orbits ($`P100`$ d). The later part of the paper will be dedicated to a discussion of the proposed explanations of these important observational findings and the possible implications for the future surveys.

A second paper in this series (Santos et al. 2003, Paper II) deals with the metallicity of the stars harbouring planets, providing accurate spectroscopic parameters for the most recent candidates, confirming known global properties and also pointing out new emerging features. In a third paper, Eggenberger et al. (2003, Paper III) examine the effect of binarity on planet formation and evolution. The presence of a close stellar companion is found to bring more massive planets close-in, in agreement with simulations by Kley (2001).

## 2 Observational emerging properties

### 2.1 No massive planets on short-period orbits

In a recent discussion of the statistical properties of massive planets versus lighter ones (Udry et al. 2002), we emphasized a shortage of massive planets on short-period orbits, based on the sample of about 80 exoplanets known at that time. This feature was also simultaneously pointed out by Zucker & Mazeh (2002). These authors furthermore verified its statistical significance and very interestingly examined the possible influence of binarity on the mass-period relation of exoplanets, an indication of potential different formation and evolution processes for planets in binaries and planets around single stars.

With now more than 100 very low-mass candidates (with $`m_2\mathrm{sin}i`$ 18 MJup), the lack of short-period massive planetary companions becomes even clearer, as seen in Fig. 1. When we neglect the multiple-star systems222As determined by dedicated adaptive optics programmes or unveiled by previous spectroscopic measurements (see Paper III; Udry et al. 2003b, c, for a review) (upper-right panel) in which the planetary formation or evolution could follow different paths (Zucker & Mazeh 2002; Udry et al. 2003b, c, Paper III), a complete void of candidates is observed in the diagram for masses larger than $``$ 2 MJup and periods smaller than $``$ 100 days333The 4 known candidates of planets orbiting an evolved star have also been discarded. They however have no influence on the results of this study. The only remaining point is HD 168443 b, member of a possible multi brown-dwarf system (Marcy et al. 2001; Udry et al. 2002).

Several processes have been proposed to explain the lack of massive planets on short-period orbits. In the context of the migration scenario, they mainly follow two different approaches: i) type II migration (after a gap opens in the disk) is shown to be less effective for massive planets i.e. massive planets stay farther out than lighter ones, or ii) when the planet reaches the central regions, some process related to planet-star interactions provokes mass transfer from the planet to the star - decreasing the mass of the former - or leads massive planets to fall into the central star. These approaches will be discussed in more detail in Sect. 3.

The possibility that multi-planet chaotic interactions send the lighter candidates in the inner regions (or out) of the system whereas the massive ones stay in the outer part may also be invoked (Rasio & Ford 1996; Weidenschilling & Marzari 1996). The frequency of planets ending very close to the central star seems however to be small (Ford et al. 2003).

### 2.2 A period valley between 10 and 100-d periods

In the lower panels of Fig. 1 and upper panel of Fig. 2, another very interesting feature is coming out of the period/separation distribution. We observe a shortage of planets with periods between roughly 10 and 100 days. This observational property is mainly due to the light candidates ($`m_2\mathrm{sin}i`$$``$ 2 MJup). As seen in the previous section, massive planets orbiting single stars are almost exclusively found on longer-period orbits. They probably form and stay further out (see discussion in Sect. 3). On the other hand, lighter planets are found at all distances from their star444except the lightest ones ($`m_2\mathrm{sin}i`$$``$$`0.75`$ MJup) at periods longer than $``$ 100 days (see next section). This mass-dependent behaviour is also clearly illustrated by the cumulative functions presented in the lower panel of Fig. 2 for two planetary mass regimes.

The features appearing in the planetary period distribution – peak at short periods and rise at intermediate periods – seem significant and not due to observational biases. The peak at short period is formed by the pile-up of migrating planets, stopped close to the central star. As seen in Fig. 1, it is almost exclusively composed of the lowest mass planets ($`m_2\mathrm{sin}i`$$``$ 0.75 MJup) and not due to observational biases as the detection limit at 10 ms-1 is at $`P`$$``$ 4000 days for $`m_2\mathrm{sin}i`$ = 0.75 MJup (with $`M_1`$ = 1 M and $`e`$ = 0. The effect of Jupiter on the Sun is 12 ms-1). The rise of the distribution at longer periods is now emerging thanks to the increase of the timebase of the radial-velocity surveys. More and more long-period planets are being detected. This rise is significant in the sense that, for a given mass range, similar planets on shorter-period orbits would have been easier to detect than the actual ones. The observed valley in the period distribution can then just be seen as a transition region between two categories of planets which suffered different migration behaviours.

### 2.3 A sharp mass transition in the migration process?

A third important feature emerging from the mass-period diagram is the apparent lack of very light planets ($`m_2\mathrm{sin}i`$$``$ 0.75 MJup) with periods larger than $``$ 100 days. This feature is already visible in the upper-right diagram of Fig. 1 (limited by the dotted line). It becomes obvious when the planetary minimum mass is displayed with a log scale (Fig. 3).

Of course this feature could be related to the observational bias inherent to the radial-velocity technique for planet search that makes the detection more difficult for distant and/or lighter planets. To obtain a quantitative idea of the effect of this bias, we have also plotted in Fig. 3 lines indicating the locations of the radial-velocity signals (semi-amplitude $`K`$ of 30, 10 and 3 ms-1), expected on a solar-mass star due to planets on circular orbits with given minimum masses and semi-major axes555The primary mass and eccentricity effects are not expected to be important as the radial-velocity semi-amplitude scales with $`(1e^2)^{1/2}M_1^{2/3}`$. The line for $`K=10`$ ms-1 gives the approximate $`3\sigma `$-limit of the present most precise surveys reaching precisions of $``$ 3 ms-1 (e.g. Vogt et al. 2000; Queloz et al. 2001). Some known candidates have been found very close to this limit, at low mass and small distance or at high mass and large distance666We can also note that these candidates “close” to the detection limit do not present special values of the eccentricity and primary or planetary masses. However, we clear see an “intermediate” region (hatched area), delimited by $`m_2\mathrm{sin}i`$ = 0.75 MJup, $`P`$ = 100 d and $`K`$ = 10 ms-1, where no planetary candidates have been found. A striking feature of this area is its very sharp limit in mass.

Is this empty region just due to the mentioned observational bias and the small-number statistics? We have in particular to worry about the observed effect coming from the shortage of candidates with periods in the 10–100 day range, as pointed out in the previous section. So, to check the statistical significance of the emptiness of the hatched area in the figure (i.e. to see whether it mainly comes from the small number of known planets), we performed Monte-Carlo simulations. The idea is to obtain $`10^6`$ realizations of the diagram, reproducing the actual observational biases of the known sample, and estimate from them the probability to have no point in the considered zone.

There are basically two ways of producing samples resembling to the observations: either i) take realistic planetary parameter distributions and model the observational biases or ii) draw the test samples from the actually observed distributions in which the acting biases are built in. As the “real” distributions are not known a priori, we have adopted the second approach:

i) We randomly select the orbital eccentricities, the primary masses, the planet minimum masses and separations from the actual observed distributions of $`e`$, $`M_1`$, $`m_2\mathrm{sin}i`$ and $`\mathrm{log}a`$ of the detected exoplanets. This ensures in particular that we are taking into account the effect of the orbital eccentricity and of the stellar mass on planet detection. It also produces distributions that reproduce the observed brown-dwarf desert and the period valley.

ii) Since the distributions of separations are different for light and massive planets (Fig. 1), we do the selection for two mass regimes independently (limit at 4 MJup).

iii) We only use the distributions for planets orbiting single dwarf stars as 1) planets in binaries may have different characteristics (Zucker & Mazeh 2002, Paper III) and 2) the post main-sequence stellar evolution may drastically change the fate of short-period planets.

iv) We reject the events giving a radial-velocity semi-amplitude smaller than $`K_{\mathrm{lim}}`$ = 10 ms-1. They would hardly be detected by the present surveys. A less restrictive limit a 15 ms-1 is also considered to check the sensitivity of the results on the detection limit.

v) Each realization consists of a selection of 86 pseudo-planets, corresponding to the actual observations (64 “light” and 22 “massive” planets orbiting single dwarf stars, as on February 2003). For each of them we count the number $`n`$ of points in the “test zone”. Note that the detection limit at $`K_{\mathrm{lim}}`$ is no more as well defined as in Fig. 3 but it now takes into account the orbital eccentricity and primary-mass values.

The normalized histogram of the distribution of the number of points $`n`$ in the test zone from our $`10^6`$ realizations is shown in Fig. 4. The probability of having no detection is only 0.00028. A bootstrap of the procedure with 100 times $`10^5`$ realizations gives a value of 0.00025 $`\pm `$ 0.00006 (Fig. 5), in agreement with the previous result. Relaxing the detection limit to 15 ms-1 yields a probability of 0.004. It shows the low sensitivity of the result on the chosen radial-velocity detection limit. We thus can conclude that the considered area is indeed empty of planets with a confidence level of about 99.97 %.

This conclusion is only valid if the observational biases are correctly taken into account. Biases related to planetary parameter distributions and detection technique are by construction included in the simulations. However, other aspects play a role. For example activity-induced radial-velocity jitter may screen planet detection. Such an effect can be neglected if we suppose that the $`m_2\mathrm{sin}i`$, $`e`$, $`M_1`$ and $`a`$ parameters are not correlated with activity. Another concern relates to the difficulty encountered when trying to actually derive the values for the orbital elements, what is different from just detecting radial-velocity variability. Observation timing and phase coverage are then also important parameters. These remarks emphasize the caution needed when considering the above described results. We should see them more like a trend emerging from the data rather than a real proof.

In the next section, we will propose a tentative explanation for this feature, based on recent results on runaway migration studied by Masset & Papaloizou (2003). We will also discuss in the last section of this paper the implications on the radial-velocity planet-search survey of the paucity of light planets with intermediate periods.

## 3 Checking the theoretical approaches

The different features pointed out above in the various presentations of the period/separation–mass diagram will allow us to check the processes proposed to explain the observed distribution of periods for extra-solar planet candidates.

### 3.1 Low migration rate for massive planets

A first approach to explain the lack of massive planets on short-period orbits invokes a low migration efficiency for higher mass planets. This view is supported by recent simulations of single-planet migration. When the mass of the planet becomes of the order of a significant fraction of the characteristic disc mass with which it interacts, the inertia of the planet becomes important and slows down the orbital evolution (e.g. Trilling et al. 1998; Nelson et al. 2000). It also takes more time to open a gap in the disk for massive planets - eventually longer than the lifetime of the disk - and initiate type II migration (Trilling et al. 2002).

Using a model assuming a simple impulse approximation for the type II migration, no ad-hoc stopping mechanism in the center and neglecting type I migration777The extremely short timescale of the type I migration represents a serious theoretical problem for planet formation. However, suggestions have been recently made that co-orbital corotation torque in thin viscous disks may decrease its dramatic effect (Artymowicz 2003) and even stop or revert it for a large set of conditions (Masset 2002). This has however to be further examined, Trilling et al. (2002) show that, for identical initial disk conditions, close-in surviving pseudo-planets have smaller masses and distant surviving pseudo-planets have larger masses. They also find that planets forming farther out migrate less rapidly and that higher-mass planets migrate on longer timescales. Thus, there should be more massive planets at intermediate and large semi-major axes than found close-in. On the contrary, the population of short-period planets should be dominated by smaller-mass planets, what is indeed observed:

– The upper panel in Fig. 2 presents the logarithmic distribution of exoplanet separations, superimposed to the Trilling et al. (2002) results. As already pointed out by these authors, if we take into account the observational bias penalizing long-period planets in the radial-velocity surveys, the theoretical and observational results agree fairly well. The observed peak at small separations is not reproduced by the simulations of Trilling et al. because they have not introduced in their study any ad-hoc process to stop the migration close to the central star.

– Furthermore, the lower panel in Fig. 2 giving the cumulative functions for 2 mass regimes (limit at 2 MJup) clearly shows that low-mass planets form the predominant population of short-period orbits whereas more massive planets are located further out.

– Moreover, the lower panels of Fig. 1 indicates that you have to go further and further out to find more and more massive planets. The hatched histogram representing planets more massive than 4 MJup rises up to 1000 days (left panel) – although such massive planets should be easier to detect at smaller periods – whereas the maximum of the rise of the distribution of lighter-mass planets is around 400 days. If we decrease down to 2 MJup the limit between the 2 mass regimes, the mentioned features almost overlap (right panel), indicating a trend for the higher masses to be found at larger distances from the central star. This higher mass trend with period is explicited in Fig. 6 which displays the mean mass (filled circle) or higher mass (average on the 3 highest values; open circles) of planets in period smoothing windows of width $`\mathrm{log}P[\mathrm{days}]`$ = 0.2. The effect becomes visible around $`P`$ = 40 d and turns to be very effective around $`P`$ = 100 d. Detection limits at 10 and 30 ms-1 are represented by the dotted lines. They show that observational biases related to the radial-velocity technique are not responsible for the observed feature. This effect is not surprising considering first that massive planets are preferentially formed in the outer regions where a larger amount of building material is available and second that higher-mass objects migrate on longer timescales. Furthermore, as the disk dissipates, the migration of more massive planets slows down earlier since their larger angular momentum produces greater resistance to disk-induced migration (Armitage et al. 2002).

### 3.2 Runaway migration and lack of light planets with longer periods

A recent study by Masset & Papaloizou (2003) brings some insights for our understanding of the observed lack of small-mass planets ($`m_2\mathrm{sin}i`$$``$ 0.75 MJup) on intermediate-period orbits ($`P`$$``$ 100 d). These authors studied the effect of co-orbital corotation torque on migrating protoplanets in the case of massive protoplanetary disks for which the planet Hill radius and the disk thickness have comparable orders of magnitude. In particular they show that, if the mass deficit created by the radial drift with the planet of the material trapped in the co-orbital region is larger than the planet mass, the migration rate undergoes a runaway which can rapidly vary the protoplanet semi-major axis by a large amount (50% over a few tens of orbits). The authors state that: “this can happen only if the planet mass is sufficient to create a gap in its surrounding region and if the surrounding disk mass is larger than the planet mass. This typically corresponds to planet masses in the sub-Saturnian to Jovian mass range embedded in massive protoplanetary disks”.

Moreover, the simulations by Masset & Papaloizou (2003) show that the limit in planetary mass for the appearance of the runaway regime is very steep (their Figs. 12-14). The threshold conditions for runaway may vary widely depending on the disk mass, viscosity or thickness. It however depends relatively weakly upon the planet mass on the high-mass side i.e. their is a critical mass relatively well-constrained ($`M_{\mathrm{crit}}`$$``$ 1 MJup) under which runaway is likely provided that the protoplanetary disk is not too light weight, and above which runaway is impossible (Masset, priv. comm.). As our observed limit at $`m_2\mathrm{sin}i`$$``$ 0.75 MJup is compatible with this critical mass, there could well be a link between the statistical observed property and runaway migration.

In this picture, planets which happen to be massive enough ($`M`$$`>`$$`M_{\mathrm{crit}}`$) before reaching the central regions would end up as the scattered population on the right of Fig. 3, whereas planets that never grow enough would remain very mobile and would eventually be flushed into the very central region.

Note that another limit exist for runaway migration on the planet low-mass side. Thus, very light planets (typically much lighter than Saturn) do not encounter such type of fast migration. That can be the case for the outer icy planets of the Solar System.

### 3.3 Observed constraints for migration: summary

In summary, putting everything together, we can try to draw a global picture for the planet migration behaviour in function of the planet mass, in single-star systems:
i) The more massive planets (typically, masses larger than $``$ 4 MJup) form preferentially in the outer regions – where there is a large-enough material reservoir – and do not migrate much. None is observed within $``$ 0.5 AU from the central star.
ii) Intermediate-mass objects migrate more easily whatever the distance where they form, the migration rate depending on the local conditions (planet mass, disk mass, viscosity, etc). They are observed at all distances.
iii) The lighter planets (masses from sub-Saturnian to Jovian) migrate easily, possibly undergoing runaway migration bringing them in the central tenth of an AU from the star if the disk is massive enough. Most of them are actually found in the very central regions. None is observed with $`P`$$``$$``$ 100 days although radial-velocity surveys are now precise enough to detect them.
iv) We know, however, from our own Solar System that this is not true for the planets that are far away or much lighter than the known exoplanets. Thus, we can speculate that for the very light planets the migration efficiency is decreasing with distance or mass. This again fits into the the runaway migration scenario, the simulations showing the existence of a transition limit in mass, on the low-mass side, for the runaway to occur. Unfortunately, the precision achieved by the present surveys does not allow us to determine the separation or mass limits at which the migration slows down significantly. Good hope to answer this question is however brought by the future surveys at higher precision (like HARPS at 1 ms-1, Pepe et al. 2002) that will increase our detection sensibility by a large factor (curve at $`K`$ = 3 ms-1 in Fig. 3).

This observational view of planet migration may be changed by the influence of a perturbing stellar companion (Paper III) that will make the massive planets migrate inwards faster or make closer-in intermediate-mass planets grow bigger (Kley 2001).

Finally, we have to stress that although theory suggests that it is possible to form massive planets closer to the stars than the ice-boundary limit (Papaloizou & Terquem 1999; Bodenheimer et al. 2000; Sasselov & Lecar 2000), within the migration scenario, there is no observational requirement for planets to form at separations smaller than a few AU.

### 3.4 Central star-planet interaction

From the central star – planet interaction point of view, several scenarios were proposed to explain the lack of massive planets on short-period orbits as alternatives to the above developed arguments in the context of the migration scenario.

Pätzold & Rauer (2002) proposed to explain the lack of massive close-in planets by tidal interactions between the planet and its central star. If the planet orbital period is smaller than the stellar rotation period, tidal friction will spin up the star whereas the semi-major axis of the planetary orbit is decreased, eventually reaching the Roche zone of the central star on a short timescale. These authors pointed out that massive planets spiral in much faster than others, explaining thus the observed trend.

Trilling et al. (1998) had already addressed this problem in also considering mass transfer between the planet and the star. Some of the pseudo-planets in their simulations reach the Roche lobe limit, start to lose mass and thus experience an outward torque from the Roche lobe overflow. They estimate the maximum mass of surviving planets to be around 2 MJup, corresponding to what is indeed observed.

Another approach could relate to the evaporation of the planet atmosphere when the planet comes close to the central star, decreasing thus its mass. Up to now, studies based on the $`T_{\mathrm{eff}}`$ of the heated close-in planets have shown that hydrogen evaporation was not efficient in the case of the known Hot Jupiters (e.g. Mayor & Udry 2000). However, new developments based on the estimate of the exospheric temperature of HD 209458 b, evaluated to be much higher than $`T_{\mathrm{eff}}`$ (Vidal-Madjar et al. 2003), call for a complete reconsideration of the question.

These approaches are very interesting to explain the period distribution of planets very close to the central stars ($`a`$$``$ 0.1 AU, with a pile up around 3 days), especially considering that the peak at very short periods is almost entirely formed by the lightest mass planets ($`m_2\mathrm{sin}i`$$``$ 0.75 MJup). However, they do not explain the extended range of periods ($`P`$$``$ 100 d) in which no planets more massive than $``$ 2 MJup are found. This large interval supports the non-migration of massive companions rather than the disappearance of close-in migrating massive planets.

## 4 Summary and concluding remarks

To summarise the main observational features pointed out in this paper and their main implications, we have that:

1) No massive planets ($`m_2\mathrm{sin}i`$$``$$``$ 2 MJup) are found on short-period orbits ($`P`$$`<`$$``$ 100 d) around single stars. This is not an observational bias as these candidates are the easiest ones to detect.

2) The maximum mass of detected planets per period interval increases with distance to the central star (Fig. 6). This is also a solid result as massive planets found at a given distance are easier to detect closer in.

3) The above points 1) and 2) suggest that the migration rate of planets decreases with increasing mass of the planetary companion. This result agrees with recent simulations of migrating planets in viscous disks (Trilling et al. 1998, 2002; Nelson et al. 2000). We thus expect a large number of massive planets to be on long-period orbits and so to be still undetected because of the limited duration of the present surveys. A large number a lower-mass planets probably also exist on long-period orbits; they are however more difficult to detect. These represent primary targets for future higher-precision surveys (like e.g. HARPS, Pepe et al. 2002) whereas the youngest among the formers are interesting targets for direct imaging of planetary-type objects. This result is also very important for on-going radial-velocity planet searches from which we now expect an increasing number of planet candidates as the survey durations will grow. The “older” the survey, the higher the expected detection rate. In this context, the first epoch measurement of the targets is an important parameter for planet detection.

4) We observe a shortage of planets with periods in the 10–100 d range. This valley in the distribution is located just in-between the peak of light planets that have migrated inwards and were stopped close to the central star and the rise of the distribution due to the increasing number of detected planets with longer periods.

5) Up to now, no planet candidates with very low masses ($`m_2\mathrm{sin}i`$$``$ 0.75 MJup) have been detected on orbits with periods longer than $``$ 100 days. This trend seems significant and not due to small-number statistics or to the particular shape of the period or separation distributions. The limit is sharp in mass, indicating a strong transition in the migration process for different mass regimes. These features are in agreement with predictions of runaway migration simulations of planets lighter than $``$ 1 MJup in massive protoplanetary disks (Masset & Papaloizou 2003).

Note however that this lack of light planets on intermediate periods may not be confirmed for much lighter planets or for planets in the considered mass range but formed far out from the central star (e.g. Neptune- and Uranus-type planets), that are still out of reach of the present observational facilities and detection techniques.

In conclusion, we have pointed out a clear dependence of migration on mass and distance from the central star, providing important constraints for models of planetary formation (migration). More quantitative constraints for the different mass regimes will require a better detection capability of the surveys, to diminish the effect of the observational bias in the distributions and thus increase the confidence we can put in the derived trends, especially for the lower-mass planets.

###### Acknowledgements.

This study has benefited from fruitful discussions with Frederic Masset and Willy Benz. We thank the Geneva University and the Swiss NSF (FNRS) for their continuous support of our planet-search projects. Support to N.S. from Fundação para a Ciência e Tecnologia (Portugal) in the form of scholarships is gratefully acknowledged.


format: text

text:
# Chapter 1 ANTARES Status Report

## Chapter 1 ANTARES Status Report

### Abstract

The ANTARES Collaboration is building a neutrino telescope 2400 m below the Mediterranean sea close to the Southern French coast. The site is already linked to the shore station by a 40 km-long electro-optical cable (EOC) which transmits power and data. A prototype line and an instrumentation line for monitoring environmental parameters have been successfully deployed and connected to the EOC via the junction box, using the IFREMER manned submarine. The Collaboration, after years of dedicated R&D and deployments of prototype lines, is now ready to deploy the detector starting in spring 2004.

### 1. The ANTARES neutrino telescope

The ANTARES (Astronomy with a Neutrino Telescope and Abyss environmental RESearch) project started in 1996 and involves physicists and engineers from France, Germany, Italy, Russia, Spain, The Netherlands and the United Kingdom. ANTARES aims to detect neutrinos with $`E_\nu 10`$ GeV in order to investigate $`\nu `$ astrophysics, dark matter in the form of weakly interacting massive particles (WIMPs), monopoles and $`\nu `$ oscillations. Cherenkov light produced by relativistic charged particles is detected by a 3D-array of optical modules (OMs), pressure resistant glass spheres containing phototubes (PMTs). Photon arrival times and PMT charge amplitudes allow track and energy reconstruction.

Current predictions and upper limits from previous generation and currently running telescopes indicate that the expected signal from cosmic neutrino sources require very large detectors to be observed. About 200 $`\nu `$-induced upward-going muons/km2/yr are expected for a diffuse flux from an isotropic distribution of optically-thin extra-galactic sources equal to the Waxman & Bahcall limit of $`4.510^8E^2`$ GeV-1 cm-2 s-1 sr-1 . Galactic source luminosities $`10^{35}`$ erg/s, achievable in the presence of compact accelerators and intense magnetic fields, are required to produce a rate larger than 10 events/yr from 100 TeV neutrinos in a km2 array. In 1 year of data-taking the ANTARES expected sensitivity will surpass that of current arrays and that expected for 600 live-days of AMANDA II . The success of this experiment, with an effective area $`>0.02`$ km2 for $`E_\nu >10`$ TeV and well reconstructed events, will be a milestone demonstrating the feasibility of an underwater $`km^3`$ detector in the Mediterranean, complementing a similar array in the South Polar ice. Two $`km^3`$-size detectors, one in each hemisphere, are needed to cover the whole sky (including the Galactic Centre which is not accessible from the South Pole using upward-going neutrinos), and guarantee a cross-check of systematic errors arising from different Cherenkov media properties.

The ANTARES site (42 50’N, 6 10’E) is well-shielded from the atmospheric muon background by 2400 m of sea water and has been selected after an intense program of sea campaigns dedicated to measurements of water transmission properties. The absorption length is about 60 m at 470 nm and mainly determines the size of the instrumented region and the PMT spacing. The effective scattering length ***Effective means that the scattered photon distribution, forward peaked in sea water, is accounted for. is more than 200 m, this is considerably larger and less depth dependent than that of ice. The optical background, which is absent in ice, is due to $`\beta `$ decays of $`{}_{}{}^{40}K`$ and to a continuous bioluminescence rate which combine to give a rate of about 60 kHz on a 10” PMT. Occasional, short bioluminescence bursts momentarily increase the rate up to MHz. These bursts induce a dead-time of 5% per PMT, however the detector dead-time from this source is far less due to the requirement for coincidences. The average light transmission loss of an OM due to bio-fouling and sedimentation is $`<2\%`$ at its equator after 1 year from deployment, this saturates with time .

Twelve lines will be deployed each with 75 OMs mounted in 25 triplets (storeys). Each OM has a 10” Hamamatsu PMT  oriented at $`45^{}`$ from the vertical. Storey separation is 14.5 m giving a total implemented height of about 350 m which starts 100 m above the sea bed. Lines are kept taut by buoys and are at an average distance of $`65`$ m from one another in an octagonal configuration. The production of the 900 OMs started in spring 2002. The readout of each PMT signal is shared by 2 Analog Ring Sampler ASICs which provide the analog signal time and charge digitization. The ARS implements a waveform (WF) shape-sensitive discrimination to distinguish single photoelectron-like pulse shapes (more than 98% of events) from larger pulses. Tests have shown that an overall time resolution of $`1`$ ns can be achieved mainly limited by the transit time spread (TTS) of the PMTs . The ARS’s together with a compass and tilt meters or hydrophones for the line positioning are located inside a titanium container which is common to each storey. Considering a singles rate of 70 kHz from $`{}_{}{}^{40}K`$ and bioluminescence rate of 2% of WF events, the typical data rate to shore will be $`7`$ MB/s/PMT; a $`100`$ PC farm on shore will process data. The DAQ system design is given in . Time calibrations are critical for the ANTARES experiment: a system of LED and laser beacons and a clock calibration system  will allow a relative time precision between OMs of $`0.5`$ ns and the absolute time will be determined with an accuracy of $`1`$ ms. The present planning indicates that the full detector will be installed by 2005.

From Nov. 1999 to June 2000 the Collaboration achieved a significant milestone in the deployment and operation of a “demonstrator line”, instrumented with 7 PMTs at a depth of 1100 m and connected to shore with a 37 km-long EOC. The shape of the atmospheric muon zenith angular distribution was reproduced, despite the small number of PMTs and their 1D-layout. The ANTARES relative and absolute positioning acoustic system of rangemeters, compasses and tilt miters was tested. Relative distances and absolute positioning were measured with an accuracy of $`5`$ cm and of $`1`$ m, respectively.

In Oct. 2001, the 40 km-long EOC for power and data transmission between the detector and the shore station in La Seyne sur Mer was deployed. Since Dec. 9, 2002, the heart of the forthcoming array, the junction box (JB), has been in communication with the shore station. It was deployed during a sea operation requiring the dredging and lifting of 2.5 km of the undersea EOC. On March 17, 2003 the first data were received from a prototype detection line equipped with 15 OMs (5 storeys corresponding to 1/5 of an ANTARES line). The line was deployed in Dec. 2002 and connected in Mar. 2003 to the underwater JB using the Nautile manned submarine of the French IFREMER oceanographic research agency. During the same mission, the Nautile connected a prototype instrumentation line (deployed on Feb. 12, 2003) incorporating instrumentation to monitor underwater environmental parameters (a pulsed laser calibration system, a deep-sea Doppler current meter, detectors to measure sound velocity, salinity and water transparency). The success of the submarine connections has proved the viability of the final detector configuration with 12 line inter-connections radiating from the JB. Data are currently being acquired from PMTs, tilt meters and compasses of the prototype line and from the instrumentation line. They are consistent with the expected single counting rate of around 60 kHz due to $`{}_{}{}^{40}K`$ $`\beta `$ decays and peaks in excess of 250 kHz from bioluminescence bursts. The status and results of the prototype lines are described in . Fig. 1 shows the layout of the future 12 line detector and underwater photographs of the prototype lines.

### 2. Expected physics performances

The ANTARES angular resolution is about 0.2 for $`E_\nu >10`$ TeV, where pointing accuracy is not limited by the $`\nu \mu `$ kinematics, but by the PMT TTS and by light scattering in water. The ANTARES sensitivity (90% c.l.) for point-like source searches to the upgoing $`\mu `$ flux induced by a typical $`E^2`$ differential $`\nu `$ flux is in the range between $`4÷5010^{16}`$ cm-2 s-1 after 1 yr. Search strategies and discovery potential are discussed in . The energy resolution and methods to reconstruct muon energy and parent neutrino spectra are discussed in . The sensitivity (90% c.l.) to $`E^2`$ diffuse differential neutrino fluxes achieved by rejecting the background with an energy cut of $`E_\mu 50`$ GeV is $`810^8`$ GeV cm-2 s-1 sr-1. The sensitivity of ANTARES in the search of WIMPs is given in . The upper limit on $`\nu `$-induced muon fluxes (90% c.l.) from neutralino annihilation in the Sun is at the level of 400 km-2 yr-1 for $`m_\chi 200`$ GeV.

### 3. References

1. Amram P. et al. 2003, Astrop. Phys. 19, 253

2. Barwick S.W. et al., astro-ph/0211269

3. Bouwhuis M. et al., A data acquisition system for the ANTARES neutrino telescope, this conference

4. Circella M. et al., Toward the ANTARES neutrino telescope: results from a prototype line, this conference

5. Heijboer A. et al., Point source searches with the ANTARES neutrino telescope, this conference

6. Hernández-Rey J. et al., Time Calibration of the ANTARES Neutrino Telescope, this conference

7. Romeyer et al., Muon energy reconstruction in ANTARES and its application to the diffuse neutrino flux, this conference

8. Thompson L. et al., Dark matter searches with the ANTARES neutrino telescope, this conference

9. Waxman E., Bahcall J.N. 1999, Phys. Rev. D59 (1999) 023002

10. Zornoza J. et al., Study of photomultiplier tubes for the ANTARES neutrino telescope, this conference


text:
# The inverse of the cumulative standard normal probability function.

## 1 Introduction

It would be difficult to overestimate the importance of the standard normal (or Gauss) distribution. It finds widespread application in almost every scientific discipline, e.g., probability theory, the theory of errors, heat conduction, biology, economics, physics, neural networks , etc. It plays a fundamental role in the financial mathematics, being part of the Black-Scholes formula , and its inverse is used in computing the implied volatility of an option . Yet, little is known about the properties of the inverse function, e.g., series expansions, asymptotic behavior, integral representations. The major work done has been in computing fast and accurate algorithms for numerical calculations .

Over the years a few articles have appeared with analytical studies of the closely related error function

$$\mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}\underset{0}{\overset{x}{}}e^{t^2}𝑑t$$

and its complement

$$\mathrm{erfc}(x)=\frac{2}{\sqrt{\pi }}\underset{x}{\overset{\mathrm{}}{}}e^{t^2}𝑑t\text{ .}$$

Philip  introduced the notation “inverfc($`x`$)” to denote the inverse of the complementary error function. He gave the first terms in the power series for inverfc($`x`$), asymptotic formulas for small $`x`$ in terms of continued logarithms, and some expressions for the derivatives and integrals. Carlitz , studied the arithmetic properties of the coefficients in the power series of inverfc($`x`$). Strecok  computed the first 200 terms in the series of inverfc($`x`$), and some expansions in series of Chebyshev polynomials. Finally, Fettis  studied inverfc($`x`$) for small $`x`$, using an iterative sequence of logarithms.

The purpose of this paper is to present some new results on the derivatives, integrals, and asymptotics of the inverse of the cumulative standard normal probability function

$$N(x)=\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^xe^{\frac{t^2}{2}}𝑑t$$

which we call $`S(x)`$. In section 2 we derive an ODE satisfied by $`S(x)`$, and solve it using a power series. We introduce a family of polynomials $`P_n`$ related to the calculation of higher derivatives of $`S(x)`$. In section 3 we study some properties of the $`P_n`$, such as relations between coefficients, recurrences, and generating functions. We also derive a general formula for $`P_n`$ using the idea of “nested derivatives”, and we compare the $`P_n`$ with the Hermite polynomials $`H_n`$.

In section 4 we extend the definition of the $`P_n`$ to $`n<0`$ and use them to calculate the integrals of $`S(x)`$. We also compute the integrals of powers of $`S(x)`$ on the interval $`[0,1]`$. Section 5 is dedicated to asymptotics of $`S(x)`$ for $`x0,x1`$ using the function Lambert W. With the help of those formulas we derive an approximation to $`S(x)`$ valid in the interval $`[0,1]`$ with error $`\epsilon `$, $`\left|\epsilon \right|<0.0023`$. Finally, appendix A contains the first 20 non-zero coefficients in the series of $`S(x)`$, and the first 10 polynomials $`P_n`$.

## 2 Derivatives

###### Definition 1

Let $`S(x)`$ denote the inverse of

$$N(x)=\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^xe^{\frac{t^2}{2}}𝑑t$$

satisfying

$$SN(x)=NS(x)=x$$
(1)
In terms of the error function $`\mathrm{erf}(x)`$,

$$N(x)=\frac{1}{2}\left[\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)+1\right]$$

###### Proposition 2

$`S(x)`$ satisfies the IVP

$$S^{\prime \prime }=S(S^{})^2$$
(2)
$$S\left(\frac{1}{2}\right)=0,\text{ }S^{}\left(\frac{1}{2}\right)=\sqrt{2\pi }$$

Proof. Since $`N(0)=1/2`$, in follows that $`S(1/2)=0`$. From (2.1) we get

$$S^{}[N(x)]=\frac{1}{N^{}(x)}=\sqrt{2\pi }e^{\frac{x^2}{2}}=\sqrt{2\pi }e^{\frac{S^2[N(x)]}{2}}$$

Substituting $`N(x)=y`$ we have

$$S^{}(y)=\sqrt{2\pi }e^{\frac{S^2(y)}{2}},S^{}(\frac{1}{2})=\sqrt{2\pi }$$

Differentiating $`\mathrm{ln}[S^{}(y)]`$, we get (2.2).

###### Proposition 3

$$S^{(n+2)}(x)=\underset{i=0}{\overset{n}{}}\underset{j=0}{\overset{i}{}}\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{i}{j}\right)S^{(ni)}(x)S^{(ij+1)}(x)S^{(j+1)}(x),n0.$$

Proof. Taking the n$`^{\text{th}}`$ derivative of (2.2), and using Leibnitz’s Theorem, we have

$`S^{(n+2)}`$ $`={\displaystyle \underset{i=0}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{i}}\right)S^{(ni)}(S^{}S^{})^{(i)}`$
$`={\displaystyle \underset{i=0}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{i}}\right)S^{(ni)}{\displaystyle \underset{j=0}{\overset{i}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)S^{(ij+1)}S^{(j+1)}`$
$`={\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{i}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)S^{(ni)}S^{(ij+1)}S^{(j+1)}.`$
###### Corollary 4

If $`D_n=\frac{d^nS}{dx^n}(\frac{1}{2})`$ , then

$$D_{2n}=0,\text{ }n0.$$

Putting $`D_n=(2\pi )^{\frac{n}{2}}C_n`$ , we can write

$$S(x)=\underset{n0}{}(2\pi )^{\frac{2n+1}{2}}\frac{C_{2n+1}}{(2n+1)!}(x\frac{1}{2})^{2n+1}$$

where

$$C_1=1,\text{ }C_3=1,\text{ }C_5=7,\text{ }C_7=127,\text{ }\mathrm{}.$$

###### Proposition 5

$$S^{(n)}=P_{n1}(S)(S^{})^nn1$$

where $`P_n(x)`$ is a polynomial of degree $`n`$ satisfying the recurrence

$$P_0(x)=1,P_n(x)=P_{n1}^{}(x)+nxP_{n1}(x),n1$$
(3)
so that

$$P_1(x)=x,P_2(x)=1+2x^2,P_3(x)=7x+6x^3,\mathrm{}.$$

Proof. We use induction on $`n`$. For $`n`$ = 2 the result follows from (2.2). If we assume the result is true for $`n`$ then

$`S^{(n+1)}`$ $`=[P_{n1}(S)(S^{})^n]^{}`$
$`=P_{n1}^{}(S)S^{}(S^{})^n+P_{n1}(S)n(S^{})^{n1}S^{\prime \prime }`$
$`=P_{n1}^{}(S)(S^{})^{n+1}+P_{n1}(S)n(S^{})^{n1}S(S^{})^2`$
$`=[P_{n1}^{}(S)+nSP_{n1}(S)](S^{})^{n+1}`$
$`=P_n(S)(S^{})^{n+1}`$
Since $`P_{n1}(x)`$ is a polynomial of degree $`n1`$ by hypothesis, is clear that

$$P_n(x)=P_{n1}^{}(x)+nxP_{n1}(x)$$

is a polynomial of degree $`n`$.

###### Corollary 6

$$C_n=P_{n1}(0)$$

## 3 The polynomials $`P_n(x)`$

###### Lemma 7

If we write

$$P_n(x)=\underset{k=0}{\overset{n}{}}Q_k^nx^k$$

we have

$`Q_0^n`$ $`=Q_1^{n1}`$
$`Q_k^n`$ $`=nQ_{k1}^{n1}+(k+1)Q_{k+1}^{n1}k=1,\mathrm{},n2`$ (4)
$`Q_k^n`$ $`=nQ_k^{n1}k=n1,n`$
Proof.

$`{\displaystyle \underset{k=0}{\overset{n}{}}}Q_k^nx^k`$ $`=P_n`$
$`={\displaystyle \frac{d}{dx}}P_{n1}+nxP_{n1}`$
$`={\displaystyle \underset{k=0}{\overset{n1}{}}}Q_k^{n1}kx^{k1}+{\displaystyle \underset{k=0}{\overset{n1}{}}}nQ_k^{n1}x^{k+1}`$
$`={\displaystyle \underset{k=0}{\overset{n2}{}}}Q_{k+1}^{n1}(k+1)x^k+{\displaystyle \underset{k=1}{\overset{n}{}}}nQ_{k1}^{n1}x^k`$
###### Corollary 8

In matrix form (3.1) reads $`A^{(n)}Q^{n1}=Q^n`$ where $`A^{(n)}`$ $``$ $`\mathrm{}^{(n+1)\times n}`$ is given by

$$A_{i,j}^{(n)}=\{\begin{array}{c}i,j=i+1,i=1,\mathrm{},n1\\ n,j=i1,i=2,\mathrm{},n+1\\ 0,ow\end{array}$$

In other words, $`A^{(n)}`$ is a rectangular matrix with zeros in the diagonal, the numbers $`1,2,3,\mathrm{},n1`$ above it, the number $`n`$ below it and zeros everywhere else

$$A^{(n)}=\left[\begin{array}{cccccccc}0& 1& 0& 0& 0& \mathrm{}& 0& 0\\ n& 0& 2& 0& 0& \mathrm{}& 0& 0\\ 0& n& 0& 3& 0& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& 0& \mathrm{}& 0& n1\\ 0& 0& 0& 0& 0& \mathrm{}& n& 0\\ 0& 0& 0& 0& 0& \mathrm{}& 0& n\end{array}\right]$$

and

$$Q^n=\left[\begin{array}{c}Q_0^n\\ Q_1^n\\ Q_2^n\\ \mathrm{}\\ Q_n^n\end{array}\right]$$

With the help of these matrices, we have an expression for the coefficients of $`P_n(x)`$

$$Q^n=\underset{k=1}{\overset{n}{}}A^{(nk+1)}=A^{(n)}A^{(n1)}\mathrm{}A^{(1)}$$

###### Proposition 9

The polynomials $`P_n(x)`$ satisfy the recurrence relation

$$P_{n+1}(x)=\underset{i=0}{\overset{n}{}}\underset{j=0}{\overset{i}{}}\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{i}{j}\right)P_{ni1}(x)P_{ij}(x)P_j(x)$$

Proof.

$`P_{n+1}(S)(S^{})^{n+2}`$
$`=`$ $`S^{(n+2)}={\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{i}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)S^{(ni)}S^{(ij+1)}S^{(j+1)}`$
$`={\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{i}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)P_{ni1}(S)(S^{})^{ni}P_{ij}(S)(S^{})^{ij+1}P_j(S)(S^{})^{j+1}`$
$`=(S^{})^{n+2}{\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \underset{j=0}{\overset{i}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{i}{j}}\right)P_{ni1}(S)P_{ij}(S)P_j(S).`$
###### Proposition 10

The exponential generating function of the polynomials $`P_n(x)`$ is

$$e^{\frac{1}{2}S^2[N(x)+tN^{}(x)]\frac{x^2}{2}}=F(x,t)=\underset{k0}{}P_k(x)\frac{t^k}{k!}$$

Proof. Since

$`F(x,t)`$ $`={\displaystyle \underset{k0}{}}P_k(x){\displaystyle \frac{t^k}{k!}}`$
$`=1+{\displaystyle \underset{k1}{}}{\displaystyle \frac{d}{dx}}P_{k1}(x){\displaystyle \frac{t^k}{k!}}+{\displaystyle \underset{k1}{}}kxP_{k1}(x){\displaystyle \frac{t^k}{k!}}`$
$`=1+{\displaystyle \underset{k0}{}}{\displaystyle \frac{d}{dx}}P_k(x){\displaystyle \frac{t^{k+1}}{(k+1)k!}}+{\displaystyle \underset{k0}{}}xP_k(x){\displaystyle \frac{t^{k+1}}{k!}}`$
$`=1+{\displaystyle _0^t}{\displaystyle \frac{}{x}}F(x,s)𝑑s+xtF(x,t)`$
it follows that $`F(x,t)`$ satisfies the differential-integral equation

$$1+(xt1)F(x,t)+\frac{}{x}_0^tF(x,s)𝑑s=0$$
(5)
Differentiating (3.2) with respect to $`t`$ we get

$$xF(x,t)+(xt1)\frac{}{t}F(x,t)+\frac{}{x}F(x,t)=0$$

whose general solution is of the form

$$F(x,t)=e^{\frac{x^2}{2}}G\left(te^{\frac{x^2}{2}}+\sqrt{2\pi }\left[N(x)\frac{1}{2}\right]\right)$$

for some function $`G(z)`$.

From (3.2) we know that $`F(x,0)=1`$, and hence

$$G\left(\sqrt{2\pi }\left[N(x)\frac{1}{2}\right]\right)=e^{\frac{x^2}{2}},$$

which implies that

$$G(z)=e^{\frac{1}{2}S^2\left(\frac{z}{\sqrt{2\pi }}+\frac{1}{2}\right)}.$$

Therefore,

$$F(x,t)=e^{\frac{1}{2}S^2[N(x)+tN^{}(x)]\frac{x^2}{2}}$$

###### Definition 11

We define the “nested derivative” $`𝔇^{(n)}`$ by

$`𝔇^{(0)}[f](x)`$ $`1`$
$`𝔇^{(n)}[f](x)`$ $`={\displaystyle \frac{d}{dx}}\left\{f(x)\times 𝔇^{(n1)}[f](x)\right\},n1`$
###### Example 12

1. $$𝔇^{(n)}\left[e^{ax}\right]=n!a^ne^{nax}$$
2. $$𝔇^{(n)}\left[x\right]=1$$
3. $$𝔇^{(n)}\left[x^2\right]=(n+1)!x^n$$

###### Proposition 13

$$P_n(x)=e^{\frac{n}{2}x^2}𝔇^{(n)}\left[e^{\frac{x^2}{2}}\right]$$

Proof. We use induction on $`n`$. For $`n=0`$ the result follows from the definition of $`𝔇^{(n)}`$. Assuming the result is true for $`n1`$

$`P_n(x)`$ $`=P_{n1}^{}(x)+nxP_{n1}(x)`$
$`={\displaystyle \frac{d}{dx}}[e^{\frac{(n1)}{2}x^2}𝔇^{(n1)}(e^{\frac{x^2}{2}})]+nxe^{\frac{(n1)}{2}x^2}𝔇^{(n1)}(e^{\frac{x^2}{2}})`$
$`=(n1)xe^{\frac{(n1)}{2}x^2}𝔇^{(n1)}(e^{\frac{x^2}{2}})+e^{\frac{(n1)}{2}x^2}{\displaystyle \frac{d}{dx}}[𝔇^{(n)}(e^{\frac{x^2}{2}})]+`$
$`+nxe^{\frac{(n1)}{2}x^2}𝔇^{(n1)}(e^{\frac{x^2}{2}})`$
$`=e^{\frac{(n1)}{2}x^2}\left[x𝔇^{(n1)}(e^{\frac{x^2}{2}})+{\displaystyle \frac{d}{dx}}𝔇^{(n1)}(e^{\frac{x^2}{2}})\right]`$
$`=e^{\frac{(n1)}{2}x^2}e^{\frac{1}{2}x^2}{\displaystyle \frac{d}{dx}}\left[e^{\frac{1}{2}x^2}𝔇^{(n1)}(e^{\frac{x^2}{2}})\right]`$
$`=e^{\frac{n}{2}x^2}𝔇^{(n)}\left[e^{\frac{x^2}{2}}\right]`$
###### Summary 14

We conclude this section by comparing the properties of $`P_n(x)`$ with the well known formulas for the Hermite polynomials $`H_n(x)`$ . Since the $`H_n`$ are deeply related with the function $`N(x)`$, we would expect to see some similarities between the $`H_n`$ and the $`P_n`$.

| $`P_n(x)`$ | $`H_n(x)`$ |
| --- | --- |
| $`P_n(x)=e^{\frac{n}{2}x^2}𝔇^{(n)}(e^{\frac{x^2}{2}})`$ | $`H_n(x)=(1)^ne^{x^2}\frac{d^n}{dx^n}(e^{x^2})`$ |
| $`\underset{k0}{}P_k(x)\frac{t^k}{k!}=e^{\frac{1}{2}S^2[N(x)+tN^{}(x)]\frac{x^2}{2}}`$ | $`\underset{k0}{}H_k(x)\frac{t^k}{k!}=e^{2xtt^2}`$ |
| $`P_n(x)=P_{n1}^{}(x)+nxP_{n1}(x)`$ | $`H_n(x)=H_{n1}^{}(x)+2xH_{n1}(x)`$ |
| $`S^{(n)}=P_{n1}(S)(S^{})^n`$ | $`N^{(n)}=\left(\frac{1}{\sqrt{2}}\right)^{n1}H_{n1}\left(\frac{x}{\sqrt{2}}\right)N^{}`$ |

## 4 Integrals of S(x)

###### Definition 15

$$S^{(n)}(x)=\underset{0}{\overset{x}{}}\underset{0}{\overset{x_1}{}}\mathrm{}\underset{0}{\overset{x_{n1}}{}}S(x_n)𝑑x_n𝑑x_{n1}\mathrm{}𝑑x_1,n1.$$

###### Lemma 16

$$P_{n1}(x)=e^{\frac{n}{2}x^2}\left[P_{n1}(0)+\underset{0}{\overset{x}{}}e^{\frac{n}{2}t^2}P_n(t)𝑑t\right]$$
(6)
Proof. It follows immediately from solving the ODE for $`P_{n1}`$ in terms of $`P_n`$.

###### Proposition 17

Using (4.1 ) to define $`P_n(x)`$ for $`n<0`$ yields

$`P_1(x)`$ $`=x`$
$`P_2(x)`$ $`=1`$ (7)
$`P_3(x)`$ $`=\sqrt{\pi }e^{x^2}N\left(\sqrt{2}x\right)`$
and the relation

$$S^{(n)}=P_{n1}(S)(S^{})^n$$

still holds.

Proof. For $`n=0`$ we have

$`P_1(x)`$ $`=P_1(0)+x`$
$`S`$ $`=S^{(0)}=P_1(S)`$
so $`P_1(0)=0`$. For $`n=1`$

$$P_2(x)=e^{\frac{x^2}{2}}\left[P_2(0)+1\right]1$$

We can calculate $`S^{(1)}`$ explicitly by

$`S^{(1)}(x)`$ $`={\displaystyle \underset{0}{\overset{x}{}}}S(t)𝑑t={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \underset{\mathrm{}}{\overset{S(x)}{}}}S[N(z)]e^{\frac{z^2}{2}}𝑑z`$
$`={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \underset{\mathrm{}}{\overset{S(x)}{}}}ze^{\frac{z^2}{2}}𝑑z`$
$`={\displaystyle \frac{1}{\sqrt{2\pi }}}e^{\frac{S(x)^2}{2}}=(S^{})^1`$
Hence, $`P_2(0)=1`$.

Finally, for $`n=2`$

$$P_3(x)=e^{x^2}\left[P_3(0)\sqrt{\pi }N(\sqrt{2}x)+\frac{\sqrt{\pi }}{2}\right]$$

A similar calculation as the one above, making a change of variables $`t=N(z)`$ in the integral of $`S^{(1)}(x)`$ yields

$$S^{(2)}(x)=\frac{1}{2\sqrt{\pi }}N[\sqrt{2}S(x)]$$

and we conclude that $`P_3(0)=\frac{\sqrt{\pi }}{2}`$.

###### Corollary 18

$$S\left[2\sqrt{\pi }S^{(2)}\right]=\sqrt{2}S$$

###### Proposition 19

$$\underset{0}{\overset{1}{}}S^n(x)𝑑x=\{\begin{array}{c}\underset{i=1}{\overset{k}{}}(2i+1),n=2k,k1\\ 0,n=2k+1,k0\end{array}$$

Proof.

$`{\displaystyle \underset{0}{\overset{1}{}}}S^n(x)𝑑x`$ $`={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}z^nN^{}(z)𝑑z`$
$`={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}z^ne^{\frac{1}{2}z^2}𝑑z`$
$`={\displaystyle \underset{i=1}{\overset{k}{}}}(2i+1),n=2k,k1`$
## 5 Asymptotics

###### Definition 20

We’ll denote by $`LW(x)`$ the function Lambert W ,

$$LW(x)e^{LW(x)}=x$$
(8)
This function has the series representation

$$LW(x)=\underset{n1}{}\frac{(n)^{n1}}{n!}x^n,$$

the derivative

$$\frac{d}{dx}LW=\frac{LW(x)}{x[1+LW(x)]}\text{if }x0,$$

and it has the asymptotic behavior

$$LW(x)\mathrm{ln}(x)\mathrm{ln}[\mathrm{ln}(x)]x\mathrm{}.$$

###### Proposition 21

$`S(x)`$ $`g_0(x)=\sqrt{LW\left({\displaystyle \frac{1}{2\pi x^2}}\right)},x0`$
$`S(x)`$ $`g_1(x)=\sqrt{LW\left({\displaystyle \frac{1}{2\pi (x1)^2}}\right)},x1`$
Both functions $`g_0(x)`$ and $`g_1(x)`$ satisfy the ODE

$$g^{\prime \prime }=g(g^{})^2\left[1+\frac{2}{g^2(1+g^2)}\right]g(g^{})^2,\text{for }\left|g\right|\mathrm{}$$

Proof.

$`N(x)`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}e^{\frac{x^2}{2}}{\displaystyle \frac{1}{x}},x\mathrm{}`$
$`t`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}e^{\frac{S(t)^2}{2}}{\displaystyle \frac{1}{S(t)}},t0`$
$`S(t)e^{\frac{S(t)^2}{2}}`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }t}},t0`$
$`S^2(t)e^{S^2(t)}`$ $`{\displaystyle \frac{1}{2\pi t^2}},t0`$
Using the definition of $`LW(x)`$ we have

$$S^2(t)LW\left(\frac{1}{2\pi t^2}\right),t0$$

or

$$S(t)\sqrt{LW\left(\frac{1}{2\pi t^2}\right)},t0$$

The case $`x1`$ is completely analogous.

###### Corollary 22

Combining the above expressions, we can get the approximation

$$S(x)g_2(x)=(2x1)\sqrt{LW\left(\frac{1}{2\pi x^2(x1)^2}\right)}$$
(9)
good through the interval $`(0,1)`$.

We can refine it even more by putting

$`S(x)`$ $`g_3(x)=Q(x)\sqrt{LW\left({\displaystyle \frac{1}{2\pi x^2(x1)^2}}\right)}`$ (10)
$`Q(x)`$ $`=1+(62\sqrt{2\pi })x+(12+6\sqrt{2\pi })x^2+(84\sqrt{2\pi })x^3`$
where $`Q(x)`$ has been chosen such that

$$Q(0)=1,Q(1)=0,Q(1/2)=0,Q^{}(1/2)=\sqrt{2\pi },Q^{\prime \prime }(1/2)=0$$

## 6 Appendix

The first 10 $`P_n(x)`$ are

$`P_0(x)`$ $`=1`$
$`P_1(x)`$ $`=x`$
$`P_2(x)`$ $`=1+2x^2`$
$`P_3(x)`$ $`=7x+6x^3`$
$`P_4(x)`$ $`=7+46x^2+24x^4`$
$`P_5(x)`$ $`=127x+326x^3+120x^5`$
$`P_6(x)`$ $`=127+1740x^2+2556x^4+720x^6`$
$`P_7(x)`$ $`=4369x+22404x^3+22212x^5+5040x^7`$
$`P_8(x)`$ $`=4369+102164x^2+290292x^4+212976x^6+40320x^8`$
$`P_9(x)`$ $`=243649x+2080644x^3+3890484x^5+2239344x^7+362880x^9`$
$`P_{10}(x)`$ $`=243649+8678422x^2+40258860x^4+54580248x^6+`$
$`+25659360x^8+3628800x^{10}`$
The first few odd $`C_n`$ are

| $`n`$ | $`C_n`$ |
| --- | --- |
| 1 | 1 |
| 3 | 1 |
| 5 | 7 |
| 7 | 127 |
| 9 | 4369 |
| 11 | 243649 |
| 13 | 20036983 |
| 15 | 2280356863 |
| 17 | 343141433761 |
| 19 | 65967241200001 |
| 21 | 15773461423793767 |
| 23 | 4591227123230945407 |
| 25 | 1598351733247609852849 |
| 27 | 655782249799531714375489 |
| 29 | 313160404864973852338669783 |
| 31 | 172201668512657346455126457343 |
| 33 | 108026349476762041127839800617281 |
| 35 | 76683701969726780307420968904733441 |
| 37 | 61154674195324330125295778531172438727 |
| 39 | 54441029530574028687402753586278549396607 |
| 41 | 53789884101606550209324949796685518122943569 |


text:
# Very Low Temperature Tunnelling Spectroscopy in the heavy fermion superconductor PrOs4Sb12

## Abstract

We present scanning tunnelling spectroscopy measurements on the heavy fermion superconductor PrOs4Sb12. Our results show that the superconducting gap opens over a large part of the Fermi surface. The deviations from isotropic BCS s-wave behavior are discussed in terms of a finite distribution of values of the superconducting gap.

The intriguing magnetic and superconducting properties of most heavy fermion metals, which include e.g. coexistence of superconductivity with ferromagnetism, non-Fermi liquid behaviors or multiple superconducting phases, represent a challenge to our current understanding of condensed matter physicsHewson ; Aoki01 ; Hasselbach89 . The route towards the discovery of new heavy fermions is to synthesize materials that maintain degrees of freedom (e.g. local magnetic or electric moments), which do not undergo a phase transition upon cooling and are coupled to the electron bath. A large entropy can then be preserved down to low temperatures, and transferred over to the electrons. Pr3+ ions offer this possibility if the Pr is in a cubic point symmetry, where the crystal electric field ground state can be a non-magnetic, non-Kramers $`\mathrm{\Gamma }_3`$ doublet, which can provide the entropy needed to create a heavy fermion ground state if a cooperative Jahn-Teller transition can be avoidedYatskar96 ; Cox87 . PrInAg2 meets these criteria, and it was shown to be the first possible example of a Pr-based heavy fermion with a very large linear specific heat coefficient of 6.5 J/mol K2Yatskar96 . Unfortunately PrInAg2 did not superconduct down to 50mK, making it hard to independently ascertain that the large linear specific heat term was indeed of electronic origin. Recently superconductivity has been found with Tc=1.85K in another Pr based heavy fermionBauer02 , PrOs4Sb12, where large electronic specific heat coefficients $`\gamma `$ between 300 and 500 mJ/mol K2 have been reportedBauer02 ; Vollmer03 ; MapleLa ; Measson ; Aoki02 . The height of the jump of the specific heat at the superconducting transition, which compares well with BCS theory, shows that the large $`\gamma `$ is indeed associated to the conduction electrons, and that, in addition, the heavy electronic bath superconducts Bauer02 ; Vollmer03 ; MapleLa ; Measson . Moreover, two superconducting transitions have been clearly resolved in the specific heat at 1.6 and 1.85 K in high quality single crystalline samples, giving strong indications for the presence of multiple superconducting phasesVollmer03 ; MapleLa ; Measson . The situation seems analogous to the only other known stochiometric superconductor that presents multiple superconducting phases, UPt3Hasselbach89 ; Joynt02 . In that case, most present theoretical and experimental scenarios associate the Cooper pairing mechanism that leads to multiple phase superconductivity to magnetic fluctuationsJoynt02 . Instead, in PrOs4Sb12, interactions with fluctuating quadrupolar (electric) moments seem the most likely mechanism that drives the system to an unconventional, multiple phase superconducting stateBauer02 ; Vollmer03 ; MapleLa ; Measson ; Aoki02 ; Kotegawa03 ; MacLaughlin02 ; Izawa02 . Note also that the isostructural compound LaOs4Sb12, which does not show heavy fermion behaviorSugawara02 , superconducts with a lower critical temperature (1K MapleLa ).

The determination of the most fundamental superconducting properties of PrOs4Sb12 is clearly needed to understand the formation of unconventional, multiple phase superconductivity. One of the first and most important points is to try to resolve the structure of the superconducting gap in the low temperature, low magnetic field phase, which occupies the largest part of the phase diagramVollmer03 . Indirect information can be obtained by thermodynamic, transport, magnetic or NMR measurements, but the experiments that have been done up to now lead to contradictory results, and are therefore not conclusive. Specific heat does not seem adequate to study the superconducting gap, because it shows a high Schottky peak at low temperaturesBauer02 ; Vollmer03 ; MapleLa ; Measson ; Aoki02 . NQR measurements show an exponential decrease of 1/T1T at low temperaturesKotegawa03  associated to a well developed gap. The London penetration depth, as measured with muon spin relaxation, also appears to decrease exponentially at low temperaturesMacLaughlin02 . On the other hand, the angular dependent thermal conductivity under magnetic fields appears to be strongly modulated due to significant changes in the superconducting gap over the Fermi surface. This result has been associated to the presence of point nodesIzawa02 .

Here we present direct measurements of the superconducting gap in PrOs4Sb12, done with high resolution tunnelling spectroscopy studies in the superconducting phase with a Scanning Tunneling Microscope (STM). We find a superconducting density of states with no low energy excitations and a well developed superconducting gap.

We use a home built STM unit and electronics installed in a partially home built dilution refrigerator. We have previously tested our experimental set-up by measuring the superconducting properties of Al, which is possibly the best known superconducting material with a critical temperature (1.2K) of the same order of magnitude as PrOs4Sb12 (1.85K). The sharpest superconducting features that we could resolve up to now are the quasiparticle peaks shown in Fig.1a, where we represent the tunnelling spectroscopy when both tip and sample are of Al (prepared following Suderow02a ). The width of these features can be taken as a measure of the resolution in energy of our spectrometer, and therefore of the effective lowest temperature of our measurements. We find a width at half maximum of 16$`\mu `$V, which corresponds to 184mK. On the other hand, the best fit to the tunnelling spectra taken between a normal tip of Au and a sample of Al (Fig.1b) using conventional isotropic BCS s-wave theory gives the same value for the temperature T=190mK (within the 5% uncertainty of the fit; we also use $`\mathrm{\Delta }=175\mu eV`$). Therefore, we can take 190mK as the lowest measuring temperature that can be achieved at present with our set-up. Our resolution is similar to the one achieved in many planar junction experimentsWolf , with the important difference that in STM measurements the tunnelling current is between 4 to 6 orders of magnitude smaller.

Single crystals of PrOs4Sb12 were grown out of a ternary melt that was rich in both Os and Sb. A starting composition of Pr2Os16Sb82 was heated to 1200 C and cooled over 100 hours to 725 C and then decanted, revealing small cubic crystals. We measured three samples in eight different cool downs by placing an Au tip on optically neat and flat faces of the single crystals. Within a given cool down, we changed the macroscopic position of the Au tip on the surface and measured many (more than 50 in total) different scanning windows (100x100nm2), using the positioning capabilities of our xy table, as in previous work Rubio01 . The measured work functions were always of several eV and the topography was reproducible upon changes of the tunnelling current. It consists of inclined planes and bumps, with typical corrugations of about 4nm (inset of Fig.2), indicating that the crystallographic direction of the surface is not well defined at the nanoscopic length scales relevant for this experiment.

Typical spectra (Fig.2) show a clean density of states with no low energy excitations, demonstrating that the superconducting gap is well developed over a large part of the Fermi surface. The best fit with isotropic BCS s-wave theory is shown by the line in Fig.2. Using T=0.19K we obtain a value for the superconducting gap $`\mathrm{\Delta }`$=270$`\mu `$eV, which gives $`2\mathrm{\Delta }/k_BT_c`$=3.4, very close to the BCS value 3.53. Note that differences between the experiment and the fit are clearly resolved thanks to the high resolution in voltage of our experiment, evidenced by the curves shown in Fig.1.

The tunnelling current depends on the overlap of the electron wavefunctions of tip and sample near the surface, being therefore a sum of the contributions from electrons coming from different sheets of the Fermi surface having k-vectors with different orientations. If the superconducting gap is anisotropic in a given sheet, or it has different values in different sheets of the Fermi surface, or both, the tunnelling spectra reflect this distribution of values of the superconducting gap by showing more broadened coherence peaks, which produce a deviation from isotropic BCS s-wave behavior (Fig.2 and e.g. Ref.Hess90 ). Note that the presence of defects, grain boundaries or other perturbations that one can figure out to occur near the surface, can lead either to a decrease of the observed anisotropy through the mixing of the electronic wavefunctions from different parts of the Fermi surface caused by strong inter and/or intraband electronic scatteringSung67 ; Golubov97 ; Martinez03 ; Martinez03a , or to pair breaking. The former was not observed, and the latter always leads to an increased residual density of states at the Fermi level, which we can indeed observe in some locations of the surface, as discussed further on. Therefore, we believe that the deviations between the experiment and the fit shown in Fig.2 are intrinsic to the superconducting density of states in this compound.

Note that the conductance begins to increase at about 120$`\mu `$V, and the highest point of the coherence peak is located at 325$`\mu `$V. As a matter of fact, the shape of the tunnelling spectra we find is similar to the form of the spectra taken in the much studied material NbSe2Hess90 , where first experiments clearly showed up the presence of a distribution of values of the superconducting gap over the Fermi surfaceHess90 . Subsequent work has identified this distribution as coming from different gap values in different sheets of the Fermi surface in that compound (multiband superconductivity)Taillefer03 ; Yokoya03 . Whereas more work is clearly needed to understand the origin of the gap distribution in PrOs4Sb12, it is noteworthy to remark that strong changes in the mass renormalization in different sheets has been found in the de Haas van Alphen experiments of Ref.Sugawara02 . These changes may also lead to the distribution of values of the superconducting gap measured in our experiment. This strengthens both the idea that the mass renormalization and superconductivity are of the same origin, i.e. that the quadrupolar fluctuations favor superconducting correlations, as well as the possible multiband character of superconductivity in this compound.

The spectra are smeared when we increase the temperature, as shown in Fig.3a, and become flat above the bulk critical temperature (1.85K). The maximum in the derivative of the conductance within the quasiparticle peaks, which gives a good estimate of the mean value of the superconducting gap Martinez03  is shown in Fig.3b and it follows well the temperature dependence of the superconducting gap within BCS s-wave theory (line in Fig.3b). The curve can be extrapolated to a critical temperature of 1.8K, which is the bulk Tc (1.85K) value within our experimental error (10%).

The superconducting properties in this compound clearly differ from the ones in magnetically mediated Ce or U heavy fermions. Although we could not find any published STM measurements in the tunnelling regime and in the superconducting phase in these materials, there is compelling evidence from many different techniques sensitive to the superconducting density of states that a large amount of low energy excitations due to a strongly anisotropic superconducting gap is found in most cases. For instance, many experiments show now that the multiphase superconductor UPt3 presents a line node along the basal plane and nodes along the c axis of its hexagonal crystalline structure, with a superconducting gap that decreases by more than an order of magnitude at the nodes Joynt02 ; SudUPt3 . In the case of PrOs4Sb12, there are at present no data pointing towards extended gapless regions on the Fermi surface, as the one caused by a line of nodes, so that the superconducting gap appears to be opened in a much larger part of the Fermi surface than in UPt3.

The observed behavior has been reproduced in different places, but it is not found over the whole surface. In general, regions with no residual density of states, as shown in Fig.4 a and b, have typical sizes of several times the coherence length ($`\xi _0`$=12nmBauer02 ) and are surrounded by regions where a finite density of states appears at the Fermi level, shown in Fig.4 c and d. We can also easily find regions with much less well defined superconducting features (not shown in the figure). When we find a finite density of states at the Fermi level, the superconducting features also disappear at a temperature smaller than critical temperature of the bulk, indicating that the physical origin for the residual density of states at the Fermi level is some kind of strong pair breaking effect appearing near the surface, easily detected with our technique.

In conclusion, we performed a direct measurement of the superconducting gap through high resolution local tunnelling spectroscopy measurements in the heavy fermion superconductor PrOs4Sb12. Typical spectra demonstrate that the superconducting gap is well developed over a large part of the Fermi surface. The presence of a finite distribution of values of the superconducting gap over the Fermi surface can be inferred from deviations between the experiment and isotropic BCS s-wave behavior.

We acknowledge discussions with J.P. Brison, M. Crespo, J. Flouquet, F. Guinea, K. Izawa Y. Kitaoka, A. Levanyuk and J.G. Rodrigo and support from the ESF programme VORTEX, from the MCyT (Spain; grant MAT-2001-1281-C02-0), and from the Comunidad Autónoma de Madrid (07N/0053/2002, Spain). The Laboratorio de Bajas Temperaturas is associated to the ICMM of the CSIC. Ames Laboratory is operated for the U. S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was supported by the Director for Energy Research, Office of Basic Energy Sciences.