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

Abstract

In data from three clear nights of échelle spectroscopy in 2000 October/November, 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 R_p as functions of its projected orbital velocity K_p≈139 sin i km s⁻¹ for various assumptions about the wavelength-dependent geometric albedo spectrum p(λ) of its atmosphere. For a grey albedo p we find , with 0.1 per cent false-alarm probability (4σ). For a Sudarsky, Burrows & Pinto Class V model atmosphere, the mean albedo in our 380–676 nm bandpass is 〈p〉∼0.42, requiring R_p<1.51R_J, whereas an (isolated) Class IV model with 〈p〉∼0.19 requires R_p<2.23R_J. The star's v_rot sin i∼10 km s⁻¹ and estimated rotation period P_rot∼10 d 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.

methods: data analysis, techniques: spectroscopic, stars: individual: ν And, planetary systems

Introduction

The discovery of ‘hot Jupiters’– giant exoplanets in 3- to 5-d 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 as 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 R_p∼1.0–1.7R_Jup for 1<M_p/M_Jup<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 mass–radius–age relation for hot Jupiters is determined largely by (i) the age at which they arrived in their present close orbits, and (ii) 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 209 458 (Charbonneau et al. 2000) provided conclusive proof of the existence and gas-giant nature of the ‘hot Jupiters’. Ultra-precise photometry of HD 209 458 b's transit profile using the Space Telescope Imaging Spectrograph (STIS) aboard the Hubble Space Telescope (HST) yields a radius about 1.35 R_Jup (Brown et al. 2001), although the planet has a mass only 0.63 M_Jup. 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 (T_eff∼1400 K) and gravities (g≥36 m s⁻²) of massive planets such as τ Boo b, possible cloud condensates include upper cloud decks of MgSiO₃ (enstatite) and Al₂O₃ (corundum) at pressures ∼5 to 10 bar, and a deeper layer of Fe grains at pressures of 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. Sudarsky et al. (1999) predict this ‘Class IV’ albedo spectrum for planets with high surface gravities. At the lower surface gravities, g≃10 m s⁻², expected of lower-mass planets like υ And b and HD 209 458 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 per cent that of Jupiter throughout most of the optical spectrum (Fig. 1).

Figure 1.

Geometric albedo spectra for the Class V (solid line), isolated (dashed) and irradiated (dot-dashed) Class IV models of Sudarsky et al. (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 10¹¹photons μm⁻¹, and gives the total photon deficit for a typical group of four exposures.

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 ε₀≃10⁻⁴. Collier Cameron et al. (1999) claimed a detection with an average ε₀=7.5×10⁻⁵ from data secured with the William Herschel Telescope (WHT) 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, ε₀<3.5×10⁻⁵, between 387.4 and 586.3 nm. The apparent tidal synchronization of τ Boo with the planet's orbit suggests a low inclination of the order of 40°. At this inclination, the combined 1998–2000 WHT observations yielded a 99.9 per cent upper limit on the geometric albedo p≤0.22, assuming a planet radius of 1.2 R_Jup, 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.

System Parameters

υ And (HR 458, HD 9826) is a late-F main-sequence star with parameters as listed in Table 1. High-precision radial-velocity studies have shown it to have a system of three 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 summarized in Table 1.

Table 1.

Physical parameters for u And and its innermost planet.

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)(R_p/a)²∼10⁻⁴ 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 ε≡f_p/f_⋆ as a function of orbital phase φ and wavelength λ. The information we aim to obtain (in order of increasing difficulty) comprises:

• (i)

  K_p, the planet's projected orbital velocity. From this we learn the orbital inclination i and hence the planet mass M_p, as M_p sin i is known from the Doppler wobble of star.

• (ii)

  ε₀, the maximum strength of the reflected starlight. From this we constrain the radius of the planet, as ε₀=p(R_p/a)², where p is the geometric albedo;

• (iii)

  p(λ), the albedo spectrum, which depends on the composition and structure of the planetary atmosphere;

• (iv)

  g(α,λ), the phase function describing the dependence of the amount of light reflected toward the observer on the star-planet-observer angle α;

• (v)

  Δν_p, the velocity width of lines in the reflected starlight. This depends mainly on the rotation of the star in the frame of the orbiting planet, and to a lesser extent on the rotation and surface winds of the planet.

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 K_p and ε₀ and hence R_p for an assumed albedo p. With sufficient data we can measure ε(λ) 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 R_p/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  

(1)

The orbital phase at time t is given by φ=(t−T₀)/P_orb, using the values of P_orb and T₀ given in Table 1.

The mass ratio q=M_p/M_⋆ is determined from the observed amplitude K_⋆ of the star's reflex orbital velocity via  

(2)

The apparent radial velocity amplitude K_p of the planet about the centre of mass of the system is  

(3)

Kepler's third law gives  

(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 orbital velocity amplitude of the planet is considerably greater than the typical widths of the photospheric absorption lines of the star. Lines in the reflected-light spectrum of the planet should therefore be Doppler-shifted well clear of their stellar counterparts, allowing clean spectral separation for most of the orbit.

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  

(5)

Low inclinations are also ruled out by the chromospheric activity levels of the star. The stellar radius and v sin i suggest an axial rotation period of 9 d. A 12±2 d 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. (1997). The rms scatter in the 30-d means is 4.7 per cent of the long-term mean flux level. The 20 per cent uncertainty in the 11.7-d 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. 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 d, 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  

(6)

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 v sin i=9.5 km s⁻¹ and estimate of its rotation period P_rot∼12±2 d, which together favour a relatively high i. The lower histogram in Fig. 2 shows the resulting probability distribution for K_p (equivalent to i).

Figure 2.

The grey-scale shows the prior joint probability density function (PDF) for projected orbital velocity K_p and the squared ratio (R_p/a)² of the planet radius to the orbit radius, based on the measured system parameters. Darker shades in the grey-scale denote greater probabilities of the planet having the corresponding combination of (R_p/a)² and K_p. 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 K_p≃135 km s⁻¹ and (R_p/a)²≃10⁻⁴. The horizontal line shows the value of (R_p/a)² for a 1−R_Jup planet in the same orbit.

Rotational broadening

The matched-filter analysis used to search for the reflected-light signal requires an a priori estimate of the rotational broadening of the spectral-line profiles in starlight reflected from the atmosphere of the planet.

The Earth-bound observer sees stellar absorption lines that are rotationally broadened by  

(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  

(8)

The Monte Carlo simulation gives the resulting relationship between the broadening v_refl of the reflected starlight and K_p, 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 km s⁻¹ and median 8.4 km s⁻¹. The planet linewidths are reduced at lower orbital inclinations, because the difference between the orbital and stellar rotation frequencies is reduced.

Figure 3.

A Monte Carlo simulation shows how the apparent stellar equatorial rotation speed v_e,refl in the starlight incident on the planet is expected to vary with the projected orbital velocity amplitude K_p. The system parameters are defined as in Section 2. The median rotational broadening of the reflected starlight is 8.4 km s⁻¹, somewhat less than the v sin i=9.5 km s⁻¹ of the direct starlight.

Albedo spectrum

The quantity that we observe is the ratio f_p/f_⋆ of the spectral-line strengths in the planet's light to those in the direct stellar spectrum as a function of wavelength:  

(9)

This depends on the fraction (R_p/2a)² 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  

(10)

The geometric albedo p(λ) is defined at α=0 as  

(11)

Here F_reflected(0,λ)=πI_reflected(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 normalized to g(0,λ)=1.

The flux received from the planet at Earth, is therefore  

(12)

As the stellar flux received at Earth, a distance D away, is  

(13)

equation (9) follows.

Phase function

The phase angle α depends on the orbital inclination i and the orbital phase φ (measured from transit or inferior conjunction):  

(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:  

(15)

Because a is tightly constrained through Kepler's third law, given the period P and the star mass M_⋆, the measurements of ε₀(λ) measure the product .

In this paper we adopt an empirical phase function similar to that of Venus and Jupiter. The expressions given in Appendix D1 (equations D4 and equation D6 for the phase function of a Lambert sphere and Venus, respectively) yield the phase-dependent flux correction factors plotted in Fig. 4.

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.

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 the albedo, we plot the relationship between g and K_p, derived from the Monte Carlo simulations, in Fig. 5. The mass of the planet 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.

Figure 5.

A Monte Carlo simulation shows how the gravitational acceleration g at the surface of the planet is expected to vary with the projected orbital velocity amplitude K_p. The system parameters are defined as in Section 2. For comparison, g_Jup=26.5 m s⁻².

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 mass of the planet. Their radiative–convective gas-giant models predict a radius between 1.2 and 1.7 R_Jup for M_p=0.6M_Jup, to between 1.15 and 1.3 R_Jup for M_p=1.15M_Jup.

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 m s⁻², is close to the ∼10 m s⁻² limit below which the Class V albedo models of Sudarsky et al. (1999) are applicable.

If Class V models apply, we can use them to estimate the brightness of the planet. Viewed fully illuminated, the reflected light from a planet with (R_p/a)²≃10⁻⁴ and Class V albedo p≃0.42 (cf. the Class V models of Sudarsky et al. 1999) is expected to be 24 000 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).

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 October 10 and November 6 and 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 charge-coupled devices (CCDs) manufactured by EEV, each with 2048×4096 13.5-μm pixels. The CCD was centred at 459.6 nm in order 124 of the 31 g mm⁻¹ echelle grating, giving complete wavelength coverage from 381.3 to 675.8 nm with minimal vignetting. The average pixel spacing was close to 1.6 km s⁻¹, and the full width at half maximum intensity of the thorium-argon arc calibration spectra was 3.5 pixels, giving an effective resolving power of R=53 000.

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 s, depending upon seeing, in order to expose the CCD to a peak count of 40 000 ADU pixel⁻¹ in the brightest parts of the image. A 450-s exposure yielded about 1.2×10⁶electrons 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 signal-to-noise ratio 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 signal-to-noise ratio in the continuum of the brightest orders is typically 1000 pixel⁻¹.

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.

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 observations of the entire year.

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-plus-sky 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.

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 24 000 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 orbital motion of the planet, so that to first order the spectrum of the planet is eliminated 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 sensitivity of the detector, thermal flexure in the spectrograph which causes 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 Data Base (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 grey-scale form as a two-dimensional ‘velocity–phase’ map in Fig. 6. This representation of the data shows a strong pattern of distortions in the residual stellar profiles within ±20 km s⁻¹ 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 10⁴ of the average continuum level. They vary too rapidly during the night to be attributable to, for example, 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 October 9 November 6 and 7. The lineweights in the deconvolution were defined assuming a grey albedo spectrum. The grey-scale runs from black at −10⁻⁴ times the mean stellar continuum level, to white at +10⁻⁴. 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 (R_p/a)²=2.57×10⁻⁴. 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-to-star flux ratio should thus be ε₀=f_p/f_⋆=1.08×10⁻⁴. We assumed an orbital inclination of 80°, giving an orbital velocity amplitude K_p=137 km s⁻¹ and a rotational broadening of the reflected starlight, v_e,refl=8±1 km s⁻¹. 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.

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.

Extracting The Planet Signature

The form of the expected planet signature in the velocity–phase diagram can be represented accurately as a travelling Gaussian of characteristic width Δv_p whose velocity varies sinusoidally around the orbit and whose strength is modulated according to the phase function g(α,λ) (see Appendix D1). The value of Δv_p is adjusted to match the expected rotational broadening of the reflected starlight. The amplitude K_p 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: because the planet 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 data set, few observations were made near quadrature so this problem does not arise.

Figure 8.

As for Fig. 6, but showing the matched filter function used to measure (R_p/a)² at K_p=137 km s⁻¹. 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 K_p.

The travelling Gaussian in this image has none the less 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⁻¹ and v≃+100 km s⁻¹ in Fig. 8.

7 Probability Maps for K_p and R_p/a

We use a sequence of such models for different values of K_p 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 K_p we construct a model H_ij(K_p), like that shown in Fig. 8, and scale it to give an optimal fit to the residual velocity–phase map D_ij using the methods presented in Appendix D. For an assumed albedo spectrum p(λ) and orbital velocity amplitude K_p, the fit of the matched filter to the data measures (R_p/a)².

The badness-of-fit statistic χ² gives a relative measure of the probability that any signal detected could have arisen by chance, and is quantified by  

(16)

Here is the variance associated with D_ij, 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 χ² measured relative to the no-planet model. For a noise-dominated signal, the estimated value of (R_p/a)² 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 χ² and is therefore the most probable, but only by a narrow margin.

Figure 9.

The upper panel shows the optimal scaling factor (R_p/a)² plotted against orbital velocity amplitude K_p, assuming a grey albedo spectrum with p=0.42. The lower panel shows the associated reduction in χ² measured relative to the fit obtained in the absence of any planet signal, i.e. for (R_p/a)²=0. Note that only positive values of (R_p/a)² are physically plausible, so the middle peak in the Δχ² plot is unambiguously a noise feature.

The relative probabilities of the fits to the data for different values of the free parameters R_p/a and K_p, given by  

(17)

are shown in Fig. 10. To make this figure, a planet signal pattern as in Fig. 8 was computed for many different K_p, and each one was scaled by (R_p/a)² to fit the residuals in Fig 6. The grey-scale 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.

Figure 10.

Relative probability map of model parameters K_p and log(ε₀/p)=log(R_p/a)², derived from the WHT/UES observations of υ And, assuming a grey albedo spectrum with p=0.42. The grey-scale denotes the probability relative to the best-fitting 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 on the signal strength derived from the bootstrap trials The uppermost contour, for example, represents the value of log(R_p/a)² that was only exceeded in three of 3000 bootstrap trials at each K_p. 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.

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 χ² by chance alone. In principle this could be done using the χ² 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 retaining 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 (R_p/a)² at each K_p are shown as contours in Figs 10 and 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 per cent false alarm probability respectively, for each fixed value of K_p.

Figure 11.

Relative probability map of model parameters K_p and log(ε₀/p)=log(R_p/a)² for a simulated planet signature with grey albedo p=0.42, K_p=137 km s⁻¹, and R_p=2R_Jup. The grey-scale 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.

Grey albedo model

At the most probable values around K_p≃135 km s⁻¹, the grey albedo model yields a 0.1 per cent bootstrap upper limit on the planet/star flux ratio ε₀<5.84×10⁻⁵. Adopting a grey albedo model with unit geometric albedo, we find that the 0.1, 1.0 and 4.6 per cent upper limits on ε₀ 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.

Table 3.

Upper limits on planet dimensions for various albedo models. The upper limits are quoted for an assumed K_p=135 kms⁻¹, near the peak of the prior probability distribution for K_p. 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 √

Two possible planet signals of comparable likelihood are seen. Their properties are listed in Table 4. The stronger, at K_p≃132 km s⁻¹(i≃80°), yields an improvement Δχ²=13.61 over the model fit obtained assuming no planet signal is present. If this feature represents a genuine planet signal, its velocity amplitude K_p=132 km s⁻¹ implies a planet mass M_p=0.74M_Jup and (for p=0.42) a planet radius R_p=1.24±0.17R_Jup. The weaker candidate has K_p≃54 km s⁻¹(i=22°) and Δχ²=11.86, and gives a larger planet radius.

Table 4.

Orbital velocity amplitude, signal strength, planet radius, Dx 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 s21. 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 ? 132km s21 candidate is found to be between 9 and 10 percent.

We used the bootstrap simulations to determine the probability that a spurious detection with Δχ²>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 per cent. These contours give the false-alarm probability only if the value of K_p 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 K_p can exceed the Δχ² of the candidate. If we assume that all values of K_p are equally likely in the range 44<K_p<152 km s⁻¹ over which our phase coverage allows signals to be detected, the false-alarm probability is found to be 26 per cent 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 K_p 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 K_p using our knowledge of the a priori probability distribution for K_p, 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 K_p=132 and K_p=54 km s⁻¹ are comparable, the latter is almost certainly spurious. The false-alarm probability for the K_p=132 km s⁻¹ feature drops to 9.4 per cent when prior knowledge of K_p is accounted for, while that for the 54 km s⁻¹ 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(R_p/a)²=2.57×10⁻⁴ and orbital velocity amplitude velocity K_p=137 km s⁻¹. This most probable combination of orbital velocity and planet radius yields an improvement Δχ²=98.6 with respect to the value obtained assuming no planet is present (i.e. R_p=0). This is far greater than the greatest Δχ²=55.3 produced at any value of K_p 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 (R_p/a)² 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.

Class V roaster model

The ‘Class V roaster’ is the most highly reflective of the models computed by Sudarsky et al. (1999). This model is characteristic of planets with T_eff≥1500 K and/or surface gravities lower than ∼10 m s⁻², 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 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 R_p=2R_jup. For the observations we used Δv_p=8 km s⁻¹ which again yielded the best fit, with scale factors and improvements in χ² as plotted in Fig. 13. The probability map for the observed signals is shown in Fig. 14.

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.

Figure 13.

As for Fig. 9, showing the optimal scaling factor (R_p/a)² and the associated improvement in χ² plotted against orbital velocity amplitude K_p, assuming a Class V albedo spectrum.

Figure 14.

Relative probability map of model parameters K_p and log(ε₀/p)=log(R_p/a)², derived from the WHT/UES observations of υ And, assuming the albedo spectrum to be that of a Class V roaster. The grey-scale and contours are defined as in Fig. 10.

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 K_p=52 and 132 km s⁻¹ are also present with this albedo model. As in the grey case, the feature at K_p=132 is the most probable, this time by a slightly wider margin. Their signal strengths and false-alarm probabilities are given in Table 4.

Isolated Class IV model

The Class IV roaster models of Sudarsky et al. (1999) have a more deeply-buried cloud deck than the Class V models. The resonance lines of Na i and K i are strongly saturated, with broad wings due to collisions with H₂ 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 temperature-pressure 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.

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.

Figure 16.

As for Fig. 9, showing the optimal scaling factor (R_p/a)² and the associated improvement in x² plotted against orbital velocity amplitude K^p, assuming an isolated Class IV albedo spectrum.

The plots of (R_p/a)² and Δχ² versus K_p (Fig. 16) 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) none the less shows the K_p=132 km s⁻¹ probability maximum to be present and again stronger than the 49 km s⁻¹ feature. The fit to the data is, however, poorer than in the grey and Class V cases, giving an improvement of only Δχ²=10.8 over the no-planet hypothesis and a planetary radius R_p=1.68±0.25R_Jup (Table 4). This is 25 per cent larger than the radius estimated with the Class V albedo model. The bootstrap test gives a false-alarm probability of 10 per cent.

Figure 17.

Relative probability map of model parameters K_p and log(ε₀/p)=log(R_p/a)², derived from the WHT/UES observations of υ And, assuming the albedo spectrum to be that of an ‘isolated’ Class IV roaster. The grey-scale and contours are defined as in Fig. 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 Class V model using the Class IV line weights and flux calibration. The radius of the planet was over-estimated by 27 per cent 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 K_p=132 km s⁻¹ 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.

Plausibility of candidate signals

Comparison of Figs 10, 14 and 17 to the joint prior probability distribution in Fig. 2 shows that the feature near K_p=50 km s⁻¹ yields a value of (R_p/a)² that is considerably greater than the radius predicted by Guillot et al. (1997). The feature at K_p=132±3 km s⁻¹, however, has ε₀=4.6×10⁻⁵. A plausible geometric albedo p=0.42 then gives (R_p/a)²=(1.09±0.30)×10⁻⁴, 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 using values of the Gaussian width parameter in the range 6.4<Δv_p<10 km s⁻¹, corresponding to 0<v_e,refl<11.3 km s⁻¹. The feature at 132 km s⁻¹ yields the greatest improvement in χ² when Δv_p=8±1 km s⁻¹, or v_e,refl=7.0±2 km s⁻¹. The 52 km s⁻¹ feature, however, grows less significant as v_e,refl is decreased to the (retrograde) 5 km s⁻¹ or so expected at the corresponding orbital inclination. The rotational broadening of the 132 km s⁻¹ 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 so as to ensure that the period and epoch were uncorrelated. The local minimum in χ² was found to be centred at K_p=133±3 km s⁻¹, P=4.621±0.005 d, and T₀=2 451 853.770±0.025. The radial-velocity ephemeris gives P=4.617 07±0.00003 d and T₀=2 451 853.777±0.012, well within the 1σ error bars derived from the back-projection. The 52 km s⁻¹ feature, on the other hand, is part of a more extended structure whose peak is located at K_p=67±3 km s⁻¹, P=4.651±0.006 d, and T₀=2 451 853.81±0.02. This bears little relation to the radial-velocity solution, suggesting a spurious origin.

We conclude that the candidate feature at K_p=132 km s⁻¹ corresponds to a local minimum of χ² 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 cannot easily dismiss this feature as being merely an extension of a larger, spurious noise-induced structure.

Conclusions

The observations of υ And presented here rule out radii for the innermost planet greater than with 0.1 per cent false-alarm probability if a grey albedo spectrum is assumed. We derive upper limits also with 0.1 per cent false-alarm probability on the planet's radius of 1.53 R_jup, assuming the albedo spectrum of a Class V roaster (low gravity, high cloud deck; Sudarsky et al. 1999), or 2.23 R_jup 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 et al. 1999) if its radius is comparable to that of HD 209 458 b. The evidence for a reflected-light signature in the observations reported here, at a projected orbital velocity amplitude K_p=132±3 km s⁻¹, 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 R_p=1.34±0.17R_Jup if a Class V spectrum is assumed. This is comparable to the radius of HD 209 458 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 per cent 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 y 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 of PPARC. We thank the referee for raising several useful points of clarification.

Appendix

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 five 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σ 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).

Figure 18.

Schematic illustration of the main steps in the construction and subtraction of the template spectrum, and the subsequent fixed-pattern noise removal.

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 t_j in each order as a function of pixel number x_j in the dispersion direction  

(A1)

and  

(A2)

The third and fourth derivatives were computed by applying the same operations to t″, for use in correcting frame-to-frame changes in the higher moments of the PSF and seeing profile.

The vector ƒ 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 least-squares 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 6-knot 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  

(A3)

where α_j is the value at x_j of the 34-knot spline and β_j, γ_j, δ_j, ε_j and η_j are the values of the corresponding 6-knot 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=ƒ−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 time-dependent ripples found in the residual spectra following subtraction of the template.

Typically the grouped spectra contained 4×10⁶ 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 s_i={s_i1,s_i2,…,s_iN} and associated variances . 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:  

Note that we subtract the inverse variance weighted mean spectrum s^ before computing the correlation matrix. The jth spectral bin of s^ is given by:  

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 spectrograph. 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 per cent 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 have the same pattern of absorption-line positions and strengths as the stellar spectrum.

Figure 19.

Schematic illustration of the main steps in the least-squares deconvolution of a time-series of spectra.

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 continuum 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 normalized 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⁻¹ per bin, which is close to the average velocity increment per pixel in the extracted spectra. The elements A_jk 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 v_k 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  

(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 lineweights 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 linelist 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 above-solar metallicity of the star.

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 et al. (1999) (Fig. 1). The theoretical albedo spectra provided an additional weighting factor in the deconvolution linelist, allowing us to place limits on the radius of the planet 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∼,  

(C2)

weighted so as to make due allowance for the observational errors σ_j associated with the individual spectral bins. Here the inverse variances associated with the N elements r_j of the residual spectrum r are incorporated via the diagonal matrix  

(C3)

The least-squares solution for z is found by solving the matrix equation  

(C4)

As the square matrix A^T·V·A is symmetric and positive-definite, 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  

(C5)

Appendix D: Matched-Filter Analysis

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.

Figure 20.

Schematic illustration of the main steps in the matched-filter analysis.

We model the velocity profile of the reflected-light signal as a sequence of Gaussians with appropriate velocities and relative amplitudes:  

(D1)

The amplitude K_p of the sinusoidal velocity variation is determined by the system inclination and stellar mass according to  

(D2)

The phase angle α can be determined at any orbital phase from  

(D3)

W_⋆ is the integrated area of the deconvolved stellar line profile, and Δv_p is the characteristic width of the reflected-light profile of the planet. The width Δv_p=9.1 km s⁻¹ of the Gaussian representing the planetary profile was determined from a fit to the deconvolved profile of υ And.

D1 Phase function

For the phase function there are two natural choices. First, the phase function of a Lambert sphere is  

(D4)

This assumes that the planetary atmosphere scatters isotropically into 2π steradians.

Secondly, it may be more realistic to adopt a phase function that resembles those for the cloud-covered surfaces of 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 visual magnitude of the planet is approximated by:  

(D5)

so that  

(D6)

D2 Attenuation during starlight subtraction

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 variance-weighted average of the travelling Gaussian. If is the variance of the ith velocity bin in the jth deconvolved profile, then the attenuated basis function becomes:  

(D7)

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 d_ij is the original data value at velocity i and phase j, then  

(D8)

gives the resulting orthogonalized data pixel value. The matched filter is orthogonalized in the same way:  

(D9)

D3 Scaling the matched filter

The attenuated basis function H(v,φ,K_p) is normalized such that when it is scaled to give an optimal fit to the entire time-series of D_ij values, the scaling factor is ε₀ and the planet/star flux ratio at α=0 is as defined in equation (15).

The optimal scalefactor ε₀ and its formal error are determined from  

(D10)

and  

(D11)

We exclude all pixels within 25 km s⁻¹ 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.

D4 Calibrating non-grey albedo models

The meaning of ε₀ in the expressions above is obvious if the albedo spectrum is grey, i.e. if ε₀ 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 w_i used for the lines in the deconvolution mask are multiplied by the albedo as described in Appendix C above, and ε₀ follows the wavelength dependence of the geometric albedo.

We circumvented this difficulty by replacing ε₀ with (R_p/a)² as the scaling factor in equation (D10).

The basis function G is rescaled to become G′=〈p〉G, where 〈p〉 is an appropriately weighted wavelength-averaged geometric albedo, then attenuated as before using equation (D7). For each candidate albedo spectrum we determine the appropriate value of 〈p〉 empirically. We inject an artificial planet signal with the required albedo spectrum and known K_p and R_p/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 linelist, to isolate the deconvolved planet signal from the noise. We then back-project the isolated planet signal, and measure the value of 〈p〉 that recovers the correct value of R_p/a at the appropriate K_p.

Appendix E: False-Alarm Probabilities for Candidate Signals

Candidate reflected-light signals are characterized by positive values of (R_p/a)² and an improvement in χ² relative to the fit obtained when (R_p/a)²=0. We determine the frequency with which (R_p/a)² exceeds a given value as a result of noise in the absence of a planet signal, using a bootstrap procedure.

In each of the 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 (R_p/a)² and the associated values of χ² as functions of K_p in each trial. The 3000 trials implicitly define empirical probability distributions of (R_p/a)² and χ² that include both the photon statistics and the effects of correlated systematic errors at each trial value of K_p.

The percentile points of the distribution of (R_p/a)² values at each K_p are used to define the 1σ, 2σ, 3σ and 4σ upper-limit contours shown in Figs 10, 11, 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 (R_p/a)²>0 that yields the greatest improvement Δχ² relative to the no-planet hypothesis. Such a peak can occur at any value of K_p. We use the distribution of the peak Δχ² value in each trial to determine how often we would expect the Δχ² of the strongest spurious noise feature giving (R_p/a)²>0 to exceed the actual Δχ² of a candidate detection in the observed data. We define this probability — summed over all possible values of K_p 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 K_p after optimizing (R_p/a)², and take no account of whether or not that value of K_p is plausible. This conditional likelihood is measured relative to the likelihood of getting the same data if no planet signal is present:  

(E1)

where we have implicitly optimized the value of (R_p/a)² using the matched-filter fit at each K_p.

We take our prior knowledge of K_p into account by modifying the likelihood to give the joint probability  

(E2)

We use the prior probability distribution for K_p as shown projected on to the K_p axis in Fig.2.

We compute P(D,K_p) for each of the 3000 bootstrap trials. To do this, we pick out the greatest value in each trial of the quantity  

(E3)

where P(K_p) is derived from the simulations described in Section 2. The cumulative distribution of maximal log P(D,K_p) values derived from the bootstrap trials gives the probability that the observed peak value of ln P(D,K_p) could be spurious.