A search for starlight reflected from υ And's innermost planet

Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 1 February 2008 (MN LATEX style file v1.4) A search for starlight reflected from υ And’s innermost planet arXiv:astro-ph/0110577v1 26 Oct 2001 Andrew Collier Cameron1 ⋆ , Keith Horne 1, Alan Penny2 and Christopher Leigh1 1 School of Physics and Astronomy, Univ. of St Andrews, St Andrews, Scotland KY16 9SS Appleton Laboratory, Chilton, Didcot OX11 0QX, United Kingdom 2 Rutherford Received; accepted 2001 ABSTRACT In data from three clear nights of a WHT/UES run in 2000 Oct/Nov, and using improved Doppler tomographic signal-analysis techniques, we have carried out a deep search for starlight reflected from the innermost of υ (upsilon) And’s three planets. We place upper limits on the planet’s radius Rp as functions of its projected orbital velocity Kp ≈ 139 sin i km s−1 for various assumptions about the wavelengthdependent √ geometric albedo spectrum p(λ) of its atmosphere. For a grey albedo p we find Rp p < 0.98 RJ with 0.1% false-alarm probability (4-σ). For a Sudarsky, Burrows & Pinto (1999) Class V model atmosphere, the mean albedo in our 380-676 nm bandpass is hpi ∼ 0.42, requiring Rp < 1.51 RJ , while an (isolated) Class IV model with hpi ∼ 0.19 requires Rp < 2.23 RJ . The star’s vrot sin i ∼ 10 km s−1 and estimated rotation period Prot ∼ 10d suggest a high orbital inclination i ∼ 70 − 80◦ . We also develop methods for assessing the false-alarm probabilities of faint candidate detections, and for extracting information about the albedo spectrum and other planetary parameters from faint reflected-light signals. Key words: Planets: extra-solar – Planets: atmospheres – Planets: radii 1 INTRODUCTION The discovery of “hot Jupiters” – giant exoplanets in 3 to 5-day orbits about their parent stars (Mayor & Queloz 1995; Marcy & Butler 1998) – has overturned conventional theories of planetary-system formation. Theorists working to explain the solar system had posited that Jupiters could form only outside 5 AU, where planetesimals retain ice mantles. Inward migration was considered possible, but deemed unlikely since the planets would then be swallowed up by the parent star. Hot Jupiters have now revived the idea of inward migration, but require mechanisms that halt some of the planets near 0.05 AU. Theoretical studies of the structure and evolution of exoplanets are yielding testable predictions. Gas-giant models yield planet radii Rp ∼ 1.0−1.7 RJ up for 1 < Mp /MJ up < 10 (Saumon et al. 1996; Guillot et al. 1996; Burrows et al. 1997). Hot Jupiters have fully-convective interiors, so the core entropy fixes the radius at a given mass (Burrows et al. 2000). Isolated Jupiters will therefore cool and contract to about Jupiter’s size in a few Gyr, but strongly irradiated giant planets find it harder to cool. Consequently the massradius-age relation for hot Jupiters is determined largely by ⋆ E-mail: [email protected] (1) the age at which they arrived in their present close orbits, and (2) their atmospheric albedos, which control the atmospheric temperature-pressure gradient and hence the rate of cooling (Burrows et al. 2000). The discovery of transits of the planet orbiting HD 209458 (Charbonneau et al. 2000) provided conclusive proof of the existence and gas-giant nature of the “hot Jupiters”. Ultra-precise HST/STIS photometry of HD 209458 b’s transit profile yields a radius about 1.35 RJ up (Brown et al. 2001), although the planet has a mass only 0.63 MJ up . The next crucial step will be to measure albedos directly. Spectral models(Seager & Sasselov 1998; Marley et al. 1999; Sudarsky, Burrows & Pinto 1999) include thermal emission in the infrared and scattered starlight in the optical. The albedo is very sensitive to the depth of cloud formation. At the expected temperatures (Teff ∼ 1400K) and gravities (g ≥ 36 m s−2 ) of massive planets such as τ Boo b, possible cloud condensates include upper cloud decks of MgSiO3 (enstatite) and Al2 O3 (corundum) at pressures ∼ 5 to 10 bar, and a deeper layer of Fe grains at pressures a factor 2 or so higher. If silicate cloud decks form this deep in the atmosphere or are absent, much of the incident starlight is absorbed in pressure-broadened absorption lines of neutral sodium and potassium, leaving only a weak Rayleigh-scattered reflection signature at blue wavelengths. 2 A.Collier Cameron, K. Horne, A. Penny and C. Leigh 0.6 0.5 0.4 p 0.3 0.2 0.1 0.4 0.45 0.5 0.55 0.6 Wavelength , microns 0.4 0.45 0.65 1 0.8 Nl 0.6 0.4 0.2 0.5 0.55 0.6 Wavelength , microns 0.65 Figure 1. Geometric albedo spectra for the Class V (solid line), isolated (dashed) and irradiated (dot-dashed) Class IV models of Sudarsky, Burrows & Pinto (1999), plotted over the wavelength range we observed. Together with a grey model having p = 0.42, the Class V and isolated Class IV models were used to probe the wavelength dependence of candidate reflected-light signals. The absorption features near 0.404 and 0.589 µm are due to potassium and sodium respectively. The geometric albedo spectrum of Jupiter (short dashes) is shown for comparison, including methane absorption bands at 0.54 and 0.62 µm. The lower panel shows the wavelength sensitivity of our detection method, expressed as the deficit of detected photons Nλ per unit wavelength due to known absorption lines in the stellar spectrum. Nλ is expressed in units of 1011 photons per µm, and gives the total photon deficit for a typical group of 4 exposures. Sudarsky, Burrows & Pinto (1999) predict this “Class IV” albedo spectrum for planets with high surface gravities. At the lower surface gravities g ≃ 10 m s−2 expected of lowermass planets like υ And b and HD 209458 b, “Class V” models predict silicate cloud decks forming higher in the atmosphere at pressures around 0.3 bar, giving a geometric albedo about 60% that of Jupiter throughout most of the optical spectrum (Fig.1). Recent attempts to detect starlight reflected from τ Boo b have produced deep upper limits on the planet’s albedo. Charbonneau et al. (1999) established an upper limit on the fully-illuminated planet/star flux ratio ǫ0 ≃ 10−4 . Collier Cameron et al. (1999) claimed a detection with an average ǫ0 = 7.5 × 10−5 from data secured in 1998 and 1999. The radius inferred for the planet was, however, implausibly large. Deeper observations made early in 2000 (Collier Cameron et al. 2001) did not confirm this candidate detection, but produced a more stringent upper limit, ǫ0 < 3.5 × 10−5 between 387.4 and 586.3 nm. The appar- ent tidal synchronisation of τ Boo with the planet’s orbit suggests a low inclination of order 40◦ . At this inclination, the combined 1998-2000 WHT observations yielded a 99.9% upper limit on the geometric albedo p ≤ 0.22, assuming a planet radius 1.2RJ up as predicted by Burrows et al. (2000). This appears to rule out the presence of a high cloud deck with the high reflectivity of the Class V models, which would give p ≃ 0.4 over most of the visible spectrum. Here we report results from a programme of observations designed to detect the spectroscopic signature of starlight reflected from the innermost of the three giant planets orbiting υ And (Butler et al. 1999). In Section 2 we use the measured system parameters to determine the a priori probabilities for the planet’s orbital velocity amplitude and the fraction of the star’s light that it intercepts. In Section 2.4 we discuss the expected albedo and phase function that the planet might display. In Sections 3 and 4 we describe the acquisition and extraction of the echelle spectra. The methods used to remove the direct-starlight signature from the data, and to extract the planetary signal by combining the profiles of the thousands of photospheric lines recorded in each echellogram, are presented in Section 5. A matched-filter method for measuring the strength of the reflected-light signal is described in Section 6. Finally in Section 6 we determine upper limits on the planet’s radius for three different models of the planet’s albedo spectrum, and discuss the plausibility of a candidate reflected-light signal that appears in the data at the most probable velocity amplitude and signal strength. 2 SYSTEM PARAMETERS υ And (HR 458, HD 9826) is a late-F main-sequence star with parameters as listed in Table 1. High-precision radialvelocity studies have shown it to have a system of 3 planets (Butler et al. 1999). For the purposes of this paper we are concerned only with the innermost of these planets, whose properties (as determined directly from radial-velocity studies or inferred using the estimated stellar parameters) are also summarised in Table 1. Starlight reflected from the orbiting planet’s atmosphere leaves a detectable signature in observed spectra in the form of faint copies of each of the stellar absorption lines, Doppler shifted by the planet’s orbital motion, and greatly diminished in brightness because the planet intercepts only a fraction (1/4)(Rp /a)2 ∼ 10−4 of the starlight and only part of that is reflected back into space. By detecting and characterizing the planetary reflected light signal, we in effect observe the planet/star flux ratio eps ≡ fp /f⋆ as a function of orbital phase φ and wavelength λ. The information we aim to obtain (in order of increasing difficulty) comprises: (i) Kp , the planet’s projected orbital velocity. From this we learn the orbital inclination i and hence the planet mass Mp , since Mp sin i is known from the star’s Doppler wobble. (ii) ǫ0 , the maximum strength of the reflected starlight. From this we constrain the planet’s radius since ǫ0 = √ Rp /a p, where p is the geometric albedo. (iii) p(λ), the albedo spectrum, which depends on the composition and structure of the planetary atmosphere. A search for starlight reflected from υ And’s innermost planet 3 Table 1. Physical parameters for υAnd and its innermost planet. Parameter Value (Uncertainty) Star: Spectral type Teff M⋆ (M⊙ ) R⋆ (R⊙ ) log g [Fe/H] π (mas) v sin i (km s−1 ) log Prot (days) Age (Gyr) F8V 6100 (80) 1.3 (0.1) 1.63 (0.06) 4.0 0.09 to 0.17 74.25 (0.70) 9.5 (0.8) 1.08 (0.08) 2.6 to 4.8 Planet: Orbital period Porb (days) Transit epoch T0 (JD) K⋆ (m s−1 ) a (AU) Mp sin i (MJ up ) 4.61707 (0.00003) 2451821.457 (0.012) 74.5 (2.3) 0.059 (0.002) 0.73 (0.04) References Gonzalez (1997); Fuhrmann, Pfeiffer & Bernkopf (1998); Gonzalez Gonzalez (1997); Fuhrmann, Pfeiffer & Bernkopf (1998); Gonzalez Fuhrmann, Pfeiffer & Bernkopf (1998); Ford, Rasio & Sills (1999) Fuhrmann, Pfeiffer & Bernkopf (1998); Ford, Rasio & Sills (1999) Gonzalez (1997); Fuhrmann, Pfeiffer & Bernkopf (1998); Gonzalez Gonzalez (1997); Fuhrmann, Pfeiffer & Bernkopf (1998); Gonzalez Perryman et al. (1997) Fuhrmann, Pfeiffer & Bernkopf (1998) Henry et al. (2000); Baliunas et al. (1997); Noyes et al. (1984) Fuhrmann, Pfeiffer & Bernkopf (1998); Ford, Rasio & Sills (1999) At this stage of the subject, our primary goal is to achieve detections of the planetary reflected light signal, thereby reaching step (2) in the above list – measuring Kp and ǫ0 and hence Rp for an assumed albedo p. With sufficient data we can measure epsilon(λ) and hence p(λ), but with present data the best we can do is to adopt p(λ) for various planetary atmospheric models and for each one measure or place an upper limit on the planet radius Rp /a. To optimize the sensitivity for detection, we restrict observations to gibbous phases of the planetary orbit, thus sacrificing information about the phase function g(α, λ). We elaborate these issues in more detail below. Radial velocity amplitude The radial velocity curve of the reflected light at orbital phase φ is Vp (φ) = Kp sin 2πφ, (1) The orbital phase at time t is given by φ = (t − T0 )/Porb , using the values of Porb and T0 given in Table 1. The mass ratio q = Mp /M⋆ is determined from the observed amplitude K⋆ of the star’s reflex orbital velocity via: q 1+q = K⋆ Porb sin i 2πa = 5.339 × 10−4 sin i  K⋆ 74.5 m s−1   M −1/3 ⋆ 1.3M⊙ (1998) (1998) Marcy, personal communication Marcy, personal communication Butler et al. (1999) Butler et al. (1999) (revised value derived in this paper) Butler et al. (1999) (revised value derived in this paper) (iv) g(α, λ), the phase function describing the dependence of the amount of light reflected toward the observer on the star-planet-observer angle α. (v) ∆, the velocity width of lines in the reflected starlight. This depends mainly on the star’s rotation in the frame of the orbiting planet, and to a lesser extent on the planet’s rotation and surface winds. 2.1 (1998) (1998) (. 2) The apparent radial velocity amplitude Kp of the planet about the system’s centre of mass is Kp = 2πa sin i sin i = 139 Porb 1 + q (1 + q)2/3  M⋆ 1.3M⊙ 1/3 km s−1 . (3) Kepler’s third law gives a = 0.0592  M⋆ 1.3M⊙ 1/3 AU. (4) The orbital inclination i is, according to the usual convention, the angle between the orbital angular momentum vector and the line of sight. For all but the lowest inclinations, the planet’s orbital velocity amplitude is considerably greater than the typical widths of the star’s photospheric absorption lines. Lines in the planet’s reflected-light spectrum should therefore be Doppler-shifted well clear of their stellar counterparts, allowing clean spectral separation for most of the orbit. 2.2 Orbital inclination Orbital inclinations i > 80◦ or so are ruled out by the absence of transits in high-precision photometry (Henry et al. 2000). The minimum inclination at which grazing transits occur is given by R⋆ + Rp . (5) a Low inclinations are also ruled out by the star’s chromospheric activity levels. The stellar radius and v sin i suggest an axial rotation period of 9 days. A 12 ± 2 day period was estimated by Baliunas et al. (1997) based on four Ca ii H and K flux measurements and the period-activity-colour relation of Noyes et al. (1984). A more comprehensive study by Henry et al. (2000), based on an additional 212 Ca ii H and K flux measurements in the 1996, 1997 and 1998 observing seasons shows the long-term average Ca ii H and K flux level to be almost identical to the original estimate by Baliunas et al. The rms scatter in the 30-day means is 4.7% of the long-term mean flux level. The 20% uncertainty in the 11.7-day period estimated from the long-term mean chromospheric flux is therefore dominated by the intrinsic scatter in of about 0.08 dex in the Noyes et al. (1984) calibration. cos imin = 4 A.Collier Cameron, K. Horne, A. Penny and C. Leigh Figure 2. The greyscale shows the prior joint probability density function (PDF) for projected orbital velocity Kp and the squared ratio (Rp /a)2 of the planet radius to the orbit radius, based on the measured system parameters. Darker shades in the greyscale denote greater probabilities of the planet having the corresponding combination of (Rp /a)2 and Kp . The one-dimensional projections of the PDF on to the two axes are shown as histograms. They show that the planet is most likely to have Kp ≃ 135 km s−1 and (Rp /a)2 ≃ 10−4 . The horizontal line shows the value of (Rp /a)2 for a 1 RJ up planet in the same orbit. Significantly faster axial rotation is needed to match the observed v sin i at orbital inclinations less than 60◦ , but is hard to reconcile with the observed chromospheric activity level. Henry et al. also point out that frequency analyses of the short-term variability in the 1996 and 1998 seasons yielded possible low-amplitude rotational modulation signals with periods of 11 and 19 days respectively, although both were classified as “weak”. These rather loose constraints on i require a Monte Carlo simulation to express our knowledge of the system parameters. For this we first generate a distribution of sin i values PCaII .v sin i (6) sin i = 2πR⋆ computed from randomly-chosen values of the stellar rotation period, radius and v sin i and rejecting those combinations that yielded sin i > 1. We assume gaussian distributions for the log of the rotation period, the measured stellar mass and radius, and the stellar reflex velocity K⋆ . High inclinations are then eliminated by the observed absence of eclipses. Low i are then given a reduced probability consistent with the measured Vrot sin i = 9.5 km s−1 and estimate of its rotation period Prot ∼ 12 ± 2d, which together favour a relatively high i. The lower histogram in Fig. 2 shows the resulting probability distribution for Kp (equivalent to i). Figure 3. A Monte Carlo simulation shows how the apparent stellar equatorial rotation speed ve,ref l in the starlight incident on the planet is expected to vary with the projected orbital velocity amplitude Kp . The system parameters are defined as in Section 2. The median rotational broadening of the reflected starlight is 8.4 km s−1 , somewhat less than the v sin i = 9.5 km s−1 of the direct starlight. 2.3 Rotational broadening The matched-filter analysis used to search for the reflectedlight signal requires an a priori estimate of the rotational broadening of the spectral-line profiles in starlight reflected from the planet’s atmosphere. The Earth-bound observer sees stellar aborption lines that are rotationally broadened by v⋆ sin i = 2πR⋆ sin i = 9.5 ± 0.8km s−1 . Prot,∗ (7) We assume that the planet orbits in the same direction as the stellar rotation and that the orbital plane and stellar equator are coincident. The light received by the orbiting planet, and subsequently reflected toward an observer in any direction, then has a rotational broadening vref l = 2πR⋆  1 Prot,∗ − 1 Porb  . (8) The Monte Carlo simulation gives the resulting relationship between the broadening vref l of the reflected starlight and Kp , as shown in Fig. 3. The distribution of the predicted rotational broadening of the reflected starlight has a mean 8.3 km s−1 , σ = 1.1 km s−1 , and median 8.4 km s−1 . The planet line widths are reduced at lower orbital inclinations, because the difference between the orbital and stellar rotation frequencies is reduced. A search for starlight reflected from υ And’s innermost planet 5 1 0.8 0.6 Flux 0.4 0.2 0 0 0.2 0.4 0.6 Orbital phase 0.8 1 Figure 4. Brightness variation of model planet with orbital phase assuming a Venus-like phase function (solid) or a Lambert sphere (dashed) for orbit inclinations i = 0, 30, 60, 90◦ . In both cases, flux is measured relative to the value for illumination phase angle α = 0. The horizontal lines show the i = 0◦ case, while i = 90◦ gives the greatest amplitude. 2.4 Albedo spectrum The quantity that we observe is the ratio fp /f⋆ of the spectral-line strengths in the planet’s light to those in the direct stellar spectrum as a function of wavelength: ǫ(α, λ) ≡ Rp2 fp (α, λ) = p(λ)g(α, λ) 2 = ǫ0 (λ)g(α, λ). f⋆ (λ) a (9) This depends on the fraction (Rp /2a)2 of the star’s light intercepted by the planet, and the fraction of this light that is reflected towards the observer. The monochromatic stellar flux illuminating the planet, orbiting at distance a, is Fincident (λ) = L⋆ (λ) . 4πa2 (10) The geometric albedo p(λ) is defined at α = 0 p(λ) = Freflected (0, λ) . Fincident (λ) (11) Here Freflected (0, λ) = πIreflected (0, λ) is the disc-averaged flux of the starlight reflected from the planet directly back toward the star, also measured at the planet’s surface. As the planet orbits the star, the observed flux varies with the star-planet-observer “phase angle” α as the product of the geometric albedo p(λ) and a “phase function” g(α, λ) which is normalised to g(0, λ) = 1. The flux received from the planet at Earth, is therefore fp (α, λ) = p(λ)g(α, λ)Fincident (λ) Rp2 . D2 (12) Since the stellar flux received at Earth, a distance D away, is f⋆ (λ) = L⋆ (λ) , 4πD2 eq. 9 follows. Figure 5. A Monte Carlo simulation shows how the gravitational acceleration g at the planet’s surface is expected to vary with the projected orbital velocity amplitude Kp . The system parameters are defined as in Section 2. For comparison, gJ up = 26.5 m s−2 . 2.5 Phase function The phase angle α depends on the orbital inclination i and the orbital phase φ (measured from transit or inferior conjunction): cos α = − sin i. cos 2πφ. (14) The observations measure ǫ(α, λ) over some range of orbital phases φ and hence phase angles α. However, the signal-to-noise ratio and orbital phase coverage of the observations is not yet adequate to define the shape of the phase function. Accordingly, current practice is to adopt a specific phase function in order to express the results in terms of the planet/star flux ratio that would be seen at phase angle zero: ǫ0 (λ) = p(λ) Rp2 . a2 (15) Since a is tightly constrained through Kepler’s third law, given the period P and the star mass p M⋆ , the measurements of ǫ0 (λ) measure the product Rp p(λ). In this paper we adopt an empirical phase function similar to that of Venus and Jupiter. The expressions given in Appendix D1 (eqs. D4 and and Eq. D6 for the phase function of a Lambert sphere and Venus respecively) yield the phase-dependent flux correction factors plotted in Fig. 4. (13) 2.6 Planet radius and surface gravity Because the surface gravity appears to play an important role in determining the height of the cloud deck and hence 6 A.Collier Cameron, K. Horne, A. Penny and C. Leigh the albedo, we plot the relationship between g and Kp derived from the Monte Carlo simulations, in Fig. 5. The planet’s mass can be determined once the orbital inclination is known. Our a priori knowledge of the radius, however, comes only from structural models for close-orbiting giant planets. In these models, the planet radius evolves with time and depends on the planet mass. For our purposes, theoretical mass-radius-age relations define a range of plausible radii at each possible value of the planet mass. The range of possible planet radii was computed specifically for υ And b by Guillot et al. (1997), allowing for uncertainties in the orbital inclination and hence the planet’s mass. Their radiativeconvective gas giant models predict a radius between 1.2 and 1.7 RJup for Mp = 0.6MJup , to between 1.15 and 1.3 RJup for Mp = 1.15MJup . When this mass-radius relation and its associated uncertainty are incorporated in the Monte Carlo model, we find that as the inclination decreases, the planet mass increases and the radius decreases. The lowest values of g are therefore found when the orbit is, as is most probable, nearly edge-on to the line of sight (Fig. 5). The median predicted value, g = 11.6 ms−2 , is close to the ∼ 10 m s−2 limit below which the Class V albedo models of Sudarsky, Burrows & Pinto (1999)) are applicable. If Class V models apply, we can use them to estimate the planet’s brightness. Viewed fully illuminated, the reflected light from a planet with (Rp /a)2 ≃ 10−4 and Class V albedo p ≃ 0.42 (cf. the Class V models of Sudarsky, Burrows & Pinto (1999)) is expected to be 24000 times fainter than the direct stellar spectrum in the visual band around 550 nm. At a phase angle of 60◦ , however, the reflected spectrum will be a factor two fainter than this (Fig.4). 3 OBSERVATIONS The strong dependence of the flux ratio on phase angle indicates that the best chance of a detection occurs shortly before and after superior conjunction. We observed υ And with the Utrecht Echelle Spectrograph on the 4.2-m William Herschel Telescope at the Roque de los Muchachos Observatory on La Palma, on the nights of 2000 Oct 10, Nov 6 and Nov 7. These nights were chosen to cover orbital phase ranges when the planet is on the far side of the star and well-illuminated, and its spectral signature is Doppler-shifted well clear of the wings of the stellar profile. The third night was partially affected by drifting cirrus cloud. Two further nights were also lost to cloud. The detector was an array of two EEV CCDs, each with 2048 × 4096 13.5-µm pixels. The CCD was centred at 459.6 nm in order 124 of the 31 g mm−1 echelle grating, giving complete wavelength coverage from 381.3 nm to 675.8 nm with minimal vignetting. The average pixel spacing was close to 1.6 km s−1 , and the full width at half maximum intensity of the thorium-argon arc calibration spectra was 3.5 pixels, giving an effective resolving power R = 53000. Table 2 lists the journal of observations for the three nights that contributed to the analysis presented in this paper. The stellar spectra were exposed for between 300 and 500 seconds, depending upon seeing, in order to expose the CCD to a peak count of 40000 ADU per pixel in the brightest parts of the image. A 450-s exposure yielded about 1.2 × 106 electrons per pixel step in wavelength in the brightest orders in typical (1 arcsec) seeing after extraction. We achieve this with the help of an autoguider procedure, which improves efficiency in good seeing by trailing the stellar image up and down the slit by ±2 arcsec during the exposure to accumulate the maximum S:N per frame attainable without risk of saturation. Note that the 450-s exposure time compares favourably with the 53-s readout time for the dual EEV CCD in terms of observing efficiency – the fraction of the time spent collecting photons is above 90 percent. Following extraction, the S:N in the continuum of the brightest orders is typically 1000 per pixel. 4 SPECTRUM EXTRACTION One-dimensional spectra were extracted from the CCD frames using an automated pipeline reduction system built around the Starlink ECHOMOP and FIGARO packages. Nightly flat-field frames were summed from 50 to 100 frames taken at the start and end of each night, using an algorithm that identified and rejected cosmic rays and other non-repeatable defects by comparing successive frames. The nightly flat fields were then added to make a master flat-field for the entire year’s observations. The initial tracing of the echelle orders on the CCD frames was performed manually on the spectrum of υ And itself, using exposures taken for this purpose without dithering the star up and down the slit. The automated extraction procedure then subtracted the bias from each frame, cropped the frame, determined the form and location of the stellar profile on each image relative to the trace, subtracted a linear fit to the scattered-light background across the spatial profile, and performed an optimal (profile and inverse variance-weighted) extraction (Horne 1986; Marsh 1989) of the orders across the full spatial extent of the object-plussky region. Flat-field balance factors were applied in the process. In all, 62 orders were extracted from each exposure. The blue CCD recorded orders 148 in the blue to 125 in the red, giving full spectral coverage from 380.7 to 461.0 nm with considerable wavelength overlap between adjacent orders. Orders 122 to 85 were recorded on the red CCD, covering the range 461.9 to 677.3 nm with good overlap. Orders 123 and 124 fell on the gap between the CCDs, but the loss in wavelength coverage was minimal. 5 STARLIGHT-SUBTRACTED VELOCITY-PHASE MAPS The maximum expected flux of starlight scattered from υ And b is, as we have seen, of order one part in 24000 of the flux received directly from υ And itself. In order to detect the planet signal, we first subtract the direct stellar component from the observed spectrum, leaving the planet signal embedded in the residual noise pattern. To achieve this, we make a model of the direct starlight by aligning and summing all the spectra of the target from all nights of the run. The contribution of the planet to this summed spectrum is blurred by the planet’s orbital motion, so that to first order the planet’s spectrum is eliminated A search for starlight reflected from υ And’s innermost planet 7 Table 2. Journal of observations. The UTC mid-times and orbital phases are shown for the first and last groups of four spectra secured on each night of observation. The number of groups is given in the final column. UTC start Phase UTC End Phase Number of groups 2001 Oct 09 22:27:51 2001 Nov 06 19:43:04 2001 Nov 07 19:25:48 0.296 0.335 0.549 2001 Oct 10 06:20:26 2001 Nov 07 04:27:15 2001 Nov 08 04:29:30 0.366 0.413 0.630 18 16 15 from the composite “template” spectrum. We then scale and distort this template spectrum to give the closest possible fit to each individual spectrum, using the spline-modulated Taylor-series method described in Appendix A. When the aligned template is subtracted from each spectrum, we find that the residual spectrum contains a spatially fixed but temporally varying pattern of ripples, superimposed on random noise. The ripples are attributable to a combination of time-dependent changes in the detector’s sensitivity, thermal flexure in the spectrograph causing faint ghost reflections to shift slightly on the detector, and changes in the strength and velocity of telluric-line absorption features in some orders. We map the form of the ripple pattern using principal-component analysis, as described in Appendix B. After subtracting this map, we are left with a planet signal consisting of faint Doppler-shifted copies of thousands of stellar absorption lines, deeply buried in noise. The positions and identifications of most of these thousands of lines are well-documented. We use the Vienna atomic-line database (VALD; Kupka et al. 1999) to compile a list of the wavelengths and relative central depths of the 3450 strongest lines falling in the observed wavelength range, for a model atmosphere appropriate to the spectral type and surface gravity of the star. We then use least-squares deconvolution (LSD; see Appendix C) to compute a composite profile which, when convolved with the line pattern, yields an optimal fit to the residual noise spectrum. This procedure is similar to cross-correlation in terms of the gain in S:N from the weighted summation of thousands of line profiles, but has the additional advantage of eliminating sidelobes caused by crosstalk with neighbouring lines. The resulting devonvolved residual profiles are presented in greyscale form as a two-dimensional “velocityphase” map in Fig. 6. This representation of the data shows a strong pattern of distortions in the residual stellar profiles within ±20 km s−1 of the stellar velocity. These undulations in the deconvolved profiles appear to be caused by high-order mismatches in the spectrum subtraction procedure. Their amplitude is typically a few parts in 104 of the average continuum level. They vary too rapidly during the night to be attributable to, e.g. stellar surface features causing time-dependent distortion of the stellar rotation profile. A more likely explanation is that the spatial profile produced by the telescope dithering procedure was not exactly repeatable from one frame to the next. Given that the slit image is slightly tilted with respect to the columns of the detector, this could give rise to small changes in the detailed shape of the line profile from one exposure to the next. Fortunately these ripples only affect a range of velocities at which the planet signature would in any case be indistinguishable from that of the star. We deliberately avoided observing between phases 0.45 and 0.55 for this reason. Figure 6. Velocity-phase map of deconvolved residual profiles derived from WHT spectra secured on 2000 Oct 09, Nov 06 and Nov 07. The line weights in the deconvolution were defined assuming a grey albedo spectrum. The greyscale runs from black at −10−4 times the mean stellar continuum level, to white at +10−4 . The velocity scale is in the reference frame of the star. Outside the residual line profile, we would expect to see a planet signature as a dark streak following a sinusoidal path that crosses the profile from positive to negative velocity at phase 0.5. No obvious planet signature is, however, visible in Fig. 6. We verified that a faint planetary signal is preserved through the analysis in the presence of realistic noise levels, by adding a simulated planetary signal to the observed spectra and performing the same sequence of operations described above. The simulation procedure simply consisted of shifting and scaling the spectrum of υ And according to the phase function, co-multiplying it by an appropriate geometric albedo spectrum, and adding it to the observed data. To ensure a strong signal we assumed the planet to have a radius 2.0 times that of Jupiter, giving (Rp /a)2 = 2.57 × 10−4 . Initially we used a wavelength-independent geometric albedo with p = 0.42, which was chosen to match the continuum albedo of a Class V model unaffected by alkali-metal absorption. When viewed at zero phase angle, the planet- 8 A.Collier Cameron, K. Horne, A. Penny and C. Leigh Figure 7. As for Fig. 6, with a simulated planet signature added at an orbital inclination of 80◦ . The model planet has a grey albedo spectrum with p = 0.42 and a radius twice that of Jupiter, and its signature appears as a dark streak crossing from positive to negative velocity at phase 0.5. Figure 8. As for Fig. 6, but showing the matched filter function used to measure (Rp /a)2 at Kp = 137 km s−1 . The grey scale runs from -0.05 (black) to 0.01(white). Scaling this function to fit Fig. 6 (and others like it) yields the strength of the planet signal for a given value of Kp . to-star flux ratio should thus be ǫ0 = fp /f∗ = 1.08 × 10−4 . We assumed an orbital inclination of 80◦ , giving an orbital velocity amplitude Kp = 137 km s−1 and a rotational broadening of the reflected starlight, ve,ref l = 8 ± 1 km s−1 . We are therefore justified in using the spectrum of υ And itself, without any modification of the rotational broadening, to approximate the reflected-light spectrum. The injected planet signal is clearly visible as a dark streak in Fig. 7. Given that no similar signal is easily seen in Fig 6, we can conclude that the planet in υ And is fainter than this one. signature is nearly stationary in this part of the orbit, some of the signal will be removed along with the stellar profile if many observations are made in this part of the orbit. The planet is detectable on a given night because its velocity changes while those of the stellar lines do not. Thus observations at quadrature are less helpful. In the present dataset, few observations were made near quadrature so this problem does not arise. The travelling gaussian in this image has nonetheless been modified to give an even closer match to the observed planet signal, by subtracting its own flux-weighted average. This procedure mimics the attenuation of the planet signal that occurs when the template spectrum is subtracted from the individual spectra. The main effect of the template subtraction is to produce the bright vertical zones seen centred at v ≃ −70 km s−1 and v ≃ +100 km s−1 in Fig. 8. 6 EXTRACTING THE PLANET SIGNATURE The form of the expected planet signature in the velocityphase diagram can be represented accurately as a travelling Gaussian of characteristic width ∆vp whose velocity varies sinusoidally around the orbit and whose strength is modulated according to the phase function g(α, λ) (see Appendix D1). The value of ∆vp is adjusted to match the expected rotational broadening of the reflected starlight. The amplitude Kp of the velocity variation and the form of the phase function both depend on the orbital inclination i. Fig. 8 shows the form of this model signal at an orbital inclination i = 80◦ . The weakening of the simulated planetary signature near quadrature (phases 0.25 and 0.75) is mostly the effect of the phase function. In some cases the signal may be further attenuated near quadrature by the way in which the templates are computed: since the planet 7 PROBABILITY MAPS FOR KP AND RP /A We use a sequence of such models for different values of Kp as matched filters to measure the strengths and velocity amplitudes of possible faint planet signals of the expected form in the velocity-phase maps. At each of a sequence of trial values of Kp we construct a model Hij (Kp ) like that shown in Fig. 8 and scale it to give an optimal fit to the residual velocity-phase map Dij using the methods presented in Appendix D. For an assumed albedo spectrum p(λ) and orbital velocity amplitude Kp , the fit of the matched filter to the data measures (Rp /a)2 . A search for starlight reflected from υ And’s innermost planet Figure 9. The upper panel shows the optimal scaling factor (Rp /a)2 plotted against orbital velocity amplitude Kp , assuming a grey albedo spectrum with p = 0.42. The lower panel shows the associated reduction in χ2 measured relative to the fit obtained in the absence of any planet signal, i.e. for (Rp /a)2 = 0. Note that only positive values of (Rp /a)2 are physically plausible, so the middle peak in the ∆χ2 plot is unambiguously a noise feature. The badness-of-fit statistic χ2 gives a relative measure of the probability that any signal detected could have arisen by chance, and is quantified by χ2 = X (Dij − (Rp /a)2 Hij (Kp ))2 i,j 2 σij . (16) 2 Here σij is the variance associated with Dij , computed from the photon statistics of the original image and propagated through the deconvolution (see Appendix C). In Fig. 9 we plot the optimum scale factor and the associated improvement in χ2 measured relative to the noplanet model. For a noise-dominated signal, the estimated value of (Rp /a)2 will be negative about as often as it is positive, and Fig. 9 bears this out. Three candidate peaks are seen in the lower panel, of which only the first and the third correspond to positive (and therefore physically plausible) reflected-light fluxes. The third (positive) peak gives the greatest improvement in χ2 and is therefore the most probable, but only by a narrow margin. The relative probabilities of the fits to the data for different values of the free parameters Rp /a and Kp , given by P (Kp , (Rp /a)2 ) = exp(−χ2 /2), (17) 9 Figure 10. Relative probability map of model parameters Kp and log(ǫ0 /p) = log(Rp /a)2 , derived from the WHT/UES observations of υ And, assuming a grey albedo spectrum with p = 0.42. The greyscale denotes the probability relative to the best-fit model, increasing from 0 for white to 1 for black. From bottom to top, the contours show the 1, 2, 3 and 4σ upper limits limits on the signal strength derived from the bootstrap trials The uppermost contour, for example, represents the value of log(Rp /a)2 that was only exceeded in 3 of 3000 bootstrap trials at each Kp . It gives a robust empirical estimate of the upper limit on the planet radius allowed by the data for the grey albedo model with p = 0.42 assumed here. are shown in Fig. 10. To make this figure, a planet signal pattern as in Fig. 8 was computed for many different Kp , and each one was scaled by (Rp /a)2 to fit the residuals in Fig 6. The greyscale is defined such that white represents the probability of a model fit where no planet signal is present, while black represents the most probable solution in the map. The curves give 1, 2, 3 and 4σ detection thresholds derived from the bootstrap procedure discussed below. To set an upper limit on the strength of the planet signal, or to assess the likelihood that a candidate detection is spurious, we need to compute the probability of obtaining such an improvement in χ2 by chance alone. In principle this could be done using the χ2 distribution for 2 degrees of freedom. In practice, however, the distribution of pixel values in the deconvolved difference profiles has extended non-gaussian tails that demand a more cautious approach. Rather than relying solely on formal variances derived from photon statistics, we use a “bootstrap” procedure to construct empirical distributions for confidence testing, using the data themselves. The bootstrap procedure, detailed in Appendix E, interchanges at random the order of the nights, and rearranges at random the order of spectra in each night, thereby scrambling planet signals while retain- 10 A.Collier Cameron, K. Horne, A. Penny and C. Leigh Table 3. Upper limits on planet dimensions for various albedo models. The upper limits are quoted for an assumed Kp ≃ 135 km s−1 , near the peak of the prior probability distribution for Kp . Note that the results for the grey albedo model are quoted for unit geometric albedo. To obtain the radii for grey models with other values of p, the radii in Column 4 should be divided √ by p. False-alarm probability (Rp /a)2 Rp /RJ up Grey (p = 1.00) 0.1% 1.0% 4.6% 0.58E-04 0.45E-04 0.34E-04 0.98 0.86 0.75 Class V 0.1% 1.0% 4.6% 1.42E-04 1.11E-04 0.87E-04 1.53 1.35 1.19 Class IV (Isolated) 0.1% 1.0% 4.6% 3.02E-04 2.20E-04 1.68E-04 2.23 1.90 1.66 Albedo model Figure 11. Relative probability map of model parameters Kp and log(ǫ0 /p) = log(Rp /a)2 for a simulated planet signature with grey albedo p = 0.42, Kp = 137 km s−1 , and Rp = 2RJ up . The greyscale and contours are defined as in Fig. 10. The synthetic planet signature is detected well above the 4σ upper limit on the strength of noise features in the absence of a planet signal. ing the same pattern of systematic errors, phase sampling, and statistical noise as in the actual data. The 68.4%, 95.4%,99.0% and 99.9% percentage points of the resulting bootstrap distribution of (Rp /a)2 at each Kp are shown as contours in Fig. 10 and Fig. 11. From bottom to top, the contours give the 1, 2, 3, and 4σ bootstrap upper limits on the strength of the planet signal. They represent the signal strengths at which spurious detections occur with 32, 5, 1, and 0.1% false alarm probability respectively, for each fixed value of Kp . 7.1 Grey albedo model At the most probable values around Kp ≃ 135 km s−1 , the grey albedo model yields a 0.1% bootstrap upper limit on the planet/star flux ratio ǫ0 < 5.84 × 10−5 . Adopting a grey albedo model with unit geometric albedo, we find that the 0.1%, 1.0% and 4.6% upper limits on ǫ0 at this velocity correspond to upper limits on the planet radius as listed in Table 3. If the orbital inclination is lower, the planet radius is less strongly constrained. Two possible planet signals of comparable likelihood are seen. Their properties are listed in Table 4. The stronger, at Kp ≃ 132 km s−1 (i ≃ 80◦ ), yields an improvement ∆χ2 = 13.61 over the model fit obtained assuming no planet signal is present. If this feature represents a genuine planet signal, its velocity amplitude Kp = 132 km s−1 implies a planet mass Mp = 0.74MJ up and (for p = 0.42) a planet radius Rp = 1.24 ± 0.17RJ up . The weaker candidate has Kp ≃ 54 km s−1 (i = 22◦ ) and ∆χ2 = 11.86, and gives a larger planet radius. We used the bootstrap simulations to determine the probability that a spurious detection with ∆χ2 > 13.61 could be produced by a chance alignment of noise features in the absence of a genuine planet signal. It is important to note that the location of the “blob” between the 2σ and 3σ bootstrap contours in Fig. 10 does not imply a false-alarm probability of ∼ 3%. These contours give the false-alarm probability only if the value of Kp is known in advance, which is not the case here. The true false-alarm probability is greater, being the frequency with which spurious peaks at any plausible value of Kp can exceed the ∆χ2 of the candidate. If we assume that all values of Kp are equally likely in the range 44 km s−1 < Kp < 152 km s−1 over which our phase coverage allows signals to be detected, the false-alarm probability is found to be 26% via the method described in Appendix E. In practice, however, we are more likely to be fooled into believing that a candidate detection near the peak of the a priori probability distribution for Kp is genuine, than would be the case if the candidate appeared at a velocity that was physically implausible given our existing knowledge of the system parameters. We can therefore refine the search range in Kp using our knowledge of the a priori probability distribution for Kp , via the method presented in Appendix E. The last two columns of Table 4 show clearly that, while the unweighted false alarm probabilities for the features near Kp = 132 and Kp = 54 km s−1 are comparable, the latter is almost certainly spurious. The false-alarm probability for the Kp = 132 km s−1 feature drops to 9.4% when prior knowledge of Kp is accounted for, while that for the 54 km s−1 feature approaches unity. For comparison, we show in Fig 11 a probability map derived by subtracting Fig. 6 from Fig. 7 to isolate the injected planet signal, then performing the matched-filter analysis. The injected planet is clearly detected as a compact, dark feature with its correct amplitude log(Rp /a)2 = 2.57×10−4 and orbital velocity amplitude velocity Kp = 137 km s−1 . This most probable combination of orbital veloc- A search for starlight reflected from υ And’s innermost planet 11 Table 4. Orbital velocity amplitude, signal strength, planet radius, ∆χ2 and false-alarm probabilities (FAP) for possible planet signals. The FAP is computed for two different priors. The first prior is uniform in Kp from 44 to 152 km s−1 . The second weights Kp in proportion to the prior probability based on Caii H & K emission and v sin i constraints and the absence of eclipses. In the latter case, the false-alarm probability for the Kp = 132 km s−1 candidate is found to be between 9 and 10 percent. Kp (km s−1 ) (Rp /a)2 (×10−4 ) Rp /RJ up ∆χ2 FAP (uniform weight) FAP (Kp prior) Grey (p = 1.00) 54 132 0.87 ± 0.25 0.39 ± 0.11 1.20 ± 0.17 0.80 ± 0.11 11.93 13.61 0.324 0.259 0.975 0.094 Class V 50 132 2.08 ± 0.64 1.09 ± 0.28 1.85 ± 0.29 1.34 ± 0.17 10.39 14.61 0.423 0.240 0.983 0.087 Class IV 49 132 3.39 ± 1.31 1.73 ± 0.52 2.36 ± 0.47 1.68 ± 0.25 6.66 10.78 0.534 0.285 0.995 0.099 Albedo model Figure 12. As for Fig. 6, showing the time-series of deconvolved profiles derived from the original observations assuming the albedo spectrum to be that of a Class V roaster. ity and planet radius yields an improvement ∆χ2 = 98.6 with respect to the value obtained assuming no planet is present (i.e. Rp = 0). This is far greater than the greatest ∆χ2 = 55.3 produced at any value of Kp in any of the bootstrap trials. The false-alarm probability for this simulated signal is substantially less than one part in 3000, and the “detection” is secure. We note that the calibration of the (Rp /a)2 scale in Fig. 10 depends on the value p = 0.42 assumed for the geometric albedo. In the next sections, we explore the relative ability of non-grey albedo models to fit the data. Figure 13. As for Fig.9, showing the optimal scaling factor (Rp /a)2 and the associated improvement in χ2 plotted against orbital velocity amplitude Kp , assuming a Class V albedo spectrum. 7.2 Class V roaster model The “Class V roaster” is the most highly reflective of the models computed by Sudarsky, Burrows & Pinto (1999). This model is characteristic of planets with Tef f ≥ 1500 K and/or surface gravities lower than ∼ 10 m s−2 , and has a silicate cloud deck located high enough in the atmosphere that the overlying column density of gaseous alkali metals is low, allowing a substantial fraction of incoming photons at most optical wavelengths to be scattered back into space. There remains, however, a substantial absorption feature 12 A.Collier Cameron, K. Horne, A. Penny and C. Leigh Figure 14. Relative probability map of model parameters Kp and log(ǫ0 /p) = log(Rp /a)2 , derived from the WHT/UES observations of υ And, assuming the albedo spectrum to be that of a Class V roaster. The greyscale and contours are defined as in Fig. 10. around the Na I D lines, as shown in Fig. 1. If υ And b has a mass and radius close to the values at the peak of the prior probability distribution, its surface gravity is close to the critical limit (Fig. 5). We carried out the deconvolution using the same line list as for the grey model, but with the line strengths attenuated using the Class V albedo spectrum (see Appendix C). We back-projected the resulting time series of deconvolved profiles (Fig. 12) as described above. We calibrated the signal strength as described in Appendix D4, by injecting an artificial planet signature consisting of the spectrum of υ And, attenuated by the Class V albedo spectrum and scaled to the signal strength expected for a planet with Rp = 2Rjup . For the observations we used ∆vp = 8 km s−1 which again yielded the best fit, with scale factors and improvements in χ2 as plotted in Fig. 13. The probability map for the observed signals is shown in Fig. 14. The form of the Class V probability map is similar to that encountered for the grey albedo spectrum. The resulting upper limits on the planet radius are listed in Table 3. The local probability maxima near Kp = 52 and 132 km s−1 are also present with this albedo model. As in the grey case, the feature at Kp = 132 is the most probable, this time by a slightly wider margin. Their signal strengths and false-alarm probabilities are given in Table 4. 7.3 Isolated Class IV model The Class IV roaster models of Sudarsky, Burrows & Pinto (1999) have a more deeply-buried cloud deck than the Class Figure 15. As for Fig. 6, showing the time-series of deconvolved profiles derived from the original observations assuming the albedo spectrum to be that of an “isolated” Class IV roaster. V models. The resonance lines of Na I and KI are strongly saturated, with broad wings due to collisions with H2 extending over much of the optical spectrum (Fig. 1). We used the procedures described above to deconvolve and back-project the data assuming an “isolated” Class IV spectrum. Although this model does not take full account of the effects of irradiation of the atmospheric temperaturepressure structure, it is a useful compromise between the Class V models and the very low albedos found with irradiated Class IV models. The resulting time-series of deconvolved spectra (Fig. 15) is noisier than the Class V and grey-albedo versions, because lines redward of 500 nm contribute little to the deconvolution. Consequently, the derived upper limits on the planet radius (Table 3) are higher than for the Class V albedo model. The plots of (Rp /a)2 and ∆χ2 versus Kp show the same mix of positive and negative signal amplitudes encountered with the grey and class V spectral models. The back-projected probability map (Fig.17) nonetheless shows the Kp = 132 km s−1 probability maximum to be present and again stronger than the 49 km/sec feature. The fit to the data is, however, poorer than in the grey and Class V cases, giving an improvement of only ∆χ2 = 10.8 over the no-planet hypothesis and a planetary radius Rp = 1.68±0.25RJ up (Table 4). This is 25% larger than the radius estimated with the Class V albedo model. The bootstrap test gives a false-alarm probability of 10%. This relatively large radius and poor fit appear to arise from a mismatch between the wavelength dependence of the putative planet signal and the Class IV model. We tested this notion by deconvolving and back-projecting a synthetic A search for starlight reflected from υ And’s innermost planet Figure 16. As for Fig.9, showing the optimal scaling factor (Rp /a)2 and the associated improvement in χ2 plotted against orbital velocity amplitude Kp , assuming an isolated Class IV albedo spectrum. Class V model using the Class IV line weights and flux calibration. The radius of the planet was over-estimated by 27% in this experiment. This occurs because the least-squares deconvolution is attempting to fit a set of lines at all wavelengths, with a deconvolution mask in which the lines at red wavelengths are too weak to fit the data properly. The best compromise is achieved (in the least-squares sense) by overestimating the depth of the deconvolved profile in order to boost the strengths of the weakened lines. We conclude that the faint signals that contribute to the Kp = 132 km s−1 feature are distributed in wavelength in a way that resembles more closely a grey or a Class V spectrum than a strongly-absorbed Class IV spectrum in which only blue light is present. 7.4 Plausibility of candidate signals Comparison of Figs.10, 14 and 17 with the joint prior probability distribution in Fig. 2 shows that the feature near Kp = 50 km s−1 yields a value of (Rp /a)2 that is considerably greater than the radius predicted by Guillot et al. (1997). The feature at Kp = 132 ± 3 km s−1 , however, has ǫ0 = 4.6 × 10−5 . A plausible geometric albedo p = 0.42 then gives (Rp /a)2 = (1.09 ± 0.30) × 10−4 , which places it very close to the peak of the joint prior probability distribution. We probed the rotational broadening of the candidate features by performing a series of backprojections us- 13 Figure 17. Relative probability map of model parameters Kp and log(ǫ0 /p) = log(Rp /a)2 , derived from the WHT/UES observations of υ And, assuming the albedo spectrum to be that of an “isolated” Class IV roaster. The greyscale and contours are defined as in Fig. 10. ing values of the gaussian width parameter in the range 6.4 < ∆vp < 10 km s−1 , corresponding to 0 < ve,ref l < 11.3 km s−1 . The feature at 132 km s−1 yields the greatest improvement in χ2 when ∆vp = 8 ± 1 km s−1 , or ve,ref l = 7.0 ± 2 km s−1 . The 52 km s−1 feature, however, grows less significant as ve,ref l is decreased to the (retrograde) 5 km s−1 or so expected at the corresponding orbital inclination. The rotational broadening of the 132 km s−1 feature is therefore consistent with the value predicted in Fig. 3. We explored the region of parameter space around both features, by carrying out back-projections on a grid of orbital periods and epochs of zero phase around the values given in Table 1, taking the epoch of transit to be near the middle of the data train to ensure that the period and epoch were uncorrelated. The local minimum in χ2 was found to be centred at Kp = 133 ± 3 km s−1 , P = 4.621 ± 0.005 d, and T0 = 2451853.770 ± 0.025. The radial-velocity ephemeris gives P = 4.61707 ± 0.00003 d and T0 = 2451853.777 ± 0.012, well within the 1σ error bars derived from the backprojection. The 52 km s−1 feature, on the other hand, is part of a more extended structure whose peak is located at Kp = 67 ± 3 km s−1 , P = 4.651 ± 0.006 d, and T0 = 2451853.81 ± 0.02. This bears little relation to the radialvelocity solution, suggesting a spurious origin. We conclude that the candidate feature at Kp = 132 km s−1 corresponds to a local minimum of χ2 with respect to the parameters of interest, whose rotational broadening, phasing and velocity amplitude are consistent with those expected of a reflected-light signature from a planet whose orbital inclination is near the most probable value. We can- 14 A.Collier Cameron, K. Horne, A. Penny and C. Leigh not easily dismiss this feature as being merely an extension of a larger, spurious noise-induced structure. 8 of PPARC. We thank the referee for raising several useful points of clarification. CONCLUSIONS The observations of υ And presented here rule out radii for √ the innermost planet greater than Rp > 0.98Rjup / p with 0.1% false-alarm probability if a grey albedo spectrum is assumed. We derive upper limits also with 0.1% false-alarm probability on the planet’s radius of 1.53 Rjup assuming the albedo spectrum of a Class V roaster (low gravity, high cloud deck, Sudarsky, Burrows & Pinto 1999), or 2.23 Rjup for an isolated Class IV model with saturated Na I D absorption from 550 nm to the red limit of our spectra. We cannot, however, rule out the possibility that the planet has an albedo spectrum as bright as a Class V roaster (Sudarsky, Burrows & Pinto 1999) if its radius is comparable to that of HD 209458 b. The evidence for a reflected-light signature in the observations reported here, at a projected orbital velocity amplitude Kp = 132 ± 3 km s−1 , is marginal but encouraging. If the orbital inclination is as close to edge-on as we infer from previously-measured system parameters, the mass and radius of the planet yield a surface gravity low enough for a Class V atmosphere to be plausible. If future observations confirm this signal, it gives a radius Rp = 1.34 ± 0.17RJ up if a Class V spectrum is assumed. This is comparable to the radius of HD 209458 b inferred from HST transit photometry (Charbonneau et al. 2000). We have explored the spectral properties of the candidate detection via the relative probabilities of the fits to the data. The Class V roaster model gives a slightly better fit to the data than a grey albedo model. Both yield planet radii and orbital velocity amplitudes close to the peak of the prior probablility density function. The isolated Class IV albedo model gives a significantly poorer fit to the data. While the possible detection discussed here gives believable physical properties for the innermost planet, it remains too marginal for us to claim a bona fide detection of reflected starlight. Our bootstrap analysis suggests a 7% to 10% likelihood that the feature could be a spurious noise feature. There is clearly a strong case for a deeper reflected-light search to be made in the υ And system. ACKNOWLEDGMENTS We thank David Sudarsky and Adam Burrows for providing us with listings of their Class IV and Class V albedo models, and Geoff Marcy for supplying us with his most recent orbital ephemeris and radial velocity data for υ And b. This paper is based on observations made with the 4.2-m William Herschel Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. We are indebted to the support staff at the ING for their assistance, and in particular Ian Skillen for doing the tricky telescope scheduling, and John Telting for considerable effort expended in implementing the telescope dithering procedure used in these observations. We acknowledge the support software and data analysis facilities provided by the Starlink Project which is run by CCLRC on behalf REFERENCES Baliunas S. L., Henry G. W., Donahue R. A., Fekel F. C., Soon W. H., 1997, ApJ, 474, L119 Brown T. M., Charbonneau D., Gilliland R. L., Noyes R. W., Burrows A., 2001, ApJ, 552, 699 Burrows A. et al., 1997, ApJ, 491, 856, Provided by the NASA Astrophysics Data System Burrows A., Guillot T., Hubbard W. B., Marley M. S., Saumon D., Lunine J. I., Sudarsky D., 2000, ApJ, 534, L97 Butler R. P., Marcy G. W., Fischer D. A., Brown T. M., Contos A. R., Korzennik S. G., Nisenson P., Noyes R. W., 1999, ApJ, 526, 916 Charbonneau D., Noyes R. W., Korzennik S., Nisenson P., Jha S., Vogt S. S., Kibrick R. I., 1999, ApJ, 522, L145 Charbonneau D., Brown T. M., Latham D. W., Mayor M., 2000, ApJ, 529, L45 Collier Cameron A., Horne K. D., Penny A. J., James D. J., 1999, Nat, 402, 751 Collier Cameron A., Horne K., James D. J., Penny A. J., Semel M., 2001, in Penny A., Artymowicz P., Lagrange A.-M., Russell S., eds, IAU Symp. 202: Planetary systems in the Universe. ASP Conference Series, San Francisco, In press: astro-ph0012186 Donati J.-F., Semel M., Carter B., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658 Ford E. B., Rasio F. A., Sills A., 1999, ApJ, 514, 411 Fuhrmann K., Pfeiffer M. J., Bernkopf J., 1998, A&A, 336, 942 Gonzalez G., 1997, MNRAS, 285, 403 Gonzalez G., 1998, A&A, 334, 221 Guillot T., Burrows A., Hubbard W. B., Lunine J. I., Saumon D., 1996, ApJ, 459, L35 Guillot T., Saumon D., Burrows A., Hubbard W. B., Lunine J. I., Marley M. S., Freedman R. S., 1997, in Cosmovici C. V., Bowyer S., Wertheimer D., eds, IAU Colloquium 161: Astronomical and biochemical origins and the search for life in the Universe. Editrice Compositori, Bologna, p. 343 Henry G. W., Baliunas S. L., Donahue R. A., Fekel F. C., Soon W., 2000, ApJ, 531, 415 Hilton J. L., 1992, in Seidelmann P. K., ed, Explanatory supplement to the Astronomical Almanac. University Science Books, Mill Valley, CA, p. 383 Horne K. D., 1986, PASP, 98, 609 Hovenier J. W., Hage J. I., 1989, A&A, 214, 391 Kupka F., Piskunov N., Ryabchikova T. A., Stempels H. C., Weiss W. W., 1999, A&AS, 138, 119 Marcy G. W., Butler R. P., 1998, ARA&A, 36, 57 Marley M. S., Gelino C., Stephens D., Lunine J. I., Freedman R., Sharp C., 1999, ApJ, 513, 879 Marsh T. R., 1989, PASP, 101, 1032 Mayor M., Queloz D., 1995, Nat, 378, 355 Noyes R. W., Hartmann L., Baliunas S. L., Duncan D. K., Vaughan A. H., 1984, ApJ, 279, 763 Perryman M. A. C. et al., 1997, A&A, 323, L49 Press W. H., Flannery B. P., Teukolsky S. A., Vetterling W. T., 1992, Numerical Recipes: The Art of Scientific Computing (2nd edition). Cambridge University Press, Cambridge Saumon D., Hubbard W. B., Burrows A., Guillot T., Lunine J. I., Chabrier G., 1996, ApJ, 460, 993 Seager S., Sasselov D. D., 1998, ApJ, 502, L157 Sudarsky D., Burrows A., Pinto P., 1999, ApJ, 538, 885 A search for starlight reflected from υ And’s innermost planet APPENDIX A: REMOVAL OF DIRECT STARLIGHT In this appendix we describe the procedure we use to subtract the stellar spectrum without introducing substantial noise and also without subtracting the planet signal. The intrinsic spectrum of the star was modelled as the co-aligned sum of all the spectra of υ And secured during the run. We began by removing cosmic-ray hits from the individual frames. Within each night’s data, each frame was divided by its predecessor and the result was median smoothed with a box covering 21 pixels in the dispersion direction and 5 adjacent orders. This smoothed frame was used to scale the predecessor frame, which was then subtracted from the frame under consideration. The difference frame was divided by the square root of the frame under consideration, to yield deviates in units of the local Poissonian noise amplitude. All pixels whose deviates were more than 7 sigma in excess of their counterparts in the scaled predecessor frame were assumed to be cosmic ray hits. Their values were replaced with values taken from the corresponding pixels in the scaled predecessor frame. The result was re-scaled and added to the scaled predecessor frame to restore the original frame with all major cosmic-ray hits removed. The cleaned frames were then shifted to co-align the stellar absorption lines to the nearest integer pixel. Finally the cleaned frames were summed to make the template frame (step 1 in Fig. A1). The individual cleaned frames were co-added in groups of four for the next stage in the procedure, which involved scaling, shifting and blurring (or sharpening) the template frame prior to subtraction from each group. The first and second derivatives of the template frame were computed from the template values tj in each order as a function of pixel number xj in the dispersion direction: t′j = tj+1 − tj−1 xj+1 − xj−1 (A1) and t′′j = 2 (tj+1 − tj )/(xj+1 − xj ) − (tj − tj−1 )/(xj − xj−1 ) .(A2) xj+1 − xj−1 The third and fourth derivatives were computed by applying the same operations to t′′ , for use in correcting frame-toframe changes in the higher moments of the PSF and seeing profile. The vector f of all observed spectrum values within each extracted echelle order was then modelled by scaling the template t along each order, the scale factors varying as a function of position x approximated by a 34-knot leastsquares spline. The use of such a high-order spline was demanded by vignetting near the edges of each order, which produced an abrupt and time-variable change in the slope of the stellar continuum towards the end of each order. The derivative frames were scaled in a similar fashion, using 6knot splines to define the scale factors. An additional 6-knot spline was added to the fit to provide a smooth background correction for any inconsistencies in the scattered-light subtraction during the extraction process. The scaled, aligned, distorted template spectrum thus had the form gj = αj tj + βj t′j + γj t′′j + δj t′′′ j + ǫj t′′′′ j + ηj (A3) where αj is the value at xj of the 34-knot spline and βj , γj , δj , ǫj and ηj are the values of the corresponding 6-knot 15 splines used to scale the derivatives and correct the background. The knot values for the various splines were determined by the method of least squares using singular-value decomposition to ensure that spurious fluctuations were suppressed in those parts of the spectrum where few lines were present. When the planet velocity is sufficient to shift the absorption lines in the reflected starlight well away from those in the direct starlight, then the residual spectrum r = f − g is expected to retain the reflected starlight signal, Doppler shifted and deeply buried in Poissonian noise, once step 2 in Fig. A1 is complete. APPENDIX B: RESIDUAL FIXED-PATTERN NOISE In this appendix we describe the principal-component analysis method used to correct the spatially fixed but timedependent ripples found in the residual spectra following subtraction of the template. Typically the grouped spectra contained 4×106 or more photons per pixel over most of the recorded spectrum. The scatter of the pixel values in the residual spectrum was compared with expectations from photon statistics, by dividing each grouped spectrum through by the square root of its computed variance. The distribution of the pixel deviates calculated over the 62 orders used for the deconvolution (see Appendix C below) was found typically to contain a few hundred extreme outliers, mostly caused by systematic errors in the polynomial fits in half a dozen regions affected by defects on the CCD. Despite clipping at ±4 times the expected local root-mean-square (rms) photon noise amplitude, we found that the distribution of the remaining values was Gaussian with an rms error between 1.4 and 1.6 times the expected value. We found a significant correlation between pixel values in successive difference frames, indicating that some fixed-pattern noise sources had not been eliminated fully by the template alignment and subtraction. The residuals were correlated over a few hours in time. The residuals in pairs of frames taken several hours apart or on different nights were also found to be correlated, though in some instances the slope of the correlation was reduced or even reversed. The fixed-pattern noise is visible as ripples in the example frame shown following step 2 in Fig. A1. We mapped these spatially fixed but time-varying patterns in the difference spectra using principal-component analysis (PCA). We start with a set of M spectra, each of length N pixels. The ith spectrum thus has elements si = {si1 , si2 , . . . , siN } and associated variances σ 2 i = 2 2 2 {σi1 , σi2 , . . . , σiN }. The spatial correlation matrix of the time variations in the individual spectral bins of this set of spectra is a real symmetric matrix of dimensions N × N , whose elements are given by: Hjk = M X (sij − ŝj )(sik − ŝk ) i=1 σij σik . Note that we subtract the inverse variance weighted mean spectrum ŝ before computing the correlation matrix. The 16 A.Collier Cameron, K. Horne, A. Penny and C. Leigh Figure A1. Schematic illustration of the main steps in the construction and subtraction of the template spectrum, and the subsequent fixed-pattern noise removal. jth spectral bin of ŝ is given by: PM 2 sij /σij ŝj = Pi=1 . M 2 i=1 1/σij The principal components are those eigenvectors of the correlation matrix accounting for the largest fraction of the variance, and in our implementation are computed using Jacobi’s method (Press et al. 1992). To keep computing time down, we processed the spectra in segments of length 100 pixels. We found that most of the unwanted additional variance was contained in the first two principal components, indicating that two independent sources of fixed-pattern noise were present. The first principal component contained a spatially fixed but temporally-varying pattern of ripples affecting all orders, which is probably attributable to a slow drift in the sensitivity pattern of the flat field and/or thermal flexure in the spectograph. The second principal component consisted mainly of imperfectly-subtracted telluric absorption features in some orders. We computed the contributions of these two fixed-pattern noise sources to each exposure and subtracted them from the difference frames. This effectively removed the correlation between successive frames, and reduced the RMS scatter in the pixel values to within ±10% of the value expected from photon statistics alone (step 3 in Fig. A1). The planet signature should not be affected by this procedure any more than it is affected by the template subtraction because, unlike the fixed-pattern noise, it is smeared out by orbital motion. APPENDIX C: LEAST-SQUARES DECONVOLUTION We extracted the reflected-light signature from the difference frames using the method of weighted least-squares deconvolution (LSD) described by Donati et al. (1997). This method, shown schematically in Fig. C1) entails taking a list of lines with relative strengths appropriate to the type of star concerned, and computing via the method of least squares a “mean” line profile which, when convolved with the line pattern, gives an optimal match to the observed spectrum. The deconvolved profile thus incorporates an average broadening function that is representative of all the lines recorded in the spectrum. Applied to the residual spectrum from which the direct stellar signal has been removed, LSD is an effective way of measuring the average line profile of the faint reflected-light signature of the planet, because the reflected light should should have the same pattern of absorption-line positions and strengths as the stellar spectrum. In terms of signal improvement, LSD is analogous to aligning and averaging (with appropriate weighting factors for line strength and local continuum signal strength on the recorded frame) the profiles of all the individual photospheric absorption lines recorded on each echellogram and included in the line list. As each individual spectral line appears in at least two adjacent echelle orders, we have 7700 images of the 3450 spectral lines listed in the observed wavelength range. If all lines were of equal strength and the continuum signal were constant over the whole frame, the signal-to-noise deconvolved profile would be proportional to the square root of the number of line images used, i.e. nearly 70 times greater than the signal-to-noise ratio of a single line in the original spectrum. In practice, the recorded con- A search for starlight reflected from υ And’s innermost planet 17 Figure C1. Schematic illustration of the main steps in the least-squares deconvolution of a time-series of spectra. tinuum is not uniform and the lines have a wide variety of depths. Even so, the signal-to-noise ratio of the deconvolved profile is found to be some 30 times greater than that of a single line in the best-exposed parts of the original spectrum. The least-squares fitting procedure has the additional advantage that neighbouring, blended lines are treated simultaneously, thereby eliminating the sidelobes that would be produced by a simpler shift-and-add procedure. We divide the observed residual spectrum by the continuum fit to the original spectrum, to obtain a normalised residual spectrum r expressed in units of the continuum level of the original spectrum. We modify the variances associated with the elements of r accordingly. We treat r as the convolution of a “mean” line profile z(v) with a set of weighted delta functions at the wavelengths of a comprehensive list of spectral lines. The profile z is defined on a linear velocity scale, with a velocity increment ∆v = 3 km s−1 per bin, which is close to the average velocity increment per pixel in the extracted spectra. The elements Ajk of the convolution matrix A are computed by summing, over all spectral lines ℓ, the fractional contribution of the element of the deconvolved profile at velocity vk to the data pixel at wavelength λj when the centre of the deconvolved profile is shifted to the wavelength λℓ of each line in turn. Hence Ajk = X ℓ wℓ Λ[(vk − c(λj /λℓ − 1))/∆v]. (C1) The triangular function Λ has the form Λ(x) = 1 + x for −1 < x ≤ 0, Λ(x) = 1 − x for 0 < x < 1, and is zero everywhere else. The line weights wℓ , incorporated in the convolution matrix A, are proportional to the central depths of the lines as computed from a Kurucz model atmosphere for the appropriate spectral type. After some experimentation we found that a line list computed for a G2 spectral type with solar abundances gave the best results. The use of a cooler template gives a better fit to the line depths than an earlier spectral type, presumably because of the star’s above-solar metallicity. Because the form of the geometric albedo spectrum of υ And b is unknown, we investigated the effects of both grey (i.e. wavelength-independent) and non-grey geometric albedo spectra. The non-grey models we employed were the Class V and “isolated” Class IV models of Sudarsky, Burrows & Pinto (1999) (Fig. 1). The theoretical albedo spectra provided an additional weighting factor in the deconvolution linelist, allowing us to place limits on the planet’s radius for any assumed albedo model, and to assess the relative goodness-of-fit to the data for different atmospheric models. The deconvolved profile z(v) has the form of a line profile normalized to unit continuum intensity, but from which the continuum has been subtracted. Determination of z(v) via least squares entails minimizing the magnitude of the misfit vector |r − A · z|, χ2 = (r − A · z)T · V · (r − A · z), (C2) weighted so as to make due allowance for the observational errors σj associated with the individual spectral bins. Here the inverse variances σj2 associated with the N elements rj of the residual spectrum r are incorporated via the diagonal matrix 2 V = Diag[1/σ12 , . . . , 1/σN ]. (C3) The least-squares solution for z is found by solving the matrix equation AT · V · A · z = AT · V · r. (C4) 18 A.Collier Cameron, K. Horne, A. Penny and C. Leigh Since the square matrix AT ·V .A is symmetric and positivedefinite, the least-squares problem can be solved using efficient methods such as Cholesky decomposition (Press et al. 1992). The deconvolved profile z is expressed in units of the weighted mean continuum level. The deconvolution procedure compensates for local line blends, and so has the advantage over cross-correlation methods that outside the region occupied by the residual stellar profile, the deconvolved spectrum is flat. The formal errors on the M points of the deconvolved profile z are obtained in the usual way from the diagonal elements of the M × M covariance matrix C = [AT · V · A]−1 . (C5) planets in our own solar system. Jupiter and Venus appear to have phase functions that are more strongly back-scattering than a Lambert sphere. Photometric studies of Jupiter at large phase angles from the Pioneer and subsequent missions have shown (Hovenier & Hage 1989) that the phase function for Jupiter is very similar to that of Venus. As a plausible alternative to the Lambert-sphere formulation, we use a polynomial approximation to the empirically determined phase function for Venus (Hilton 1992). The phase-dependent correction to the planet’s visual magnitude is approximated by: ∆m(α) = 0.09(α/100◦ )+2.39(α/100◦ )2 −0.65(α/100◦ )3 .(D5) so that g(α) = 10−0.4∆m(α) . (D6) APPENDIX D: MATCHED-FILTER ANALYSIS D2 To detect the planet we must co-align and add its signal from profiles at many different orbit phases. Given an assumed orbital inclination, we can compute the sinusoidal path that the planet should follow through velocity space, together with the attendant changes in signal strength. To search for a pattern of faint features displaying the expected behaviour, we construct a set of matched filters that can be scaled to fit the data for each member of a set of appropriate orbital inclinations, as shown schematically in Fig. D1. We model the velocity profile of the reflected-light signal as a sequence of gaussians with appropriate velocities and relative amplitudes: The template spectrum that is subtracted in turn from each spectrum contains a blurred planet signal, as described in Appendix A. Our matched filter must therefore mimic accurately the effects of subtracting a template constructed from the spectra themselves, so we subtract the inverse variance2 weighted average of the travelling gaussian. If σij is the variance of the ith velocity bin in the jth deconvolved profile, then the attenuated basis function becomes: G(v, φ|Kp ) = W⋆ √ g(φ, i) × ∆vp 2π " 1 exp − 2  v − Kp sin φ ∆vp 2 # . (D1) The amplitude Kp of the sinusoidal velocity variation is determined by the system inclination and stellar mass according to Kp ≃ 139  M⋆ 1.3M⊙ 1/3 sin ikm s−1 . (D2) The phase angle α can be determined at any orbital phase from cos α = − sin i. cos 2πφ. (D3) W⋆ is the integrated area of the deconvolved stellar line profile, and ∆vp is the characteristic width of the planet’s reflected-light profile. The width ∆vp = 9.1 km s−1 of the gaussian representing the planetary profile was determined from a fit to the deconvolved profile of υ And. Attenuation during starlight subtraction hij (Kp ) = G(vi , φj |Kp ) − Dij = dij − g(α) = (sin α + (π − α) cos α)/π. P This assumes that the planetary atmosphere scatters isotropically into 2π steradians. Second, it may be more realistic to adopt a phase function that resembles those for the cloud-covered surfaces of P j 2 dij /σij 2 1/σij i Hij = hij − D3 P 2 1/σij . (D7) (D8) 2 hij /σij . 2 1/σij i (D9) Scaling the matched filter The attenuated basis function H(v, φ, Kp ) is normalised such that when it is scaled to give an optimal fit to the entire time-series of Dij values, the scaling factor is ǫ0 , the planet/star flux ratio at α = 0 as defined in eq. 15. The optimal scale factor ǫ0 and its formal error are determined from 2 X Dij Hij (Kp )/σij i,j (D4) 2 G(vi , φj |Kp )/σij gives the resulting orthogonalised data pixel value. The matched filter is orthogonalised in the same way: Phase function For the phase function there are two natural choices. First, the phase function of a Lambert sphere is j To allow for any systematic errors in the continuum level following template subtraction and deconvolution, we subtract the inverse variance-weighted mean value of each of the deconvolved profiles at each orbital phase. If dij is the original data value at velocity i and phase j, then ǫ0 (Kp ) = D1 P and Var(ǫ0 ) = X i,j 2 2 Hij (vi , φj , Kp )/σij 1 . 2 2 Hij (Kp )/σij (D10) (D11) We exclude all pixels within 25 km s−1 of the stellar line core from the summations, to eliminate spurious effects arising from the ripple pattern in the core of the residual deconvolved stellar profile. A search for starlight reflected from υ And’s innermost planet 19 Figure D1. Schematic illustration of the main steps in the matched-filter analysis. D4 Calibrating non-grey albedo models. The meaning of ǫ0 in the expressions above is obvious if the albedo spectrum is grey, i.e. if ǫ0 is independent of wavelength. Problems arise, however, when we wish to test how well a given non-grey albedo spectrum fits the data. In this case, the weighting factors wi used for the lines in the deconvolution mask mask are multiplied by the albedo as described in Appendix C above, and ǫ0 follows the wavelength dependence of the geometric albedo. We circumvented this difficulty by replacing ǫ0 with (Rp /a)2 as the scaling factor in eq. D10. The basis function G is rescaled to become G′ = hpi G, where hpi is an appropriately weighted wavelength-averaged geometric albedo, then attenuated as before using eq. D7. For each candidate albedo spectrum we determine the appropriate value of hpi empirically. We inject an artificial planet signal with the required albedo spectrum and known Kp and Rp /a into the data, and deconvolve the synthetic data with the albedo-weighted line list. We subtract the profiles of the observed spectra, following deconvolution with the same line list, to isolate the deconvolved planet signal from the noise. We then back-project the isolated planet signal, and measure the value of hpi that recovers the correct value of Rp /a at the appropriate Kp . APPENDIX E: FALSE-ALARM PROBABILITIES FOR CANDIDATE SIGNALS Candidate reflected-light signals are characterised by positive values of (Rp /a)2 and an improvement in χ2 relative to the fit obtained when (Rp /a)2 = 0. We determine the frequency with which (Rp /a)2 exceeds a given value due to noise in the absence of a planet signal, using a bootstrap procedure. In each of 3000 trials, we randomize the order in which the three nights of observations were secured, then we randomise the order in which the observations were secured within each night. The re-ordered observations are then associated with the original sequence of dates and times. This ensures that any contiguous blocks of spectra containing similar systematic errors remain together, but appear at a new phase. Any genuine planet signal present in the data is, however, completely scrambled in phase. The re-ordered data are therefore as capable as the original data of producing spurious detections through chance alignments of blocks of systematic errors along a single sinusoidal path through the data. We record the least-squares estimates of (Rp /a)2 and the associated values of χ2 as functions of Kp in each trial. The 3000 trials implicitly define empirical probability distributions of (Rp /a)2 and χ2 that include both the photon statistics and the effects of correlated systematic errors at each trial value of Kp . The percentile points of the distribution of (Rp /a)2 values at each Kp are used to define the 1σ, 2σ, 3σ and 4σ upper-limit contours shown in Figs. 10, 11, 10, 14 and 17. To assess the false-alarm probability for a candidate detection, however, we need to examine the likelihood of the fit to the data. In each bootstrap trial, the most likely candidate is the one among those with (Rp /a)2 > 0 that yields the greatest improvement ∆χ2 relative to the no-planet hypothesis. Such a peak can occur at any value of Kp . We use the distribution of the peak ∆χ2 value in each trial to determine how often we would expect the ∆χ2 of the strongest spurious noise feature giving (Rp /a)2 > 0 to exceed the actual ∆χ2 of a candidate detection in the observed data. We define this probability – summed over all possible values of Kp with some weighting according to prior probability estimates – to be the “false-alarm probability” for a candidate signal. The unweighted false-alarm probabilities are based on the likelihood of obtaining the data D at a given Kp after 20 A.Collier Cameron, K. Horne, A. Penny and C. Leigh optimizing (Rp /a)2 , and take no account of whether or not that value of Kp is plausible. This conditional likelihood is measured relative to the likelihood of getting the same data if no planet signal is present: 2 2 P (D|Kp ) = e−(χ (Kp )−χ (No planet))/2 , P (D|No planet) (E1) where we have implicitly optimised the value of (Rp /a)2 using the matched-filter fit at each Kp . We take our prior knowledge of Kp into account by modifying the likelihood to give the joint probability P (D, Kp ) = P (D|Kp )P (Kp ). (E2) We use the prior probability distribution for Kp as shown projected on to the Kp axis in Fig.2. We compute P (D, Kp ) for each of the 3000 bootstrap trials. To do this, we pick out the greatest value in each trial of the quantity ln P (D, Kp ) = ln P (Kp ) + ∆χ2 2 (E3) where P (Kp ) is derived from the simulations described in Section 2. The cumulative distribution of maximal log P (D, Kp ) values derived from the bootstrap trials gives the probability that the observed peak value of ln P (D, Kp ) could be spurious.