A&A 431, 1157-1166 (2005)
DOI: 10.1051/0004-6361:20041640

On measuring planetary winds using high-resolution spectroscopy in visible wavelengths

T. Civeit ¹ - T. Appourchaux ^2,3 - J.-P. Lebreton ¹ - D. Luz ^4,5 - R. Courtin ⁴ - C. Neiner ^1,6 - O. Witasse ¹ - D. Gautier ⁴

1 - Research and Scientific Support Department, ESA/ESTEC, PO Box 299, 2200 AG Noordwijk, The Netherlands
2 - Payloads and Advanced Concepts Office, ESA/ESTEC, PO Box 299, 2200 AG Noordwijk, The Netherlands
3 - Institut d'Astrophysique Spatiale, Université Paris-Sud, Bât. 121, 91405 Orsay Cedex, France
4 - Observatoire de Paris, LESIA, 5 place Jules Janssen, 92195 Meudon Cedex, France
5 - Observatório Astronómico de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal
6 - GEPI, UMR 8111 du CNRS, Observatoire de Paris-Meudon, 92195 Meudon Cedex, France

Received 12 July 2004 / Accepted 11 October 2004

Abstract
We present a new method that uses high-resolution spectroscopy in the visible wavelength domain to measure planetary winds. A rotating atmosphere illuminated by the Sun induces a Doppler shift in the back-scattered solar light. Its analysis with a cross-dispersed echelle spectrometer allows the direct determination of both the wind speed and direction. We describe in this paper the image and data processing algorithms used in the method and implemented in a data reduction and analysis package. Since the velocity of planetary winds leads to Doppler shifts smaller than the width of the solar lines, accurate Doppler measurements are performed by running the algorithm proposed by Connes (1985), which is an optimum technique using the full available spectral information. We apply the method to Io as a test case, a small slowly-rotating body with no atmosphere, to measure its solid rotation. The observations span wavelengths from 414 to 621 nm and were carried out with the Ultraviolet and Visual Echelle Spectrograph (UVES) on the 8.2 m Kueyen unit at the Very Large Telescope (VLT - ESO). The results we obtain for Io validate the principle of the method by optimally measuring the well-known surface rotational velocity of this moon, with an uncertainty smaller than 2 m s^-1. However, the analysis of the set of observations shows that systematic errors are large and one needs to consider the retrieved velocity as a lower limit.

Key words: instrumentation: spectrographs - methods: data analysis - techniques: spectroscopic - planets and satellites: individual: Io - planets and satellites: individual: Titan

1 Introduction

The aim of this study is to improve our knowledge of the wind direction and velocity in the stratosphere of Titan. Titan, Saturn's largest moon, will be the target of the Huygens probe (Lebreton & Matson 2002) on January 14, 2005. One of the numerous scientific objectives of this mission is to measure the zonal wind, from an altitude of about 160 km down to the surface. Up to now, few measurements of Titan's winds have been made. The stratospheric winds are generally considered to be zonal, with a velocity in the range 100-200 m s^-1 (Flasar et al. 1981; Kostiuk et al. 2001; Hubbard et al. 1993). The knowledge of winds is very important for the understanding of the atmospheric composition, dynamics and energetics. For the Huygens mission and from an engineering point of view, winds are also important, especially for the probe performance and for the reconstruction of the trajectory.

With a high-resolution spectrometer, we analyse in the visible wavelength domain the solar light back-scattered by a planetary atmosphere. A rotating atmosphere induces a Doppler shift in the reflected solar spectrum, whose analysis allows the direct determination of both the wind speed and direction. One of the techniques involved in deriving the wind velocity is based on Doppler measurements which use the Absolute Accelerometry proposed by Connes (1985). The Connes technique has been succesfully applied to detect the first evidence of oscillation modes in solar-type stars other than the Sun (Martic et al. 1999). Moreover, its applicability to planetary winds has been demonstrated on Venus (Martic et al. 2001) by using the ELODIE fibre-fed echelle spectrograph installed at the Observatoire de Haute-Provence (OHP); a mean retrograde flow of the order of 100 m s^-1 was measured.

Observations of Titan were performed in 2002 and 2003 with the Ultraviolet and Visual Echelle Spectrograph (UVES) at the Very Large Telescope (VLT - ESO). However, the data reduction and analysis using the UVES pipeline led to results difficult to interpret. Therefore, it was decided to implement new procedures and to test them on a rotating body with no atmosphere, namely Jupiter's moon Io, whose apparent diameter and surface rotational velocity are similar to the quantities involved in the case of Titan. Since the solar light is reflected by Io's surface, this test should allow us to determine the well-known rotational velocity of Io, before applying the procedure to Titan. A detailed analysis of the Titan observations is presented in a separate article by Luz et al. (2005).

In this paper, we present new data reduction and analysis procedures for the retrieval of the wind velocity in a planetary atmosphere, and we apply them to Io. Sections 2 and 3 give the method overview and describe the Io observations and instrument settings. The following two sections review all the data reduction steps and carefully address the random and systematic errors inherent to the method. The results for Io are discussed in Sect. 6.

   
2 Method overview

The induced Doppler shift in the reflected solar spectrum depends on both the motion of the planet (relative to the Earth and the Sun) and its rotational velocity, which results from both the rotation of the planet's surface and the wind flow. Assuming that the rotational velocity is purely zonal and depends only on latitude, the Doppler shift induced by rotation is maximal at the limbs of the planet for an observer at infinity, and null at the planet's center. On the other hand the Doppler shift induced by the relative motion of the body, resulting from a velocity offset, is constant for the whole planet. Thus, detecting a differential Doppler shift from two regions of the atmosphere observed simultaneously (not necessarily opposite limbs) will allow the measurement of the shift induced by rotation only. Then the wind velocity is derived by subtracting the planet's surface rotation velocity, which is well known.

Given the order of magnitude of the planetary wind velocity (a few hundred m s^-1) a detection in visible wavelengths of the Doppler shift on individual solar lines is not possible since it would require a spectral resolution greater than 2 million. However, current high-resolution spectrometers such as UVES or the HIgh Resolution Echelle Spectrometer (HIRES) at the Keck observatory reach a maximal spectral resolution of about 100 000, and at this resolution the solar lines (whose typical width is 150 mÅ) are well resolved. The wide spectral range covered by the instrument and the many lines it contains can therefore be used to apply the Absolute Accelerometry retrieval scheme of Connes (1985), which allows one to measure Doppler shifts much smaller than the line width with a precision, limited by photon noise, that can reach about 1 m s^-1. A comparison between current methods for measuring small changes in radial velocities is done in Sect. 4.7.

Grating echelle spectrometers are well suited to such a technique. The spectrometer's entry slit can be aligned with the planet's equator, or set perpendicular to the planet's spin axis at a given latitude (Fig. 1). Decreasing the slit width increases the spectral resolution, reduces the range of observed latitudes, but obviously reduces the incoming flux. The whole slit image is then dispersed by the echelle grating that builds sequences of monochromatic slit images called an echellogram, recorded on a pixel array detector. In this way each region of the slit image (i.e. each region of the planet's atmosphere along the slit) yields its own spectrum, the slit image sampling depending on the spectrometer design. It should be noted that this step requires an accurate two-dimensional spectrum extraction, which is not correctly performed by the standard UVES reduction package, as specified in the UVES Quality Control report available at the ESO website.

	Figure 1: Observational geometry in the test case of Io, where a thin slit is aligned with the equator of the planet. The relative motion of the body induces a velocity offset V0. Measuring differential Doppler shifts allows the retrieval of the rotational velocity v only.
Open with DEXTER

Finally, the rotational velocity is retrieved from the set of measured Doppler shifts, assuming that the latter must be symmetrical in relation to the planet's center. The extracted velocity value is a lower limit of the wind velocity when systematic errors are considerable, as in our case. To better understand the measurement technique, results can be compared to a model we developed, that takes the Earth's atmosphere degradation and the instrumental effects into account. This step is particularly useful when the apparent diameter of the target is close to the seeing value.

   
3 Observations

Four observations of Io were carried out on the same night, the 18th of February 2003 (see Table 1), with the UVES instrument mounted on the Nasmyth platform of the 8.2 m Kueyen unit at the VLT, Cerro Paranal. The observational geometry is shown in Fig. 1, with the east-west direction defined as viewed on the sky.

Table 1: Summary of the observations of Io. The last column shows Io's apparent range-rate relative to the Paranal Observatory.

A thin slit (dimension 0.3 $^{\prime\prime}$ $\times$ 8.9 $^{\prime\prime}$) was aligned with the equator of Io by setting the UVES derotator angle to 109$^\circ$, according to Io's known north pole position angle with respect to the direction of the celestial north pole. The UVES Red Arm and the cross-disperser #3 were used with the central wavelength set to 520 nm, allowing one to record in a single exposure a wavelength domain from 414 to 621 nm at a resolution of about 100 000. The set of echelle orders were recorded on two 2048 $\times$ 4096 CCDs, referred to as EEV (blue part) and MIT (red part), with a loss of two half-orders ($\sim$5 nm) in the gap between the CCDs. With this setting, the pixel scale in the cross-dispersion direction is about 0.18 $^{\prime\prime}$.

The seeing was estimated from surrounding stars by the Differential Image Motion Monitor (DIMM) at the Paranal Observatory to be contained between 0.70 $^{\prime\prime}$ and 0.86 $^{\prime\prime}$ at 500 nm, which is relatively large compared to the 1.2 $^{\prime\prime}$ apparent diameter of Io. The expected pointing accuracy of the VLT as measured from stars should reach 0.1 $^{\prime\prime}$. The target is centered automatically on the slit and the tracking of the telescope is corrected by an autoguiding facility. It is noteworthy that the procedure has been designed primarily for stars. Since the planets of the solar sytem are moving targets and extended sources, we expect larger pointing or tracking errors.

During the night, Io's topocentric velocity increases substantially with time (up to 5 km s^-1). Consequently we cannot sum the 4 exposures, since the induced Doppler shifts will be very different. If there is no interruption during an exposure the target acceleration stretches the solar lines slightly, however, the differential Doppler shift induced by rotation remains the same. Thus, each exposure should lead to the well-known rotational velocity of 74.6 m s^-1, which is not significantly modified by the inclination of the equatorial plane since the sub-Earth point latitude is almost null ($\sim$0.13 $^\circ$).

In addition to the observations of Io a thorium-argon lamp exposure was taken immediately after the science exposure to provide an accurate wavelength calibration frame, minimizing the effects of changing temperature, pressure or of a small earthquake. Such caution is mandatory since we expect shifts smaller than 0.2 pixel.

   
4 Data reduction

4.1 Objectives

An echelle order is a two-dimensional spectral image built from a sequence of one-dimensional monochromatic slit images (Fig. 2). Standard algorithms are designed to extract a one-dimensional object spectrum from the two-dimensional spectral image of a point source, usually by using the optimal extraction algorithm introduced by Horne (1986). The latter applies a non-uniform pixel weight to sum the monochromatic slit image in order to maximize the signal-to-noise ratio and eliminate bad pixels. In the case of an extended object, by definition incompatible with optimal extraction, the extraction of a two-dimensional object spectrum does not allow one to maximize the signal-to-noise ratio and raises specific issues which are described in the next sections.

Note that the whole data reduction procedure described in this section can be performed by the WIND MATLAB package we developed. The source code can be obtained on request by contacting the first author.

	Figure 2: Each order is a spectral image built from a sequence of monochromatic slit images. The curvature of both the spectral and slit images needs to be corrected.
Open with DEXTER

4.2 Required data

To run the reduction procedure, a set of calibration data is required in addition to the target's spectrum, hereafter "science frame'', preferably in FITS (Flexible Image Transport System) format. This includes bias frames (non-zero signal at zero flux), flat field frames (white lamp exposure), and a wavelength calibration frame (e.g. a thorium-argon lamp exposure). Optional data may include an order definition frame (flat field obtained with a thin slit) which we use to detect possible drifts of the target during an observation, and a slit viewer image taken at the beginning of spectroscopic exposures, which allows an estimation of the position of the slit on the planet's disc.

Since our own wavelength calibration does not improve the one performed by the UVES pipeline (Ballester et al. 2003), we use the wavelength calibration table provided by ESO. Improving the line identification does not allow us to reach a sub-pixel scale accuracy; to complete this task we use the statistical power of the Connes algorithm.

4.3 Preliminary processing

Images are imported by a parametric procedure which cuts unwanted regions in the echellogram and sets its orientation with wavelength increasing left to right in the dispersion direction, and bottom to top in the cross-dispersion direction. Considering that there is a wavelength overlap between adjacent orders of a few nm and that the signal-to-noise ratio is very low at the beginning of the orders, we reject about 10% of each order to improve the robustness of the overall data reduction. The intensity is converted from ADU (Analog Digital Unit) to photoelectrons by using the gain value written in the FITS header ($\sim$1.5 e^-/ADU).

A median bias frame is computed from the available bias frames (usually 5 per night), which is subtracted from all non-bias images. A median flat field frame is computed in the same way. Pixel-by-pixel median filtering is used in order to eliminate cosmic hits.

4.4 Mapping order locations

Echelle orders are curved and tilted in both the dispersion and cross-dispersion directions (Fig. 2). In the case of UVES, the inclination of the spectral image is a few degrees and the slit image curvature reaches a maximum value of one pixel ($\sim$1 km s^-1) at the extremes (Ballester & Rosa 1997). A correction of the curvature with a sub-pixel sensitivity is therefore required, which will be performed later by using the Connes algorithm (see Sect. 5.2).

The algorithm that determines spectral image curvature and location consists of four steps. First, the procedure finds edges in the image where the gradient of intensity is maximum along the cross-dispersion direction, and returns a binary mask equal to unity where maxima occur. Then the Radon transform of the binary mask is computed, which is maximal when the projection line is aligned with an order. The coordinates of the projection lines are used as first guesses by a threshold filter set to half of the maximal intensity, which estimates for each monochromatic slit image the top of the order spatial profile (Fig. 3). The estimate is the mean position of the pixels with intensity greater than the threshold. The spatial profile is very close to a Gaussian function and we will characterize it in the next sections by its full width at half maximum (FWHM).

	Figure 3: Spatial profile of an order. The latter is close to a Gaussian function in the case of Io and Titan. By convention, the region corresponding to the centroid of the intensity distribution is r=0.
Open with DEXTER

We define the target image center (i.e. the region recorded on pixel zero) where the slit image intensity is maximal. It should be noted that usually the slit image center does not coincide with the target image center. For this reason, order location mapping should be performed with the science frame and not with the order definition frame. Finally, for each order the set of spatial profile centers is fitted by a polynomial of tunable order (set to 3 for a good estimation), yielding a polynomial description of order center locations.

4.5 Scattered light correction

Science and flat field frames both must be corrected for scattered light and sky emission. Since the spectral image itself includes the background, it precludes a direct measurement beneath the orders. However, inter-order measurements and interpolation provide an adequate estimate of the background.

This step is performed by using the scheme given by Piskunov & Valenti (2002). The procedure computes the position of background regions (about 10 pixels wide) using the order polynomial description. For each column (cross-dispersion direction), medians are calculated above and below each order. The resulting set of background values is softly filtered in the dispersion direction by performing a convolution with a thin Gaussian function. Finally, linear interpolation in the cross-dispersion direction is used to estimate the scattered light component of the spectral image.

The scattered light intensity is about 300 e^- per pixel for an exposure time of 140 s. It is very small compared to the spectral intensity at the center of the planet (where the flux intensity reaches 5 $\times$ 10⁴ e^-), but the ratio (ruled by the Gaussian-like spatial profile) becomes lower and lower when considering pixels far from the disc center, until the background intensity becomes larger than the signal from the target. Moreover, since the background correction does not remove the noise induced by scattered light (photon noise), it is therefore important to take the background noise into account, particularly for regions far from the planet's center.

4.6 Spectrum extraction

For each order, a two-dimensional spectrum is extracted from the echellogram by reading the intensity value of the spectral image on both sides of the previously determined target image centers. It should be noted that target image centers are real numbers which most often do not match the pixel CCD grid. Moreover, assuming the rough sampling of the slit image in the cross-dispersion direction, the Gaussian-like spatial profile (in the case of Io or Titan) and the order inclination of a few degrees, extracting the nearest recorded pixel value leads to periodic breaks in the spectrum, corresponding to the changes of row (cross-dispersion direction) as shown in Fig. 4. This effect is similar to fitting a tilted line with a step function and depends on the derivative of the spatial profile (hence it is maximal in the Gaussian wings). We checked with Fourier analysis that the periodicity of the breaks is equal to the cotangent of the order slope, a slope of 2$^\circ$ leading to a period of about 30 pixels.

	Figure 4: Breaks due to the CCD grid periodically appear in the extracted spectra when no interpolation is used.
Open with DEXTER

In order to avoid such features and to keep the accurate position estimate of the target image center, the extraction algorithm performs an interpolation of the spectrum values. The interpolation method is chosen according to the spatial profile, which depends, among other parameters, on both the source type and the slit dimension. Thus, we use a spline interpolation for the spectrum extraction of Io (Gaussian-like profile), and a linear interpolation for the extraction of the flat field and the thorium-argon lamp (gate-like profile).

During the interpolation, known CCD bad pixels are ignored with a mask, which can also take cosmic hits into account. Finally, in order to correct the pixel-dependent CCD quantum efficiency, each extracted science spectrum is divided by the flat field spectrum (both previously corrected for scattered light). From this step on, the flux is no longer given in photoelectrons. Thus, we compute the correct expected photon noise before the flat field correction by considering the science spectrum and scattered light intensity in photoelectrons.

Note that merging spectra would be useless for our purpose since the Connes technique is not a cross-correlation-like algorithm. Doppler shift will be measured for each wavelength (i.e. pixel), independently of the other values. Moreover, merging orders implies a loss of information since "wavelength overlap'' is not equivalent to "information overlap'', the photons relative to adjacent orders being independent.

   
4.7 Doppler shift measurement

Thanks to the search for extra-solar planets and stellar seismology, much progress has been accomplished in techniques dedicated to measuring small changes in radial velocities. The existing methods mainly differ depending on the wavelength calibration scheme and on the procedure used to estimate the variation of velocity. The calibration spectrum is obtained either separately from the target's spectrum (but simultaneously to the observation), or superposed on it by placing a iodine cell at the spectrometer's entrance slit, see Baranne et al. (1996) and Butler et al. (1996) for details. On the other hand, the velocity changes can be estimated by a procedure based on a cross-correlation technique, or by using a procedure that computes an optimum weight, such as the Connes algorithm or a weighted analysis in the spectral Fourier space as described by Chelli (2000).

A thorough comparison between these methods has been done by Bouchy et al. (2001). They found that the optimum weight procedure is a factor 1.6 more efficient than the cross-correlation technique (by running a Monte Carlo noise simulation), and that the uncertainty on velocity measurements done by the iodine cell technique is greater by a factor of 2 to 2.4 compared to the separate calibration method. According to these conclusions, we use the algorithm proposed by Connes (1985) and a thorium-argon lamp for the calibration, performed subsequently to the target exposure since the UVES design precludes us from obtaining two spectra simultaneously.

The Connes technique uses the full spectral information available to compute Doppler shifts much smaller than the line width. Let us consider a digitized and calibrated reference spectrum, ideally noise free. For a given pixel i along the dispersion direction, the intensity  $I_{\rm ref}(i)$ and the wavelength  $\lambda(i)$ are known. Let us also consider the same spectrum slightly shifted whose intensity is I(i). Assuming that the Doppler shift  $\delta \lambda$ is small compared to the line width, the velocity change $\delta v$ obtained at the pixel i can be expressed by

$${\delta v}_i = \frac{c}{\lambda(i)}~ {\delta \lambda}_i = \... ...\rm ref}}{\partial \lambda}\right)_{\lambda=\lambda(i)}}\cdot$$ (1)

Consequently, the change of velocity is measured through a change of intensity. To increase the sensitivity of the measurement the whole available spectral range must be used. Hence, the contribution of all wavelengths (i.e. all pixels) is summed considering an optimum weight w_i equal to the inverse square of the standard deviation, which gives

$$\delta v = \frac{\sum {\delta v}_i~w_i}{\sum w_i} \quad \textrm{where} \quad w_i = \frac{1}{{\sigma[{\delta v}_i]}^2}\cdot$$ (2)

It should be noted that the spectrum $\bar{I}(i,r)$, extracted from the region r along the slit (r=0 at center), has a continuum level smaller than the one of the planet's center spectrum (Figs. 3 and 5), thus we need to normalize the set of spectra in order to have

$$I(i,r) \sim I_{\rm ref}(i) = K(r) ~\times~ \bar{I}(i,r).$$ (3)

For each order, the normalization function K(r) is equal to the inverse of the mean spatial profile. Thus, assuming Eqs. ((1)-(3)) the optimum weight for the spectrum extracted from region r is given by

$$w_i = \frac{{\left( \frac{\lambda(i)}{c} \right)}^2 {\left(... ...f}(i) + ~K^2(r)~B(i) + \left(1 + K^2(r)\right)\sigma_{\rm R}}$$ (4)

where B(i) is the background intensity estimated at pixel i and  $\sigma_{\rm R}$ is the detector read-out noise equal to a few electrons. We consider that both the reference and the shifted spectrum are noisy. The noise results from the quadratic sum of the detector noise and the photon noise, the latter being equal to the square root of the intensity in photoelectrons. Moreover, since I and  $I_{\rm ref}$ have the same intensity level and the Doppler shift is small, we set  $I(i)=I_{\rm ref}(i)$ in Eq. (4).

	Figure 5: Segment of a two-dimensional spectrum. Each pair corresponds to the eastern and western regions of the target, intensity decreasing with distance from the planet's center.
Open with DEXTER

We set the reference spectrum as the sum of all the extracted spectra weighted by the spatial profile of the order. In this way the reference spectrum is mainly given by the spectrum from the planet's center and we improve the signal-to-noise ratio.

Finally, the weighting function is modified to eliminate unwanted lines, such as telluric absorption lines. The Earth's telluric line mask was obtained by observing the nearly featureless spectrum of a rapidly rotating B star (Martic et al. 2001). However, the use of such a mask did not lead to perceptible differences when computing the Doppler shifts. On the other hand, in the case of Io it was necessary to remove sodium emission line contamination, since sodium D-lines at 589.0 and 589.6 nm are strongly present all along the slit, hence still detectable at seven Io's radii.

5 Errors

5.1 Doppler shift uncertainty

The error estimate of the Doppler shift value retrieved with the Connes algorithm can be computed. According to Eq. (2) the uncertainty on the velocity change is given by

$$\sigma[\delta v] = \frac{1}{\sqrt{\sum w_i}}\cdot$$ (5)

The error depends on the photon noise and the spectral line richness of the spectrum. It is the main source of random error. In the case of the observations of Io, we achieve an accuracy smaller than 2 m s^-1 for the regions near the center of the disc.

   
5.2 Slit image curvature

In a grating echelle spectrometer, let us consider the chief ray from the object point to the vertex of the grating surface, and let $\gamma$ be the angle it makes with the plane xz perpendicular to the y grating grooves axis. The projection on the xz plane of the angle $\gamma'$ of the refracted chief ray makes an angle $\beta$ with the z axis. Due to the dependence of $\beta$ on $\gamma$ at constant wavelength, the monochromatic image of a long straight slit parallel to the y axis is curved (Schroeder 1987,  p. 261). The slit image curvature, equal to tenths or hundredths of the pixel width, induces errors affecting the wavelength calibration at sub-pixel scale. Such errors cannot be corrected with standard algorithms based on line identification. Nevertheless, if uncorrected, the curvature leads to errors larger than 50 m s^-1, as shown in Fig. 6.

	Figure 6: The slit image curvature is computed by using the Connes algorithm. The curvature is almost the same for each order, as shown by the dotted lines.
Open with DEXTER

In order to correct the artificial shift of spectra due to the slit image curvature, we extract the two-dimensional spectrum of the wavelength calibration lamp and we run the Connes algorithm as for the science frame. Such a technique works well if the curvature is almost constant, since the Connes algorithm is able to detect a constant Doppler shift, but our results show that the variation is very small. For each order and both CCDs, the curvature is almost the same (dotted lines in Fig. 6). Thus, we compute the sum of the slit image curvatures of each order weighted by the error estimate. The achieved accuracy is about 1 m s^-1.

Finally, the pseudo Doppler shift value resulting from the image curvature and given by the calibration lamp is subtracted from the Doppler shift value computed from the science observation. In this way the lamp, that is obviously motionless in relation to the slit, has a measured velocity equal to zero.

5.3 Earth's atmosphere

Depending on the apparent diameter of the target and atmospheric conditions, the degradation due to the Earth's atmosphere can have important consequences on the measurement by mixing incoming photons. We thus have to consider atmosphere turbulence and refraction, which are both wavelength dependent.

	Figure 7: The FWHM variation of the spatial profile mainly depends on the seeing and possible drifts of the target during the exposure.
Open with DEXTER

The atmospheric turbulence is given by the seeing value, which is equal to about 0.75 $^{\prime\prime}$ during the observations of Io. Note that this value is computed from stars (point-like sources) and may not be fully applicable to an extended source. Assuming the standard model (Fried 1965) the image FWHM has a weak  $\lambda^{-0.2}$ dependence on wavelength, the seeing should therefore decrease when the wavelength increases. This trend is observed in Fig. 7 on the two CCDs by measuring the FWHM (in pixels) of the spatial profile of each order. The shift between the CCDs is due to the small variation of pixel scale for each CCD. Actually we measured an EEV pixel scale about 1.7% wider than the MIT pixel scale, leading to a continuous decrease of the seeing. Moreover, the abnormal increase of the seeing that occurs for the last MIT orders of the Io829 and Io767 observations is linked to the saturation of the CCD. The exposure time was slightly overestimated and the orders from 104 to 100 with the highest flux intensity start to saturate, thus yielding an erroneous FWHM value when performing a Gaussian fit. It is noteworthy that the slit image FWHM depends also on possible target drift (as discussed in the next section) and, therefore, we cannot deduce the seeing from the FWHM value.

	Figure 8: Relative position of the orders. The centering variation along the slit length is larger than the pixel scale for two consecutive observations.
Open with DEXTER

The Earth's atmospheric refraction is another source of image degradation. However, the UVES Atmospheric Dispersion Corrector unit (ADC) is inserted in the light path in order to compensate for this effect. Along the slit length direction, the refraction effect can be measured by computing the residual distance between the slit image center (given by the order definition frame) and the target image center. If the image position depends on wavelength (i.e. order number), the residual distance will not be constant. Figure 8 shows that along this direction and after the ADC correction, the measured refraction effect is smaller than 0.05 $^{\prime\prime}$. This value is consistent with the expected ADC maximal residuals, which is equal to about 0.07 $^{\prime\prime}$ (along both directions) at the central wavelength of 520 nm.

   
5.4 Pointing accuracy

The following discussion assumes that the ephemerides used by the telescope pointing facility are perfect. However, that is not the case since our calculation method was not optimal. Therefore, the quality of the ephemerides could also be responsible for the pointing issues.

Both the study of the order spatial profile and the determination of the target image position show evidence of target drift during exposures. Moreover, the slit images taken by the UVES slit viewer show that the target is not well centered on the slit in both directions, but, since the exposure time is short (1 ms), it could also result from a seeing effect.

Nevertheless, the mean position of the target center from the slit center during an exposure can be determined by measuring the residual distance of orders from the ones given by the order definition frame (Fig. 8). It shows that the centering variation along the slit length reaches 0.22 $^{\prime\prime}$ for two consecutive observations, which is larger than the pixel scale.

In addition, the spatial profile variation (Fig. 7) shows that the orders for the Io615 observation are much too wide compared to other observations. Moreover, despite the fact that the Io767 observation has the worst seeing value with 0.86 $^{\prime\prime}$ (and seeing "spreads'' the spatial profile), it is one of the observations with the thinnest orders. These effects very likely result from a target drift.

	Figure 9: A set of Doppler shifts is computed for each order (dotted lines). We sum independently the orders from the blue part (triangles) and red part (squares) of the spectral range. The slope in the central part corresponds to the direct rotation of Io, as expected.
Open with DEXTER

Thus, a pointing accuracy ideally smaller than the pixel scale is not easily obtained. We expect variations in both target latitude and longitude directions. A latitude variation (slit width direction) can be an issue since the maximal velocity at the limbs and the apparent diameter of the target depends on the pointed latitude. A longitude variation (slit length direction) is less problematic since the target image center is determined where the flux is maximal. Finally, a possible drift during an exposure has consequences that are poorly predictable since the drift direction is unknown. The degradation of the target image by such an effect may lead, in the extreme case, to a non-existant Doppler shift between the extracted spectra.

   
6 Data analysis

   
6.1 Results

The 4 observations of Io have been reduced according to the procedure described in Sect. 4. The result for the Io829 observation is plotted in Fig. 9. Doppler shift measurements are performed independently for each order (dotted lines in the figure), to check the consistency of the results obtained with the Connes algorithm and to detect a possible wavelength dependence. Such a dependence is visible for the two CCDs, corresponding to the blue and red part of the spectral range. Thus, we sum independently the set of orders for each CCD (squares and triangles in Fig. 9), assuming a weight equal to the inverse square of the standard deviation.

The computed Doppler shift curves exhibit a similar shape for the 4 observations, but each has a different slope in the central part and a different maximal Doppler shift value, the Io829 observation having the largest value. Moreover, the blue part of the spectral range systematically gives shifts smaller than those obtained from the red part. Finally, for all the observations, the direct rotation of Io is undoubtedly obtained by the measurement.

It should be noted that the window extraction we use is very wide (30 pixels is almost 5 times Io's diameter). The signal-to-noise ratio decreases quickly from the center of Io and we do not expect to retrieve reliable information too far from the central region. Doppler shift values measured outside the target disc may seem to be unexpected at first glance and we therefore developed a model (Sect. 6.2) in order to have a correct understanding of the measurement.

Moreover, we checked that the experimental standard deviation is consistent with the expected error bars. For this, we compared the error computed by the Connes algorithm for an individual order at a given pixel with the standard deviation of the set of values given by all the orders. The equality of the two quantities is verified up to about the seventh pixel from the center; unexpected effects might occur beyond this limit.

   
6.2 Measurement model

Let o(x) be the flux distribution of the source along the slit length, the source having a radius R centered at the origin, and let P(x) be the point spread function (PSF) relative to the atmospheric degradation and instrumental response. In the region r (i.e. pixel r along the slit), the number of photons coming from the region x of the source is

$$n_r(x) = \int_{\Delta_r} P(x-t)~o(x)~{\rm d}t \quad\quad x \in [-R,R]~\textrm{in arcsec}$$ (6)

where the integration range $\Delta_r$ (in arcsec), equal to the pixel scale and centered at the pixel r, is given by

$$\Delta_r = [\textrm{pixel scale}\times(r - 0.5),~\textrm{pixel scale}\times(r + 0.5)].$$ (7)

Equation (6) gives the spatial origin of the photons. Note that n_r(x) is not constant for each pixel since the number of photon depends on the flux distribution which tends to zero at the limbs. Nevertheless, in order to compare spectra with each other we normalize the spectrum (see Eq. (3)), which is equivalent to considering the normalized distribution

$$\bar{n}_r(x) = \frac{n_r(x)}{\int_{\Lambda} n_r(x)~{\rm d}x} \quad\quad \textrm{with}\quad \Lambda = [-R,R].$$ (8)

The mean Doppler shift $\langle s(r) \rangle$ measured at pixel r depends on the velocity distribution v(x) of the source, hence

$$\langle s(r) \rangle = \int_{\Lambda} \bar{n}_r(x)~v(x)~{\r... ...{\Lambda} \int_{\Delta_r} P(x-t)~o(x)~{\rm d}t~{\rm d}x}\cdot$$ (9)

Note that this result is not inherent to the Connes algorithm and is valid for any method used to measure Doppler shifts.

	Figure 10: The measured Doppler shift ideally tends towards a maximal velocity value at infinity, smaller than the target's rotational velocity. An artificial decrease may occur due to normalization error.
Open with DEXTER

The model can be applied in the case of Io. We use a cosine-like flux distribution equal to zero outside Io's disc, and a PSF proportional to a 5/3 power exponential, assuming the standard seeing model. We need to take a double Doppler shift effect into account since the target is both a light receiver in relation to the Sun and a light transmitter in relation to the Earth. Let $\theta$ be the longitude relative to the center of the disc as seen from the Earth, at $\theta$ the velocity relative to Earth is V. Let $\theta'$ and V' be theses quantities as seen from the Sun. The induced Doppler shift results from

$$V + V' = (V_0 + v~\sin \theta~\cos \delta) + (V_0' + v~\sin \theta'~\cos \delta')$$ (10)

where $\delta$ (resp. $\delta'$) is the sub-Earth (resp. sub-Solar) point latitude, V₀ (resp. V₀') is the velocity induced by the motion in relation to the Earth (resp. the Sun), and v is the rotational velocity at the limbs. Assuming $\delta \sim \delta'$ and $\theta = \phi + \theta'$ where $\phi$ is the Earth-Target-Sun phase angle, the measured differential Doppler shift relative to the disc center results from

$$\Delta V = v~\cos \delta \left( \sin \theta + \sin \phi + \sin (\theta - \phi) \right).$$ (11)

Note that $\Delta V$ depends only on the rotational velocity as desired, and when the $\phi$ and $\delta$ angles are close to zero, as in the case of Io, Eq. (11) leads to a maximal Doppler shift simply equal to  $2 \times v$ at the limbs.

The results for three different simulations are plotted in Fig. 10. The figure shows that the slope of the curve in the central region is much smaller than the one given by the velocity distribution. The slope decreases when the seeing increases. Moreover, the measured Doppler shift ideally tends towards the maximal rotational velocity at infinity, without reaching it. This set of effects is consistent with the mix of photons which becomes more important when the seeing increases. Furthermore, since the normalization is performed by an algorithm, a residual error is possible. The limited accuracy of this step induces an artificial decrease of the measured Doppler shift value for pixels far from the center, where the normalization factor is large, as shown by the simulation (circles in Fig. 10).

However, despite the similar trends given by the model and by the measurement, experimentally the seeing effect is not as strong as predicted, particularly for regions outside the disc of the target. It seems that an accurate prediction would very likely require a more accurate seeing model, including the large spatial scales.

   
6.3 Retrieved velocity

We extract a distinct velocity value for the spectral range of each CCD, and for each observation. For this, we first determine the best point of symmetry (in the least square sense) of the velocity curve, since usually the latter is not precisely centered on the frame. Then we determine the domain where the shifts are symmetrical in relation to the estimated planet's center, assuming a tolerance proportional to the error estimate of the shift ($\sim$$5\sigma$). The velocity is finally given by the maximal Doppler shift in this range (Fig. 11), which usually occurs close to the end of the domain. The velocity is therefore measured consistently for a continuous set of values from the target's center, considering non-symmetrical values as spurious. In this way the retrieved velocity is a lower limit since atmospheric degradation or target drift, with a sufficient magnitude, may reduce the maximal detectable Doppler shift. It is noteworthy that the Doppler shift values beyond the retrieved velocity present a low signal-to-noise ratio which becomes quickly smaller than 5 (dashed line in Fig. 11). Moreover, the set of consistent values belong to the range where error bars have been validated (as explained in Sect. 6.1). The $4 \times 2$ extracted velocities, corrected for the double Doppler effect given by Eq. (11), are plotted in Fig. 12.

	Figure 11: Velocity retrieval scheme. The extracted velocity is the maximal Doppler shift value measured in the range where the shifts are symmetrical in relation to the estimated target's center.
Open with DEXTER

	Figure 12: Summary of the results for Io. Velocities extracted from the blue and red spectral ranges are plotted with triangles and squares. The spurious variation results from systematic errors, such as the uncertainty on the pointed latitude, drifts, or the degradation inherent to the Earth's atmosphere.
Open with DEXTER

The red part of the Io829 and Io615 observations give the expected rotational velocity of Io with an uncertainty smaller than 2 m s^-1, the blue part gives values 7% and 20% smaller. The wavelength dependence of the seeing can explain the variation for Io829, but an indeterminate additional effect is probably involved regarding Io615. The Io315 observation shows evidence of pointing errors (as expected in Sect. 5.4) and the low velocity could be explained by an observed latitude different from Io's equator (e.g. 55$^\circ$ instead of 0$^\circ$). A comparison with the Io829 observation shows that seeing would permit a correct measurement at least in the red spectral range. In the case of the Io767 observation, both the seeing and pointing errors could explain the velocity value 15% smaller than the one expected at the equator.

On the one hand random errors are consistent with our expectation and are fully understood. On the other hand some systematic errors are too large, probably due to pointing issues and the atmospheric turbulence. In the best cases (2 times), we measure very accurately the known surface velocity of Io. Thus, the principle of the method is validated, but a large number of measurements should preferably be carried out in order to identify spurious variations resulting from systematic errors.

7 Conclusion

We presented a new method that uses high-resolution spectroscopy in the visible to measure planetary winds. The zonal wind flow in a planetary atmosphere reflecting the solar light induces spatially dependent Doppler shifts in the reflected solar spectrum. Using a cross-dispersed echelle spectrometer and a dedicated spectral analysis allows a direct determination of both the wind speed and direction.

Since the UVES pipeline procedures were not applicable, we developed a dedicated package to perform the data reduction. The latter consists of extracting a two-dimensional spectrum from the target image, and then comparing Doppler shifts induced by different regions of the target's atmosphere. Since the velocities involved are in the range 100-200 m s^-1, Doppler measurements are performed with the Connes algorithm, which is an optimum technique using the full available spectral information.

To validate the method we applied it to Jupiter's moon Io, whose apparent diameter and surface rotational velocity are similar to the quantities involved for Titan. In the best case,we retrieved the well-known surface velocity of Io with an uncertainty smaller than 2 m s^-1, consistent with the expected random errors. Nevertheless, this test case showed that taking the Earth's atmospheric degradation into account is essential and that some systematic errors are too large and very likely due to pointing issues. Therefore, one needs to consider the retrieved velocity as a lower limit.

Finally, the method has been applied to Titan and the results are discussed in a separate paper by Luz et al. (2005).

Acknowledgements

T.C. thanks ESA/RSSD for having initiated this work; he is currently working at LESIA under ESA financial support (Huygens funding). We are grateful to N. Piskunov for providing his REDUCE package. We thank H. Dekker and A. Kaufer who answered our questions about the technical aspects of VLT/UVES. The authors are also grateful to E. Lellouch for his helpful comments, and to M. Martic for sharing the Earth's telluric lines table. We acknowledge M. Reyniers for useful discussions at the beginning of this work. DL acknowledges financial support from Fundação para a Ciência e a Technologia.

References

• Ballester, P., Boitquin, O., Modigliani, A., & Wolf, S. 2003, in UVES Pipeline User's Manual. Issue 5 (V2.0.0)
• Ballester, P., & Rosa, M. R. 1997, A&AS, 126, 563 [NASA ADS]
• Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS]
• Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733 [EDP Sciences] [NASA ADS]
• Butler, R. P., Marcy, G. W., Williams, E., et al. 1996, PASP, 108, 500 [NASA ADS] [CrossRef]
• Chelli, A. 2000, A&A, 358, L59 [NASA ADS]
• Connes, P. 1985, Ap&SS, 110, 211 [NASA ADS]
• Flasar, F. M., Samuelson, R. E., & Conrath, B. J. 1981, Nature, 292, 693 [NASA ADS]
• Fried, D. L. 1965, Opt. Soc. Am. J., 55, 1427 [NASA ADS] [MathSciNet]
• Horne, K. 1986, PASP, 98, 609 [NASA ADS] [CrossRef]
• Hubbard, W. B., Sicardy, B., Miles, R., et al. 1993, A&A, 269, 541 [NASA ADS]
• Kostiuk, T., Fast, K. E., Livengood, T. A., et al. 2001, Geophys. Res. Lett., 28, 2361 [NASA ADS]
• Lebreton, J.-P., & Matson, D. L. 2002, Space Sci. Rev., 104, 59 [NASA ADS]
• Luz, D., Civeit, T., Courtin, R., et al. 2005, in Characterization of the zonal wind flow in the stratosphere of Titan with UVES, Icarus, submitted
• Martic, M., Lebrun, J.-C., Schmitt, J., Lebreton, J.-P., & Appourchaux, T. 2001, in Feasibility demonstration of Titan wind measurement technique at OHP using absolute accelerometry method, Internal ESA report
• Martic, M., Schmitt, J., Lebrun, J.-C., et al. 1999, A&A, 351, 993 [NASA ADS]
• Piskunov, N. E., & Valenti, J. A. 2002, A&A, 385, 1095 [EDP Sciences] [NASA ADS]
• Schroeder, D. J. 1987, Astronomical optics (San Diego, CA: Academic Press, Inc), 363