M Gaudesius

Y.-C Zhang

T Pohl

R Kaiser

G Labeyrie

Phase diagram of spatiotemporal instabilities in a large magneto-optical trap

Large clouds of cold atoms prepared in a magneto-optical trap are known to present spatiotemporal instabilities when the frequency of the trapping lasers is brought close to the atomic resonance. This system bears similarities with trapped plasmas where the role of the Coulomb interaction is played by the exchange of scattered photons, and astrophysical objects, such as stars whose size is dependent on radiative forces. We present in this paper a study of the phase space of such instabilities and reveal different dynamical regimes. Three-dimensional simulations of the highly nonlinear atomic dynamics permit a detailed analysis of the experimental observations.

I. INTRODUCTION

Since its first realization in 1987 [START_REF] Raab | Trapping of Neutral Sodium Atoms with Radiation Pressure[END_REF], the magneto-optical trap (MOT) has been a very successful device for the preparation of dilute and cold atomic ensembles and is nowadays present in many research laboratories throughout the world. In addition, in the so-called "multiple scattering regime" of the MOT [START_REF] Walker | Collective Behavior of Optically Trapped Neutral Atoms[END_REF], the nonlinear dynamics arising from the light-induced atomic interactions has attracted broad fundamental interest [START_REF] Labeyrie | Self-Sustained Oscillation in a Large Cloud of Cold Atoms[END_REF][START_REF] Pohl | Self-Driven Nonlinear Dynamics in Magneto-Optical Traps[END_REF][START_REF] Mendonça | Collective oscillations in ultracold atomic gas[END_REF][START_REF] Mendonça | Photon Bubbles in Ultracold Matter[END_REF] and bears intriguing analogies to one-component plasmas [START_REF] Hansen | Statistical mechanics of dense ionized matter. III. Dynamical properties of the classical one-component plasma[END_REF]. Compared to other related systems, such as Cepheid stars [START_REF] Cox | Theory of Stellar Pulsation[END_REF], the MOT presents the advantage to be highly tunable: The number of trapped particles can be adjusted as well as the parameters that determine the magnitude of the effective charge.

Spatiotemporal instabilities in various MOT geometries have been studied for 2 decades [START_REF] Labeyrie | Self-Sustained Oscillation in a Large Cloud of Cold Atoms[END_REF][START_REF] Wilkowski | Instabilities in a Magneto-Optical Trap: Noise-Induced Dynamics in an Atomic System[END_REF][START_REF] Di Stephano | Deterministic instabilities in the magneto-optical trap[END_REF]. We investigate in this Letter the instabilities of the "balanced" MOT where the three pairs of laser beams used to generate the trapping and cooling force in three spatial dimensions are produced with six independent intensity-balanced laser beams (as opposed to three retroreflected beams). In a previous work [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF], we specifically addressed the issue of the instability threshold which was experimentally studied both at fixed trapped atom number N and as a function of N. The observed behaviors were compared to the predictions of different models. In the present Letter, we are interested in analyzing the cloud's deformation modes in the unstable regime. We construct a phase-space diagram of these modes by scanning two of the main MOT parameters, the laser detuning, and the magneticfield gradient. We identify several distinct regimes for the unstable MOT with a particularly interesting regime occurring at low values of the magnetic-field gradient that generates the MOT's confinement. The observed behaviors are compared to three-dimensional numerical simulations.

II. EXPERIMENTAL SETUP AND NUMERICAL SIMULATIONS

The characteristics of our 87 Rb MOT have been described elsewhere [START_REF] Camara | Behavior of a very large magneto-optical trap[END_REF]. Its main feature is the ability to trap large atom numbers (up to 10 11 ), which enables us to reach a regime where collective effects are very strong. Producing such large cold atom samples is achieved by combining a MOT loading from an ambient room-temperature vapor with a large trapping volume. This volume is defined by the intersection of six independent laser beams of waist 3.4 cm that are intensity balanced. The typical peak intensity is 5 mW/cm 2 per beam with a total power of ≈ 550 mW in the MOT. A pair of anti-Helmoltz coils produces the magnetic-field gradient needed to generate the restoring force of the trap, whose typical value is 7 G/cm along the axis of the coils. Three pairs of compensation coils allow to finely tune the position of the zero of the magnetic trap, defining the MOT center. The typical loading time constant of the MOT is 2 s.

Our aim is to study the instability phase diagram when two of the MOT parameters, the laser detuning δ and the magnetic-field gradient ∇B, are varied whereas the others (the laser intensity I and trapped atom number N) are kept constant. Since the steady-state population in the trap depends on δ and ∇B (see, e.g., Ref. [START_REF] Camara | Behavior of a very large magneto-optical trap[END_REF]), we have to use the temporal sequence described in Ref. [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF]. It consists of two successive phases, a "loading" phase where the parameters (here δ and the loading time) are adjusted to obtain the desired N, and an "instability" phase where δ is rapidly changed to a target value. The magnetic-field gradient is set to the desired value and kept constant during the whole sequence. A fluorescence image is recorded at the end of this phase using the MOT laser beams with a fixed detuning of -8 (where is the natural width of the employed F = 2 → F = 3 transition) to minimize multiple scattering. The duration of the instability phase is 100 ms, long enough to decorrelate the cloud dynamics from the initial kick and short enough so that the MOT loading during the instability phase can be neglected. We probe the dynamics of the cloud in the instability phase by repeating this temporal sequence many times. Several analytical models have been proposed in the past to predict the threshold of spatial instabilities in MOTs [START_REF] Labeyrie | Self-Sustained Oscillation in a Large Cloud of Cold Atoms[END_REF][START_REF] Mendonça | Photon Bubbles in Ultracold Matter[END_REF][START_REF] Terças | Driven collective instabilities in magneto-optical traps: A fluid-dynamical approach[END_REF]. They all rely on different assumptions using a macroscopic ansatz and do not predict the same thresholds. In order to be less dependent on such simplifying assumptions, we have developed a numerical simulation, which, in principle, should include the various analytical predictions. In contrast to a first numerical approach described in Ref. [START_REF] Pohl | Self-Driven Nonlinear Dynamics in Magneto-Optical Traps[END_REF] which confirmed one of the analytical models, we here use a three-dimensional simulation based on microscopic ingredients of atom-light interactions [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF]. It includes the usual Doppler trapping and cooling force plus the two collective forces that are known to play an important role in the multiple scattering regime of the MOT [START_REF] Walker | Collective Behavior of Optically Trapped Neutral Atoms[END_REF]. The first one is the compressive force due to the attenuation of the MOT's laser beams through the cloud (the so-called "shadow effect") [START_REF] Walker | Collective Behavior of Optically Trapped Neutral Atoms[END_REF][START_REF] Dalibard | Laser cooling of an optically thick gas: The simplest radiation pressure trap?[END_REF]. The second one is the repulsive, effective Coulombian interaction force arising from exchange of scattered photons between trapped atoms [START_REF] Walker | Collective Behavior of Optically Trapped Neutral Atoms[END_REF]. It is the interplay between these two collective forces of opposite sign that gives rise to the unstable behavior of the trapped atomic cloud. Some details about the model and simulation are provided in Sec. III B of Ref. [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF]. In this Letter, the reported numerical results were in qualitative agreement with experimentally observed scaling laws, which gave us some confidence that the relevant physical ingredients were correctly included.

III. EXPERIMENTAL OBSERVATIONS

We report in Fig. 1 The dashed line materializes the position of the instability threshold [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF]. The cloud is stable in the region to the left of this line and unstable to the right. We observe three broad instability regimes, approximately delineated by the bold boxes. Even though one can clearly distinguish qualitatively different instabilities, the transition between these regimes is not abrupt, and a precise boundary is, thus, delicate to locate. We stress that although the details of the dynamics in the different regimes may vary with experimental conditions (beam alignment,...), the general features and the separation of phase space in three regions (four with the stable regime) is quite robust and reproducible.

The most striking new regime is observed in the weak ∇B limit. Here the envelope of the cloud remains roughly stationary but with small-scale structures inside undergoing a complex dynamics. As this regime appears to exhibit several small-scale structures, we call it the turbulent regime (TR), noting that a comprehensive study in terms of weak turbulence would require a systematic study of scaling laws, going beyond the scope of this Letter.

Another regime occurs at large ∇B ( 8 G/cm) and close to resonance. There the cloud shows large anisotropic deformations but in random directions. We, thus, call this the statistically isotropic regime (SIR). Everywhere else in the unstable region of the diagram, the cloud experiences deformations with a significant proportion along the MOT beams directions. This is the anisotropic regime (AR). Some examples of single-shot images illustrating the cloud's dynamics in these three regimes and their features described above are provided in Fig. 2.

Although these three regimes of cloud deformation are qualitatively quite different, it seems useful to provide a more quantitative characterization. A well-known technique is principal component analysis [START_REF] Jolliffe | Principal Component Analysis[END_REF], which gives access to a set of uncorrelated spatial patterns involved in the dynamics. In our case, however, the number of such patterns is usually large, and the information is too rich to characterize the different regimes in a simple way. We, thus, employ a cruder method, based on the division of each of the 50 images in a given data set into eight angular sectors as illustrated in Fig. 3. Note that one could, in principle, choose a smaller number of sectors but that would result in a resolution too poor for our purpose. The center of each image is chosen as the center of mass of the average of the 50 realizations. We then compute the histograms of the intensities contained in each sector. The horizontal axis in the histograms shown in Fig. 3 is the intensity in the corresponding sector, normalized to the total intensity in the image (sum of the eight sectors). The vertical axis is the number of occurrences for a given intensity. The total number of occurrences is equal to the number of images (50). This method of analysis provides some information on the magnitude of the cloud's deformation and its anisotropy. Note, however, that this approach does not detect shape-preserving oscillations, whose histograms would not be distinguishable from that of stable clouds.

The data shown in Fig. 3 correspond to the same sets of parameters as in Fig. 2. In (a) (TR), the histograms are similar in all sectors, bell-shaped with approximately the same mean and width, which is rather small (at most 30% of the total in one sector). These observations are consistent with limited deformations of the cloud's envelope, equivalent in all directions. The complex dynamics of the small-scale structures seen in Fig. 2(a) is averaged by the relatively large area of the sectors. In (b) (SIR), we observe that the histograms are also more or less similar in all sectors. However, their shape and width are very different from case (a). It is seen that the histograms correspond to a decreasing probability law with the most probable intensity in each sector being almost zero. This means that the cloud undergoes large deformations associated with large excursions of its center of mass. Concurrently, the large width of the histograms with instances where 90% of the cloud is in one sector, indicates very anisotropic deformations. In case (c) (AR), the situation is again quite different. Here, the shapes and widths are different in all sectors. In particular, we observe double-peaked histograms, such as in sector 2 (arrows). This feature indicates that the cloud jumps repeatedly between two preferred states with different shapes, one of which is highly elongated along the direction of one (or several) MOT beams. Overall, the cloud deformations are quite large with up to 50% of the cloud in one sector. We, thus, see that based on the characteristics of the collection of eight histograms, namely, the dispersion of histograms widths and value of the average width, we can establish a distinction between the different instability regimes. This is confirmed by the color maps shown in (d) with the mean width displayed on the left and the dispersion of the widths on the right. The bold boxes delineating the different regimes, obtained by visual inspection of the sets of 50 images, are the same as in Fig. 1. As can be seen, the result of the eight sectors analysis is in fair agreement with the visual inspection: In the TR both mean and dispersion are small for the IR the mean is large and the dispersion small, and in the AR the mean is intermediate but most importantly the dispersion is large.

IV. NUMERICAL RESULTS

We now investigate the deformation regimes observed with the three-dimensional (3D) simulation scheme, aiming at a direct comparison between our microscopic model and the experiment. As will be seen, the modeling in three spatial dimensions appears essential as it allows for deformation regimes that were not accessible in previous studies (see, e.g., Ref. [START_REF] Pohl | Self-Driven Nonlinear Dynamics in Magneto-Optical Traps[END_REF]).

The numerical simulation proceeds as follows. We start at t = 0 with 7 × 10 3 "superparticles" (see Ref. [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF]) with random positions within a Gaussian density distribution. We then let the cloud evolve under the action of the forces described in Sec. II. In the unstable regime, we observe that the spontaneous cloud oscillations typically occur after a transient time t tr (see Fig. 3 of Ref. [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF]) which varies with the parameters, here δ and ∇B. We then randomly choose 50 time values >t tr and integrate the corresponding 3D density distributions along the z axis to allow for a direct comparison with experimental images. We adopt this approach where we probe different times for a single realization (whereas in the experiment we probe different realizations) because the numerical simulations are too time consuming to sample many realizations. By comparing the results of these different protocols, we implicitly assume ergodicity.

We report in Fig. 4 some examples of instability regimes obtained in the simulations that are reminiscent of those observed in the experiment. In particular, the typical density distributions as well as the histograms obtained with the eight sectors method match fairly well the experimental ones. The first regime (a) also occurring for low magnetic-field gradients is characterized by small-sized structures fluctuating inside a roughly fixed envelope. The structures are not as small and numerous as in the experiment, which we attribute to the coarse graining due to the limited number of superparticles in the simulations. In the second regime (b), observed in a limited range of detunings and for high B-field gradients, the cloud undergoes large deformations in arbitrary directions similar to the SIR of Fig. Thus, in addition to being able to reproduce the scaling laws observed for the instability threshold [START_REF] Gaudesius | Instability threshold in a large balanced magneto-optical trap[END_REF], our model and numerical simulations do predict some instability regimes that are similar to what is observed experimentally. We note, however, that the regions of phase space where these regimes occur are not in perfect match with the experiment. The simulations also predict modes of instabilities that could not be observed in the current experiment (see Fig. 5). The regime shown in Fig. 5(a) occurs for all values of ∇B when the laser frequency is brought close enough to resonance. Typically, the cloud in this regime exhibits two components: A central core which shows spatial and temporal fluctuations and a nearly static "cross" of atoms accumulating on the axes of the laser beams. There is a continuous transition from the anisotropic regime of Fig. 4(c) to this regime when |δ| is decreased. The regime shown in Fig. 5(b) occurs at relatively large ∇B and |δ|. It is characterized by a "twisted" elliptical shape oriented along the axis of one pair of MOT beams. Different simulation runs can yield different orientations. In this regime, the spatial fluctuations of the cloud are reduced. Finally, the regime in Fig. 5(c) occurs at intermediate to large ∇B for detunings just above threshold. There, the cloud assumes an almost elliptical shape, horizontal or vertical. In this regime, the cloud is nearly static, but the trapping symmetry is spontaneously broken. The precise origin of this different behavior in the simulations and experiments remains to be clarified in future numerical work to explore a broader parameter space and scrutinize some of the remaining approximations employed in the simulations.

V. CONCLUSION

We presented in this Letter an experimental study of the spatiotemporal instabilities taking place in a balanced magneto-optical trap containing a large number of atoms. By scanning the laser detuning and magnetic-field gradient at fixed atom number, we identified three broad instability regimes. These observations, thus, constitute a finer description of unstable MOTs using a characterization based on spatial symmetries. We have identified new regimes, beyond previous predictions. These we also observed in extensive three-dimensional numerical simulations. In the most spectacular instability regime observed in the experiments at low ∇B, the cloud breaks into small-scale filaments undergoing complex motions reminiscent of turbulence, e.g., in flames. We plan to investigate this regime in more detail to determine if turbulence is indeed present in this system and of what nature. Making use of the high tunability of the MOT may allow us to induce transitions between different types of turbulence. This ab initio study of optomechanical instabilities contributes to the wider field of collective effects in light-matter coupling, and we speculate that some of the features identified in our experiments and simulations can be present in plasmas and astrophysical systems.

FIG. 1 .

 1 FIG. 1. (δ, ∇B) instability phase diagram (experiment). The atom number is fixed at N = 1.5 × 10 10 . Each image corresponds to the average of 50 single-shot fluorescence images with intensity displayed in logarithmic scale (the intensity range is 2 decades). The field of view is 6.5 × 6.5 cm 2 . The dashed line materializes the instability threshold. The three different instability regimes (see the text) are roughly delineated by the bold boxes.

  FIG. 2. Different instability regimes (experiment). We show single-shot images illustrating the cloud's dynamics in the three regimes identified in Fig. 1. (a) TR (∇B = 1.2 G/cm, δ = -). (b) SIR (∇B = 12 G/cm, δ = -). (c) AR (∇B = 7.2 G/cm, δ = -1.8 ). The dotted lines in the bottom right image show the directions of two pairs of MOT beams.

FIG. 4 .

 4 FIG. 4. Different instability regimes (simulations). We illustrate three of the instability regimes obtained in the simulations that are similar to those observed in the experiment (see Fig. 2). The images in each line correspond to different times in the simulation. (a) ∇B = 1.2 G/cm, δ = -1.3 . (b) ∇B = 12 G/cm, δ = -2.6 . (c) ∇B = 4.8 G/cm, δ = -1.9 .

2 .

 2 In the third regime (c), the cloud's deformations are mostly occurring along the directions of the MOT's laser beams (the diagonals of the images) with the cloud jumping from side to side in a manner very similar to what is experimentally observed in the AR [Figs. 2(c), 3(c) and 4(c)]. This is an important observation, which shows that this kind of spontaneous symmetry breaking leading to large center of mass excursions can occur in a perfectly centrosymmetric trapping potential in the absence of experimental imperfections, such as local intensity imbalance or beam misalignment.

FIG. 5 .

 5 FIG. 5. Other instability regimes (simulations). These instability regimes are not observed in the experiment. (a) ∇B = 7.2 G/cm, δ = -1.1 . (b) ∇B = 9.6 G/cm, δ = -2.9 . (c) ∇B = 6 G/cm, δ = -3.3 .

ACKNOWLEDGMENTS

This work was performed in the framework of the European Training Network ColOpt, which is funded by the European Union (EU) Horizon 2020 Program under the Marie Skłodowska-Curie action Grant Agreement No. 721465, and in the framework of the European Project ANDLICA, ERC Advanced Grant No. 832219. It has also been supported by the Danish National Research Foundation through a Niels Bohr Professorship to T.P. and through the Center of Excellence "CCQ" (Grant Agreement No. DNRF156).


Patrick Gaulme
email: gaulme@mps.mpg.de

François-Xavier Schmider

Ivan Gonçalves

Measuring Planetary Atmospheric Dynamics with Doppler Spectroscopy

Keywords: Methods, observational -Techniques, imaging spectroscopy -Techniques, radial velocities -Planets and satellites, atmospheres -Planets and satellites, individual, Venus -Planets and satellites, individual, Jupiter

Doppler imaging spectroscopy is the most reliable way to directly measure wind speeds of planetary atmospheres of the Solar system. However, most knowledge about atmospheric dynamics has been obtained with cloud-tracking technique, which consists of tracking visible features from images taken at different dates. Doppler imaging is as challenging -motions can be less than 100 m s -1 -as appealing because it measures the speed of cloud particles instead of large cloud structures. Significant difference is expected in case of atmospheric waves interfering with cloud structures. The purpose of this paper is to provide a theoretical basis for conducting accurate Doppler measurements of planetary atmospheres, especially from the ground with reflected solar absorption lines. We focus on three aspects which lead to significant biases. Firstly, we fully review the Young effect, which is an artificial radial velocity field caused by the solar rotation that mimics a retrograde planetary rotation. Secondly, we extensively study the impact of atmospheric seeing and show that it modifies the apparent location of the planet in the sky whenever the planet is not observed at full phase (opposition). Besides, the seeing convolves regions of variable radial velocity and photometry, which biases radial-velocity measurements, by reducing the apparent amplitude of atmospheric motions. Finally, we propose a method to interpret data, i.e., how to retrieve zonal, meridional, vertical, and subsolar-to-antisolar circulation from radial velocity maps, by optimizing the signal to noise ratio.

Introduction

For dense atmospheres in the solar system, wind measurements mostly come from the cloud tracking technique. The method consists of following cloud features at specific wavelengths taken on image pairs obtained at different times, to retrieve the wind speed before cloud structures evolve or disappear (Sánchez-Lavega et al. 2008;Peralta et al. 2008;Moissl et al. 2009;Peralta et al. 2012). This method has also provided results with spacecrafts that flew-by or orbited planets (e.g. Choi et al. 2007). Cloud tracking indicates the motion of large cloud structures (limited by spatial resolution), which is rather an indication of the speed of iso-pressure regions, than the speed of the actual cloud particles. A complementary solution to access direct wind speed measurement is Doppler spectrometry, because it measures the actual speed of cloud particles. So far, most observational efforts with Doppler measurements of dense planetary atmospheres have been dedicated to Venus and Jupiter, for two distinct reasons. Doppler spectrometry was envisioned for Venus to understand its atmospheric circulation as it is a mostly featureless planet in the visible domain. First attempts based on visible spectrometry were performed in the 1970s, but did not lead to robust results (Traub & Carleton 1975;Young et al. 1979). In 2007, a significant international effort was organized to support the atmospheric observations of Venus by ESA mission Venus Express (VEx) (Lellouch & Witasse 2008). The objective was to mea-sure the atmospheric circulation using different spectral ranges, to probe different altitudes in the Venus mesosphere. Significant results on the upper mesospheric dynamics were obtained using mid-infrared heterodyne spectroscopy (Sornig et al. 2008(Sornig et al. , 2012) ) and millimeter and submillimeter wave spectroscopy (Clancy et al. 2008(Clancy et al. , 2012;;Lellouch et al. 2008;Moullet et al. 2012) but Doppler spectroscopy is more challenging at shorter wavelength. Visible observations of solar Fraunhofer lines scattered by Venus clouds were performed by Widemann et al. (2007Widemann et al. ( , 2008)); Gabsi et al. (2008); Gaulme et al. (2008b); Machado et al. (2012Machado et al. ( , 2014Machado et al. ( , 2017)). Regarding other planets, similar measurements were performed in the mm/sub-mm domain for Mars (Lellouch et al. 1991;Moreno et al. 2009) and Titan (Moreno et al. 2005), as well as in the infrared with the 10-µm heterodyne observations of Titan by Kostiuk et al. (2001Kostiuk et al. ( , 2005Kostiuk et al. ( , 2006Kostiuk et al. ( , 2010)).

The use of Doppler spectrometry for giant planets, Jupiter in particular, was inspired by its success with helioseismology to probe their deep internal structures (e.g., Appourchaux & Grundahl 2013). The principle relies on the detection of global oscillation modes, whose properties are function of the internal density profile. Giant planets being mostly fluid and convective, their seismology is much closer to that of solar-like stars than that of terrestrial planets. The basic principle relies on monitoring the position of a spectral line that probes an atmospheric level where the amplitude of acoustic modes is maximum. So far, resonant cell (Schmider et al. 1991;[START_REF] Cacciani | [END_REF] and Fourier transform spectrometry (FTS, e.g., Mosser et al. 1993;Schmider et al. 2007;[START_REF] Gonçalves | Proc. SPIE[END_REF]) have been considered. Two dedicated instruments were subsequently designed: the SYMPA instrument was an FTS based on a Mach-Zehnder interferometer (Schmider et al. 2007) which has provided so far the clearest observational evidence of Jupiter's oscillations (Gaulme et al. 2011). Since then, the new instrument JOVIAL/JIVE, inherited from SYMPA, has been developed and tested to perform both atmospheric dynamics of dense atmospheres and seismic measurements of Jupiter and Saturn (Gonçalves et al. 2016, and submitted).

In this paper, we first study the influence of two possibles sources of biases regarding Doppler spectroscopic measurements of planetary atmospheres. On one hand, we review and extensively revise the impact of the solar rotation on radial velocity measurements, which was originally studied by Young (1975) in the case of Venus (Sect. 2). On the other hand, we consider the impact of atmospheric seeing on both positioning on the planet and on radial velocity maps (Sect. 3). In a third part, we focus on how to retrieve zonal and meridional components of atmospheric circulation from radial velocity maps (Sect. 4).

Effect of the apparent diameter of the Sun

2.1. The Young effect Young (1975) was the first to point out that the apparent diameter of the Sun as seen from Venus causes an artificial radial velocity field on it. Indeed, rays from different parts of the Sun, which display different radial velocities, reach the planet with (slightly) different incidence angles. Regions of the Sun that are closer to Venus' horizon contribute less to the reflected solar spectrum than regions closer to zenith. Thus, the radial velocity integrated over the whole solar disk is not zero in a given point of Venus. In other words, even if Venus would not be rotating, we would still measure a Doppler shift near the Venus terminator, and that Doppler shift would mimic a retrograde rotation because the solar rotation is prograde.

Despite its pioneering aspect, Young (1975)'s analytical description of the artificial Doppler shift ∆V Y on Venus was limited to the equator. Out of the equator, the expression ∆V Y ≈ 3.2 tan γ, where γ is the solar-zenith angle, is no more valid. Nonetheless, most authors who hitherto worked with Doppler spectroscopic data of Venus in the visible used Young's formulation, without noticing it was valid at the equator only, and certainly got biases in their results and conclusions (Widemann et al. 2007(Widemann et al. , 2008;;Gabsi et al. 2008;Machado et al. 2012Machado et al. , 2014Machado et al. , 2017)). We propose a revised analytical description of the "Young" effect with the same approximations used by Young (1975), as well as a more comprehensive numerical simulation.

Mathematical description

In this section, we focus on Venus, as this effect is almost negligible for external planets, even though our analysis is valid for any planet of the solar system. We do not select Mercury as the test case, as such a kind of measurements is useful only to monitor dynamics of optically thick atmospheres and not solid surfaces covered by fixed contrasted features. To express the artificial Doppler shift ∆V Y , known as Young effect, we assume that the uppermost cloud layer of the considered planet is a flat surface. Within this assumption, at a given point on Venus, the intensity of a ray coming from the Sun is proportional to the cosine of the incidence angle γ.

γ θ O M S A θ = -π/2 θ = -π/4 θ = 0 θ = π/2 θ = π/4 θ = -3π/4 θ = 3π/4
Fig. 1. Reference frame used to express the Young effect. The points S and A are the sub-solar and anti-solar points respectively, O is the center of the coordinates and M a given point on the planet. In this frame, the angles of incidence γ and horizon inclination θ are uniform along parallels and meridians respectively. The sub-pictures represent the radial velocity map of the Sun a seen with θ = [0, -π/4, -π/2]. The artificial radial velocity is not observed along the meridian θ = 0, while it is maximum along θ = ±π/2. Two factors drive the Young effect, which were both mentioned by Young himself but only the first one was implemented. Firstly, the variation of incidence angle from different parts of the Sun results in averaging regions with different radial velocities and intensities. Besides, the center part of the Sun appears brighter than the edge -limb darkening effect -which slightly counterbalances the artificial Doppler signal, by giving more importance to regions with low radial velocities. Secondly, the inclination of the solar rotation axis with respect to the horizon modulates the Young effect. For example, near Venus poles, where the solar rotation axis is perpendicular to the horizon, both the blue and red sides of the Sun reach the surface in a symmetrical fashion and average out. To the contrary on Venus morning terminator, the blue side of Venus rises first and the artificial Doppler signal is maximum (Fig. 1). This second factor was not taken into account in Young (1975)'s equation, because he focused on describing the amplitude of the phenomenon along the equator, where it is maximum.

In other words, the parameters that drive the Young effect are: 1) the angle of solar incidence γ (solar-zenith angle); 2) the inclination θ of the solar rotation axis with respect to the local horizon. Let us consider the coordinate system that connects the sub-solar to the anti-solar points of Venus (Fig. 1). The incidence angle γ is uniform along parallels of this system of coordinates, while the tilt of the solar spin axis with respect to the horizon θ is uniform along the meridians. In a given point on Venus of solar incidence γ and rotation axis inclination θ, the artificial Doppler effect is the integral of the radial velocity map weighted by the differential photometry:

∆V Y (γ, θ) = ρ -ρ √ ρ 2 -y 2 - √ ρ 2 -y 2 I LD (x, y) cos(γ -y)∆V (x, y, θ) dx dy ρ -ρ √ ρ 2 -y 2 - √ ρ 2 -y 2 I LD (x, y) cos(γ -y) dx dy (1)
where ρ is the solar apparent diameter as seen from Venus, I LD (x, y) its limb-darkening law, and ∆V (x, y, θ) is the radial velocity field of the solar surface, expressed in the frame (x, y) centered on the solar disk (Fig. 2). 2. Schematic view of the photometric aspect of the Sun as seen from Venus. The angle γ is the solar-zenith angle, i.e., the angle between the zenith and the center of the Sun. The system of coordinates (x, y) is centered on the Sun; the y-axis points to the zenith. Hence, the ordinate of the solar region closest to the horizon is γ-y. All variables are angles.

γ 0 y 2 -y 2 2 -y 2 - ρ ¤ ρ ¤ ρ ¤ Fig.

Approximated analytical solution

In this section, we demonstrate the analytical expression of the Young effect, by assuming the same approximations that were used by Young (1975). The first approximation consists of considering the Sun as a solid rotator, i.e., with uniform rotation period as function of latitude, which is actually quite far from the actual Sun, for which differential is significant (24 days at equator and 38 in polar regions, according to Snodgrass & Ulrich 1990). The second approximation consists of neglecting any solar limb-darkening, which is not accurate either (e.g., [START_REF] Allen | Astrophysical quantities Appourchaux[END_REF]. The third approximation consists of considering Venus at top-cloud level as a Lambert scattering surface, i.e. iso-scattering by surface element. Notations for angles are described in Fig. 3.

The radial velocity map of a solid rotator as seen from its equatorial plane is a linear function of the coordinates on the considered rotator. Let us consider the (x , y ) system centered on the Sun and tilted accordingly to its rotation axis (Fig. 3), the radial velocity field ∆V (x , y ) of the solar surface is:

∆V (x , y ) = V eq ρ x , (2)
where V eq si the solar equatorial velocity, and x , y , and ρ are expressed in arcsec. By considering that the Sun is tilted by the angle θ, it re-expresses in the local system (x, y), in which x is parallel to the horizon:

∆V (x, y, θ) = V eq ρ (x cos θ + y sin θ) . (3)
Since we neglect the solar limb-darkening, Eq. 1 gets simplified:

∆V Y (γ, θ) = V eq ρ ρ -ρ √ ρ 2 -y 2 - √ ρ 2 -y 2 cos(γ -y) (x cos θ + y sin θ) dx dy ρ -ρ √ ρ 2 -y 2 - √ ρ 2 -y 2 cos(γ -y) dx dy (4) = V eq ρ sin θ ρ -ρ y ρ 2 -y 2 cos(γ -y) dy ρ -ρ
ρ 2 -y 2 cos(γ -y) dy

(5) The solar rotation axis is tilted by the angle θ. The system of coordinates (x, y) is the frame in which we compute the integral in Eq. 1, while (x , y ) is the system in which the solar rotation is described.

By assuming that y γ, which is true except for the subsolar area, we can develop cos(γ -y) into cos γ + y sin γ and analytically solve the integral:

∆V Y (γ, θ) = V eq ρ sin θ sin γ cos γ ρ -ρ y 2 ρ 2 -y 2 dy ρ -ρ ρ 2 -y 2 dy (6)
The integration of the latter expression leads to a revised version of the Young effect:

∆V Y (γ, θ) = V eq ρ 4 sin θ tan γ (7)
By assuming a solar equatorial velocity of 2 km s -1 , and by considering that the apparent solar radius as seen from Venus is 6.43 10 -3 rad, V eq ρ /4 ≈ 3.2 m s -1 . Our expression differs from Young's only by the introduction of the sin θ dependence. The numerical simulation that is described in the next section allows us to adjust the amplitude factor -here 3.2 m s -1 -when we take into account the solar limb darkening and differential rotation.

As pointed out by previous authors (e.g., Machado et al. 2012), the expression of the Young effect diverges at the terminator, where γ = π/2. This is because the above expression is a good approximation as long as the Sun is entirely visible from a given point on Venus, i.e., γ < π/2-ρ . Once the Sun is partially below the horizon, Eq. 7 is no more valid, and Eq. 5 must be integrated from the horizon line (ordinate y = γπ/2) to ρ (see Fig. 2). In that case, the Young effect near terminator becomes:

∆V Y (γ, θ) = V eq ρ 4 sin θ Φ sin γ + X cos γ ρ 2 4 X sin γ + Ψ cos γ (8)
where

Φ(γ) = π 2 -sin -1 γ -π 2 ρ - γ -π 2 ρ 2 ρ 2 -γ - π 2 2 1 2 + 2 γ -π 2 ρ 4 ρ 2 -γ - π 2 2 3 2 (9) X(γ) = 8 3 1 ρ 4 ρ 2 -γ - π 2 2 3 2
(10)

A&A proofs: manuscript no. Gaulme_Planet_Doppler_Spectro ). The plain red line corresponds with the case of the rising Sun (Eq. 8). The vertical line indicates the first contact of the Sun with the horizon. The minimum value of ∆V Y is -2, 000 m s -1 corresponds with the Sun's equatorial velocity V eq used for the simulation.

Ψ(γ) = π 2 -sin -1 γ -π 2 ρ - γ -π 2 ρ 2 ρ 2 -γ - π 2 2 1 2 (11)
which is valid as long as γ ∈ π 2 -ρ , π 2 + ρ . Note that we include the expression of the Young effect in this configuration for completeness only. Actually, it does not make much sense from a physical point of view since other effects interfere, as the atmospheric refraction. In addition, as we indicate in Sect. 3.3 (Fig. 12), the Young effect is drastically diminished by the bias on radial velocities caused by instrumental and/or atmospheric Point Spread Function (PSF).

Numerical simulations

It is impossible to analytically integrate Eq. (1) without significant approximations, such as considering y γ. We opt for a numerical simulation based on the observational knowledge of the Sun's limb darkening and differential rotation. The simulations are built into two steps. First, we make a Venus 3dimension model that we project into a 2-dimension frame, on which the angles γ and θ are known for each point of Venus surface (see Appendix A). Secondly, for each point on Venus, we simulate the rotating Sun and integrate Eq. (1).

From simple trigonometric considerations, the two angles γ and θ can be expressed in the planet's coordinate system:

γ = cos -1 cos(λ -λ S ) cos(φ -φ S ) (12) θ = tan -1 - sin(φ -φ S ) tan(λ -λ S ) (13)
where λ and φ are the latitude and longitude of the planet, and λ S and φ S are the coordinates of the sub-solar point in this frame. Latitudes run from 0 at the equator to π/2 at the poles in both hemispheres, and longitudes run from -π/2 at the evening terminator to π/2 at the morning terminator.

-1000 In each position on Venus we consider both the solar photometric map, i.e. a disk dominated by the limb darknening law weighted by cos(γ -y), and the radial velocity map, which includes the solar differential rotation and the inclination θ with respect to the local horizon. According to [START_REF] Allen | Astrophysical quantities Appourchaux[END_REF], the limb darkening law can be described as a quadratic law:

-1 0 0 -1 0 0 -2 5 -2 5 -10 -1 0 -1 0 -10 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 10
I(µ) = I(0) 1 -u 2 (1 -µ) -v 2 (1 -µ 2 ) , (14)
where µ = cos γ , and where γ is the angle between the normal to the solar surface and the direction pointing at the observer. A point located at a radius r from the solar center is function of the γ angle as r = R cos γ . Hence, µ = cos sin -1 (r/R) , i.e., µ = 1 -(r/R) 2 . Therefore, Eq. ( 14) can be expressed in the (x, y) coordinate system:

I LD (x, y) = 1 -u 2        1 -1 - x 2 + y 2 R 2        -v 2 x 2 + y 2 R 2 (15)
where we set I(0) = 1. From [START_REF] Allen | Astrophysical quantities Appourchaux[END_REF], the solar limb darkening coefficients are u 2 and v 2 are given in a table that covers the whole visible range.

According to Snodgrass & Ulrich (1990), the angular velocity of the solar differential rotation can be described with the function:

Ω = A + B sin 2 λ + C sin 4 λ ( 16
)
where λ is the latitude on the Sun, and A = 2.97 ± 0.01, B = -0.48±0.04, and B = -0.36±0.05 µrad s -1 . The rotation velocity map is then R Ω cos λ , and the radial velocity map ∆V is the projection of it towards Venus:

∆V = R A + B sin 2 λ + C sin 4 λ cos λ cos φ , ( 17
)
where R is the solar radius in meters, and φ the solar longitude as seen from Venus, ranging from 0 to π from West to East. The numerical integration of Eq. ( 1) is therefore the sum of the solar radial velocity map ∆V times the photometric contribution I LD cos(γ -y), then normalized by the sum of the photometric term. We simulated the Sun with a 256×256 element grid, which is enough in terms of precision. Note that next to the terminator, we take into account that the Sun is partially below the horizon.

The maximum value of the Young effect corresponds to the configuration where only a thin portion of solar equator is visible, i.e., about 2 km s -1 . We represent the Young effect both for a 2Dmap with θ ranging from -180 to 180 • and γ from 0 to 90 • (Fig. 5), and also for a example of Venus figure as seen from the Earth. Simulations were done with a limb darkening law corresponding to 550 nm optical wavelength. The two figures highlight the fact that this effect is significant only near the terminator, as it is less than ±10 m s -1 at incidence lower than ≈ 75 • , whereas zonal winds are expected to exceed 100 m s -1 .

To determine the empirical coefficient of the Young effect, we computed the ratio of the numerical Young-effect map to tan γ sin θ. We also considered limb-darkening laws that covers the whole visible range (Fig. 7). The coefficient appears to range from 2.7 to 3 from 350 to 800 nm, which is logically smaller than the 3.2 value, obtained by neglecting the Sun's limb-darkening.

We suggest to use the following expression for the Young effect:

∆V Y (γ, θ) = Y(Λ) sin θ tan γ ( 18
)
where Λ indicates the optical wavelength, and the coefficient Y(Λ) is indicated in Table 1 for the major planets of the Solar system. Even it is useless for measuring atmospheric dynamics, we include Mercury and Mars in Table 1 for completeness and also because detecting the Young effect in a solid-rotating body would be a good test for our ability to lead high precise radial velocity measurements. Again, near the terminator, the expression of the Young effect must take into account that the Sun is partially below the horizon:

∆V Y (γ, θ) = Y(Λ) sin θ Φ sin γ + X cos γ ρ 2 4 X sin γ + Ψ cos γ , (19)
with the same notations as in Eq. 8.

Impact of atmospheric seeing

Biases caused by the point spread function

Instrumental PSF is another major cause of biases in planetary observations. Since we mainly consider earth-based observations in the visible, the PSF is dominated by the atmospheric seeing, so we hereafter refer to the effects of atmospheric seeing instead of instrumental PSF. However, the effects that we describe are also valid for any PSF, even in the case of space-borne measurements, although the biases we identify are much smaller. Firstly, seeing convolves the planetary shape and photometry by the PSF, which is usually described as a Gaussian. For external planets at their opposition, the seeing simply blurs the planetary image, which is not a major issue, except for defining the edge. Indeed, it is not alway possible to perfectly know the pixel field-of-view (FOV) of an instrument mounted at a telescope focus, and convert the ephemeris apparent diameter from arcsec to pixels 1 . For planets that are out of opposition, especially the inner planets Mercury and Venus, important phase angles (Sunplanet-observer) lead to significant mislocation of the planet because the blurring effect caused by the atmospheric seeing is not the same at limb or terminator edges. Secondly, besides positioning issues, atmospheric seeing blends regions with different photometry and radial velocity, which biases radial velocity measurements, more importantly towards the planetary edges, where photometric variations are steep. We insist on the term bias because atmospheric seeing is often seen as a phenomenon that smooths out radial velocities, while it actually introduces systematic biases in the measurements, which lead to wrong estimates of the wind circulation even with excellent SNR data (see Sect. 3.3, and Fig. 11).

Positioning on the planet

A fundamental step in measuring winds of a planet with a mostly featureless dense atmosphere is to make sure to know where the telescope is slewing. If ephemeris are very accurate for the main solar system bodies, pointing errors of any telescope are typically of about 10 arcsec, which is not enough for an extended object of that kind of apparent size. The use of a guider camera then allows the observer to position the spectrometer's pinhole/slit on the planetary disk. However, this step is not as accurate as one could think.

There are commonly two ways of determining the location on an extended object as a planet. A first approach is based on the calculation of the photometric centroid of the image, by assuming the background to be negligible. This statement is rather uncertain because the background signal is generally not uniform and scattering is often present around any bright object. This could bias the determination of the photocenter. Moreover, some other factors may bias the photocenter with time. For instance, in the case of Jupiter, visible atmospheric features, moon transits and their cast shadows alter the measured photocenter. The second approach consists of identifying the planetary edge on the image. If the pixel FOV is perfectly known, the observer applies a photometric threshold on the actual image until the portion of the planet above the threshold matches the expected size of the planet on the detector. If the pixel FOV is not perfectly known, a photometric threshold on the image refines it. For example an observer can assume that the planet is all what is above 5 % of the maximum intensity, and adjust the pixel FOV accordingly. Since this is rather arbitrary, we assume that the pixel size can be known (it can be calibrated using double stars or motion applied to the telescope). In the following, we only consider the second approach. We show that the size of the PSF leads to misestimating the planet's position on the detector, and we determine what correction must be applied. Note that the photometric centroid does not suffer from such a bias related to atmospheric seeing, but as said earlier, it could be affected by other uncalibrated variations.

We simulate a Lambert sphere seen at different phase angles and degrade it with typical seeing conditions ranging from 0.5 to 4 arcsec. A four-arcsec seeing can sound quite large but is not uncommon for daytime observations, as it often happens for Venus (e.g., Gaulme et al., in prep.). In Fig. 8, we present the cut of a Lambert sphere along the equator as well as the same observed with atmospheric seeing of 1 and 3 arcsec. Phase angles range from 0 to 135 • . At zero phase angle, atmospheric seeing impacts only the apparent diameter of the planet because the blurring effect is symmetrical around the edge. Then, as the phase angle increases, the planet seems to get shifted towards the Sun's direction. This is due to the fact that the limb edge is steeper than the terminator, making the seeing-caused blurring stronger on limb side than on the terminator. In Fig. 8, we also represent two examples of images, without and with seeing, which show that the planet seems to drift away from its actual position.

To quantify how much seeing alters the apparent position of a planet, we simulate a Lambert sphere in the conditions Venus and Jupiter are observed. For this, we download the ephemeris of the two planets from the NASA Horizons web interface (https://ssd.jpl.nasa.gov/horizons.cgi) on long time periods to get an average relationship connecting the phase angle and apparent diameter (Fig. 9). For Venus, we choose to present the results as function of the phase angle as it is a critical criterion when observing it. Indeed, winds are impossible to measure when a small fraction of the disk is visible, and very low or very large phase angle correspond to bad observing conditions near Patrick Gaulme et al.: Measuring Planetary Atmospheric Dynamics with Doppler Spectroscopy solar conjunctions. As regards Jupiter, for a given phase angle, two possible apparent diameters are possible as it is an outer planet. This causes issues to produce a clear plot, as two results should be presented at each value of the phase. This is why we choose to present Jupiter simulations as function of apparent diameter.

We consider positive phase angles, which means the Sun is on the left side of the simulated planet. As indicated in the introduction of this section, we assume the pixel FOV and planet diameter to be known and determine the location of the planet from the photometric profile along the equator, as done in Fig. 8. In Fig. 10, we display the apparent shift of "Lambertian" Venus and Jupiter along the x-axis from the position of the equator, expressed as a fraction of the planet's apparent radius. The common characteristic to both planets is that a zero phase angle leads to no apparent displacement of the planet, because the seeing symmetrically alters both sides.

As regards Venus, Fig. 9 firstly reminds us how small Venus is relatively to Jupiter if observations aim at measuring winds, i.e. if more than half of the dayside is visible: its apparent diameter is less than 25 arcsec. Therefore, the seeing impacts it more than Jupiter. Then, we notice than beyond about 55 • , the apparent displacement of the planet gets smaller to reach about zero at 180 • . This is due to the fact that photometric profile near terminator is becoming steeper, as can be seen in Fig. 8 for α = 90 and 135 • . Next to α = 180 • the dayside crescent is so thin that seeing degradation is practically symmetrical on limb and terminator edges, making the apparent displacement to be none. Overall, for Venus, the apparent shift of the planet is large at it reaches from 4 and 17 % of its radius within seeing conditions ranging from 0.5 to 4 arcsec. For a typical observation with 2arcsec seeing and phase angle of 70 degrees the planet appears shifted of about 11 % of its radius.

As regards Jupiter, the maximum impact of seeing on positioning occurs at Jupiter's maximum phase angle, which is about 11 • for an apparent diameter of 38 arcsec. In such a case the apparent shift of the planet ranges from 0.8 to 3.6 %, which sounds almost negligible. However, Jupiter's equatorial velocity is about 12.5 km s -1 . Remembering that the Doppler signal of a reflected solar line is enhanced by a factor (1 + cos α), i.e. almost 2 for α ≤ 11 • , the Doppler signal varies of about 50 km s -1 along Jupiter's equator. A shift of 3 % along the equator is equivalent to a 1.5-km s -1 radial velocity bias. It means that seeing impact on positioning needs to be considered too for a planet as Jupiter. We further develop some observation strategy to circumvent this positioning bias in Sect. 3.5.

Radial velocities

As originally pointed out by Civeit et al. (2005), atmospheric seeing alters radial velocity measurements as it blends regions with non-uniform radial velocity and photometry. The resulting radial velocity is thus the convolution of the radial velocity signal with the photometric map of the considered object, including its degradation by seeing. The mean Doppler ∆V measured in a given pixel (x, y) on the detector can be expressed as:

∆V (x, y) = (∆V F * P)(x, y) (F * P)(x, y) (20)
= ∆V(u, w) F(u, w) P(xu, yw) du dw

F(u, w) P(x -u, y -w) du dw (21)
where ∆V and F are the "real" radial-velocity and photometric maps prior to seeing degradation, P is the PSF of the atmo- For both panels, the plain red curves indicate moving averages of the black dots, which were then used to associate an (two for Jupiter) angular diameter to a given phase angle in Fig. 10. spheric seeing, and the asterisk sign * indicates the convolution product. This effect produces an artificial Doppler shift in opposition to the actual planetary rotation radial velocity. Indeed, because of the seeing, regions of larger photometric intensity and low Doppler signal contribute more to the Doppler map than the faint edges with high Doppler signal. Civeit et al. (2005) originally studied this question for attempting to measure the rotation of Io, one of Jupiter's satellites. It was then used to measure the rotation of Saturn's moon Titan (Luz et al. 2005(Luz et al. , 2006)). In both cases, the objects apparent sizes are smaller than the seeing PSF, which motivated the authors to take the seeing bias on radial velocities into account. The only visible Doppler-spectroscopic based observations of other planets regarded Venus, and none of the authors considered the seeing bias (Widemann et al. 2007(Widemann et al. , 2008;;Gabsi et al. 2008;Gaulme et al. 2008b;Machado et al. 2012Machado et al. , 2014Machado et al. , 2017)). The reason is likely that since the apparent diameter is significantly larger than the seeing PSF, most authors considered it to be negligible. In this section, we first show that such an assump- tion is erroneous, and second, that the Young effect turns out to be much altered. For a zonal circulation, the atmospheric seeing reduces the amplitude of radial velocities because the regions where the Doppler shift is maximum display lower photometry. In Fig. 11, we present simulations of a Lambert sphere that rotates as a solid-body, with 1-m s -1 amplitude at equator, 15-arcsec apparent diameter, and phase angles (0, 45, 90, 135 • ) as in Fig. 8. We assume radial velocities to be measured on solar Fraunhofer lines, so the Doppler shift is enhanced by (1 + cos α). Radial velocities along the equator are compared to what is measured by assuming a 1-and 3-arcsec seeing. At zero phase angle, a 1arcsec seeing reduces the measured radial velocity of about 6 % at the edge, whereas a 3-arcsec seeing reduces it of about 18 %. The larger is the phase angle, the larger is the bias on radial velocities.

As regards the Young effect, since its amplitude varies very steeply near the terminator, it is obviously much altered by at- mospheric seeing. In Fig. 12, we present a map of the Young effect in the case of Venus, in the same conditions as in Fig. 6 (α = 40 • , D app = 11.6 arcsec), as it appears with 2-arcsec atmospheric seeing conditions. As expected, where photometric variations are slow, the same order of magnitude is preserved, while towards the morning terminator, the maximum amplitude is reduced to 17 m s -1 instead of about 2,000 m s -1 . This tells that if the Young effect should not be neglected when searching for winds of about 100 m s -1 , as its amplitude can be over 10 m s -1 in the early morning area, it also never reaches values as large as the solar rotation velocity, even with infinitely small pixels.

Combined effect of seeing photometry and radial velocities

In the previous two sections, we show that atmospheric turbulence: 1) generates an apparent shift of a planet seen with a nonzero phase angle; 2) biases radial velocities, especially close to the edges, where the photometric gradient is steep 2 . In this paragraph, we comment whether one or both of these effects lead to over-or underestimating the amplitude of a global atmospheric circulation pattern, especially in the zonal case.

The effect of a PSF on estimating a wind circulation is nonlinear and impossible to fully predict. Nonetheless, we identify different possible scenarios that depend on the observing configuration. us consider of an inner planet, as Venus, where the atmospheric circulation is purely zonal. Let us also neglect the Young effect, given it is almost canceled by the PSF. In the case of a correct positioning of the spectrometer on the planet -apparent shift is taken into account -the reduction of the radial velocity field towards the edges leads to systematically underestimating the zonal wind (Fig. 13).

To the contrary, when the positioning on the planet is erroneous -apparent shift is not taken into account -, the effect of the seeing can lead to either case. We remind that the apparent displacement of the planet caused by atmospheric seeing is towards 2 Note that in the absence of a photometric gradient, the seeing would simply smooth the radial velocity map. Because of spatial photometric variability, such as limb darkening, the atmospheric seeing introduces a systematic bias in radial velocity measurements of any extended object. the Sun, wherever the phase angle is positive or negative (Fig. 13). If the spectroscopic measurements are away from the limb and terminator regions, as it is often the case for the Venus observations led by Widemann and Machado, the apparent displacement of the planet leads to overestimating the wind amplitude because the radial velocity profile gets monotonically steeper towards the limb. In the case spectroscopic measurements cover the whole planet, the wind estimation is the result of the competing phenomena: underestimation caused by reduced radial velocities close to the edges, and overestimation by apparent shift. The outcome is unpredictable a priori. In any case, at some point, if the atmospheric seeing is very large -enough to alter the whole radial-velocity profile, as with the 3-arcsec case in Fig. 11 -underestimating the wind amplitude dominates. Interestingly, the above comments could explain in part why zonal-wind measurements of Venus by Machado et al. (2012Machado et al. ( , 2014Machado et al. ( , 2017) ) provide very consistent results -both internally and with respect to cloud tracking -, with average zonal winds estimated in between 117 and 123 m s -1 , whereas Widemann et al. (2007); Gabsi et al. (2008) report much lower values (≈ 75 m s -1 ). The former results were obtained either with the CFHT or the VLT telescopes which are located on sites where the seeing is very good (maximum seeing reported is 1.4 arcsec in Machado et al. 2012). To the contrary, both papers Widemann et al. (2007) and Gabsi et al. (2008) report bad observing conditions at the Observatoire de Haute Provence, which could explain why they obtain such low wind estimates.

How to deal with such kind of data

Except for space-borne observations, atmospheric seeing alters radial-velocity measurements. Adaptive optics could considerably reduce the seeing effect, but it has never been used so far for this type of observations, because adaptive optics systems are relatively rare on high-resolution spectrometric devices.

Adaptive optics is planned as part of the second stage of the JOVIAL/JIVE project [START_REF] Gonçalves | Proc. SPIE[END_REF]. In this section, we establish a data analysis protocol to optimize the reliability of such kind of observations. The first step is to know where "we are" on the planet. Having a partial or full image of the planet in parallel to spectroscopic data is fundamental. It could be a guiding camera image as in Widemann-or Machado-led papers, or a slice of planet as in Gaulme et al. (2008b, and in prep). The image is the result of the apparent diameter of the planet, the telescope optics and the atmospheric turbulence. Let suppose the photometric limbdarkening profile of the planet to be known, one could retrieve an estimate of both the atmospheric seeing and the pixel FOV, by fitting a photometric profile degraded by a PSF. This is the approach used by Gaulme et al. (in prep.) for Venus, as the observations were taken mostly during daytime so there was no way to monitor the atmospheric seeing. Different levels of refinements are possible for this approach. Firstly, the pixel FOV can be calibrated by observing apparently close stars (e.g., wide multiple-stars systems), which reduces the number of free parameters. Secondly, the limb darkening law of the planet can be very simple (e.g., Lambertian), quite refined by including improved description of light scattering, or even based on spaceborne data (e.g., Mayorga et al. 2016). From Gaulme et al. (in prep.) and Gonçalves et al. (submitted) about Venus and Jupiter respectively, a Lambert is usually enough. Indeed, Venus' smooth in mid-visible and small apparent diameter (5 to 10 times a typical seeing value) does not justify models that are too complex to be distinguished from a Lambert sphere given the observation quality. As for Jupiter, its apparent diameter is usually much larger than the seeing, and it is observed at low phase angle, where the seeing effect is radially symmetrical. The ideal situation involves the presence of a Differential Image Motion Monitor (DIMM, Vernin & Munoz-Tunon 1995) at the same observatory, which points at the same location on the sky (atmospheric seeing can substantially vary according to azimuth and altitude).

In Fig. 14, we show a photometric profile of Venus obtained at the Teide observatory together with a fit of it by a Lambert limb-darkening law. This profile is part of a larger dataset which will be presented in Gaulme et al. (in prep.). Fitting was performed with a by a dedicated Monte-Carlo and Markov Chain (MCMC) routine, based on the Metropolis-Hasting algorithm and parallel tempering, and Bayesian inference. Free parameters are the apparent diameter R, the seeing, as well as the maximum photometric flux H and the position δX of the center of the planet along the x-axis with respect to the center of the detector. Note that for computational reasons, we fitted the apparent radius of the planet on the detector instead of the pixel FOV. Obviously, the actual apparent size of the planet in the sky is perfectly known from ephemeris databases. Note also that we need to fit H even after normalizing the photometry to the maximum value (as in Fig. 14) because the maximum of intensity onto the detector's discretized grid is not the actual maximum of the photometric profile. The error-bars on the free parameters were obtained from their posterior density function, sampled by the MCMC. This illustrates that even with a simple Lambertian model, the fitted model is very close to the measurements, except at the very edge. 5.9

6.2 6.5 0.9

1 1.1 1.2 1.3 -0.7 -0.4 -0.1 0.2
0 1 2 3 4 5.3 5.6 5.9 6.2 6.5 0.9 1 1.1 1.2 1.3 Fig. 14. Modeling Venus photometric profile with a Lambert law. Top panel, sub-frame: image of Venus obtained with the THEMIS solar telescope (Gaulme et al., in prep.). Top panel, main frame: the black line is the projection of the enclosed image along the x-axis, and the red line is the best-fit profile of it. Bottom panel: the above fitting was obtained by maximizing the likelihood of a Lambert-sphere model to the THEMIS observations, with an MCMC algorithm. Each sub-panel represents the two-parameter joint posterior distributions of all free parameters. The 68.3 %, 95.5 %, and 99.7 % confidence regions are denoted by three different gray levels. Probability Density Functions (PSFs) are plotted at the bottom and the right of the panels. The dotted lines mark the maximum probability value of the PDF of each parameter. After 200, 000 iterations, the fitted values are: 2.3 +0.5 -0.6 arcsec for the seeing, 5.89 +0.12 -0.23 arcsec for the apparent radius, 0.98 +0.07 -0.05 photometric amplitude (here normalized to the maximum value), and -0.15 +0.11 -0.16 arcsec for the planet's decentering with respect to the center of the detector. The red curve on the upper panel is the result of the fitted values.

Inversion of Doppler spectroscopic observations

Projection factors

In this section, we assume the above-described issues about location caused by the seeing are taken into account. We assume that whatever the technique (pinhole, long-slit, atomic cell, or Fourier transform spectrometry) we have a partial or full radial velocity map of the planet.

Spectral lines may be Doppler-shifted by planetary atmospheres for four possible causes: zonal "z", meridional "m", vertical "v", and/or subsolar-to-antisolar "SA" motions. To invert radial velocity maps, we need the projection factors which trans-Patrick Gaulme et al.: Measuring Planetary Atmospheric Dynamics with Doppler Spectroscopy form any motion into a radial velocity field:

∆V i (λ, φ) = f i (λ, φ) V i (λ, φ) (22)
where the subscript "i" is either "z", "m", "v", or "SA", and the coordinates (λ, φ) indicate the latitude and longitude in the planet system of coordinates. Details about projection factors are presented in Appendix C. In the case of reflected solar lines, the projection coefficient are expressed as follows:

f i = -u i .O -u i .S (23)
where the vectors O and S point towards the observer and the Sun, and u i is the vector describing the considered wind pattern of unitary amplitude, all expressed in the planetary frame (see Appendix C). Note that in the case we work with planetary emission or absorption lines, instead of solar reflected lines, the radial velocity with respect to the Sun is zero, i.e., f i = -u i .O.

Deprojection versus global modeling

At any given time, a radial velocity map is possibly the result of zonal, meridional, vertical or subsolar-to-antisolar motions.

It is obviously impossible to retrieve a value of each circulation component at each point of the planet from a single map. In other words, some approximations are needed to get quantitative information on atmospheric circulation. The choice of approximation depends on the type of dominant circulation in the considered atmosphere, as well as the level of detail we seek for. Global properties of atmospheric circulation can be obtained by globally fitting an ensemble of measurements, whereas deprojection -dividing an observed map by a projection factor -can be performed either when only one circulation component is significant, or on small scales in case of high-resolution observations from a spacecraft, for example. We now review what is appropriate for the planets of the solar system. Venus is the planet that has been studied the most with high-resolution Doppler spectroscopy in the visible domain. We know the circulation is dominated by a strong zonal wind, and a weaker equator-to-pole meridional circulation (e.g., Machado et al. 2017). In such a case, it is not accurate to simply divide a radial-velocity map by f z to estimate the zonal wind at any latitude and longitude, as it completely neglects the meridional circulation. It is more appropriate to assume a circulation model composed of a zonal and a meridional component, as the vertical flow is likely to be negligible on global scale. Then, different hypothesis can be tested and compared in terms of model likelihood. For example, a zonal component can be assumed as cylindrical (uniform speed with latitude), solid-body (uniform angular rotation), or more complex, as assuming an individual amplitude coefficients per latitude. Note that subsolar-to-antisolar winds were measured in Venus mesosphere (e.g., Lellouch et al. 2008), above 90 km. A possibly weaker return branch (antisolarto-subsolar) deeper could possibly affect atmospheric dynamics measurements at cloud-top level where visible observations probe (≈ 70 km). Like Venus, Titan has an atmospheric superrotation and similar modeling could be used. So far, its wind circulation was measured by Luz et al. (2005) in the visible, Kostiuk et al. (2001Kostiuk et al. ( , 2005Kostiuk et al. ( , 2006Kostiuk et al. ( , 2010) ) at 10-µm, Moreno et al. (2005) in the mm range. These were ground-based observations before or in support of the Huygens Probe measurements, and only zonal circulation was considered.

The case of Jupiter and Saturn are very different as these two planets are dominated by significant differential rotation: once removed the rotation of their interiors, their atmospheres are structured in bands and zones that rotate in opposite directions. Regarding meridional flows, because the band/zone alternate structure are small versions of Hadley cells, we do not expect global meridional circulation. Some intra band (or zone) meridional flows may exist but they are likely not measurable from Earth because of the relatively low spatial resolution3 . If local vertical motion of several tens of m s -1 are expected in some jovian planet thunderstorms (Hueso et al. 2002), vertical motions on large scales are a priori negligible with respect to zonal circulation. In such a configuration, where only zonal circulation is considered, deprojection can be applied. Each map can be divided by f z (λ, φ) (see details in next sub-section) and a value of the zonal wind is obtained in each point a any time. In practice, no observation is able to provide a reliable estimate from a single map: several hours are needed to be stacked to make the signal get out of the noise (see Gonçalves et al., submitted). A complete planisphere of zonal wind can be obtained by averaging data from several consecutive nights. This implicitly assumes that the zonal circulation does not evolve over a few days, and that there are no variations as function of local time either. This assumption is reasonable for Jupiter and Saturn, thanks to their fast rotation (∼ 10 hours), which allows for complete longitude coverage in two to three consecutive nights. At last, Uranus and Neptune are not easy targets to be observed from the ground, as their apparent diameters are barely larger than typical seeing conditions (about 3.7 and 2.2 arcsec respectively). No such band/zone structure is visible, even though jets are observed at several latitudes for both planets.

Implementing wind inversion

As it is not possible to obtain all the wind components from a single image, we describe the best possible method for interpreting a radial-velocity map, by taking into account the peculiarities of the atmospheric dynamics of the considered planet. On one hand, if we have good scientific reasons to expect the wind structure to be composed of simple functions of local coordinates, we can express the model as a function a limited number of parameters, and try to estimate these parameters. That way, it is possible to monitor possible temporal evolution based on observations taken at different times. On the other hand, for a one-component wind circulation, the best approach consists of deprojecting the radialvelocity map because the result is not limited by model assumptions.

Deprojection. Let us consider the case of a one-component global circulation pattern, where zonal circulation dominates. In principle, the recipe to deproject a radial velocity map into a zonal map is simple, as it consists of dividing the measured radial velocity map ∆V obs by the geometrical projection factor f z :

V z = ∆V obs f z , (24)
which was done, for example, by Machado et al. (2012) with observations of Venus with the UVES spectrometer of the Very Large Telescope (VLT). However, this is absolutely true only in absence of seeing degradation. In practice, the actual radial velocity map is described by Eq. ( 20). In the simple case of pure zonal circulation, and by expressing ∆V = f z V z , Eq. ( 20) rewrites as:

∆V obs = ( f z V z F) * P F * P (25)
where we dropped the (x, y) dependence to ease the reading. By assuming the angular rotation to be constant within the size of the PSF, we can take the term V z out of the parenthesis and get:

∆V obs ≈ V z ( f z F) * P F * P ; (26)
then:

V z ≈ ∆V obs f z C where C ≡ f z F * P ( f z F) * P . ( 27
)
The introduction of the C factor is the way to take into account the biases introduced by the atmospheric seeing, otherwise velocity maps are biased. Our ability to deproject is a matter of knowing the three parameters f z , F, and P. The term F is the photometric model, i.e. the actual image deconvolved from seeing degradation. The term P is the PSF and is estimated in the data processing steps.

The retrieved velocity map V z tends to infinite where the projection factor is close to 0. Therefore, a threshold needs to be applied to keep only values of the radial velocity that are reasonable. This is done by computing the noise map associated to the velocity map:

σ z = σ ∆V | f z | C (28)
The above equation shows that deprojection is not a good way to retrieve the properties of wind circulation because uncertainties tend to infinite where radial velocities tend to zero. We can circumvent this limitation whenever the data consists of a cube of images taken at different orientations of the planet. Within the assumption that the wind only depends on local coordinates of the planet and can be regarded as constant during all the observations, we can average several deprojected maps with different orientations and avoid to have regions with high levels of noise. To do so, individual deprojected maps are converted into spherical coordinates, and images with different rotation angle, i.e., looking at different longitudes, are averaged with a weighted mean, to discard data with larger noise level from the final map. Estimation of Jovian zonal wind maps obtained in this way will be presented in a future paper by Gonçalves et al. (in prep.). Global modeling. When atmospheric circulation needs to be described as the sum of several components (e.g., zonal and meridional), a global atmospheric model can be fitted on the whole map at once. The first step consists of minimizing χ 2 , the sum of the quadratic difference of the observations relatively to the model:

χ 2 = 1 n n ∆V obs -a z ∆V mod,z -a m ∆V mod,m -κ 2 , (29)
where "obs" refers to the observed radial velocity map, "mod" to the model, "z" to zonal, and "m" to meridional; n is the number of pixels on the planetary disk, a z and a m are the amplitudes of radial velocity components, and κ is a global offset in case of drift of the spectrometer. To take into account the atmospheric seeing effect, the models need to be priorly convolved with the PSF in order to minimize the χ 2 .

We just considered a very simple model with only three parameters, but the principle is the same for more complex models, as long as they linearly depend on the free parameters, in the form:

∆V obs = Hα + N, ( 30
)
where α is the vector containing the m parameters of the model to be estimated, and n is the number of measurements, where m < n. The matrix H is of rank m with a m × n size. The random noise N has independent components and null mean. An unbiased estimate of the parameters is obtained by minimizing the least square. Since the convolution by the seeing's PSF is a linear operation, this approach is still valid even when taking it into account. The convolution can be included in the model. If the noise is not uniform across the image -it is certainly the case when photon noise dominates as flux is not uniform -a weighting of the least square by the noise can be applied before minimization. The vector of estimate parameters α is retrieved by computing:

α = (H t H) -1 H t ∆V obs , (31)
or in case of weighted least square:

α = (H t WH) -1 H t W ∆V obs ( 32
)
where the weighting function W is a diagonal matrix where the diagonal terms are the inverse values of the noise of the data. The estimate α of the parameter vector α is truly unbiased, however, as the sensitivity of Doppler signal to the different components is null at some locations, some caution has to be taken to define the model that can effectively be recovered from the data. For instance, a model with a non-zero meridional component at the equator is not only unlikely for physical reasons, but also impossible to estimate for a mathematical point of view. In practice, we should take care that the Gramian matrix H t H is not close to be singular, because it would result in a strong increase of the uncertainties on the estimated parameters. As the same procedure provides an estimate of these uncertainties, it is possible to check the goodness of the fit. A poor normalized value of χ 2 indicates an inadequate model.

Conclusion and prospects

Despite the Doppler effect has been known for over 175 years, its application to planetary atmospheres in the visible has remained limited because of difficult implementation. Significant efforts in that matter have been led in the past decade, on one hand thanks to the ground-based follow-up organized to support the Venus Express mission (Lellouch & Witasse 2008;Widemann et al. 2007Widemann et al. , 2008;;Gaulme et al. 2008b;Gabsi et al. 2008;Machado et al. 2012Machado et al. , 2014Machado et al. , 2017)), on the other hand for the ground-based projects SYMPA and JOVIAL/JIVE dedicated to Jupiter seismic observations based on with visible Doppler imaging (Schmider et al. 2007;Gaulme et al. 2008aGaulme et al. , 2011;;[START_REF] Gonçalves | Proc. SPIE[END_REF]).

As it is somehow a new domain for planetary observations, the purpose of this paper was to consolidate the theory for rigorous Doppler spectro-imaging observations. We first revised and modified the theoretical prediction of the artificial rotation caused by the Sun's rotation, which was originally proposed and incompletely expressed by Young (1975). We revised the analytical approach, and produced numerical simulations which refine the estimated amplitude of the phenomenon by including the Sun's limb darkening and differential rotation. We also explored the impact of atmospheric seeing Patrick Gaulme et al.: Measuring Planetary Atmospheric Dynamics with Doppler Spectroscopy on planetary Doppler observations. It arises that for planets that are seen at non-zero phase angles, the convolution of the actual planetary image with the seeing PSF shifts it on the detector, which leads to significant biases when associating latitude and longitude reference frames to images of a planet. In addition, we studied the impact of atmospheric seeing on radial velocity measurements, which was pointed out by Civeit et al. (2005) for almost point-source objects (Titan, Io), and we showed it is significant as well for larger objects such as Venus or Jupiter. In the last section, we propose a roadmap for interpreting radial velocity maps, by including the Young and atmospheric seeing effects. We expressed the projection factors or any zonal, meridional, vertical, or subsolar-to-antisolar motion, and presented two possible approaches to retrieve the wind circulation: direct deprojection or global modeling of radial velocity maps.

This paper is the first of a series of three, the other two being in the submission process. Gaulme et al. (in prep.) present observations of Venus atmospheric dynamics obtained in 2009 with the MTR long-slit spectrometer of the THEMIS solar telescope located in Tenerife (Spain). Gonçalves et al. (submitted) presents observations of Jupiter with the JOVIAL/JIVE Doppler spectro-imager [START_REF] Gonçalves | Proc. SPIE[END_REF], which was specially designed for studying atmospheric dynamics of the planets of the solar system as well as performing seismology of Jupiter and Saturn. They present average zonal wind profiles of Jupiter obtained in 2015 and 2016 at the Calern observatory (France), and compare them with simultaneous Hubble-space-telescope cloud tracking measurements [START_REF] Johnson | Planetary and Space Science[END_REF]. sin λ cos λ cos φ

        (C.5)
A sub-solar to anti-solar circulation is described as in the frame centered on the sub-solar point, in which (λ S , φ S ) = (0, 0):

u SA, =         cos γ sin θ cos γ cos θ -sin γ         (C.6)
In the planet system of coordinates, it is expressed as:

u SA =        
cos λ S sin φ Ssin λ S cos λ S cos φ S cos φ S 0 sin φ S sin λ S sin φ S cos λ S sin λ S cos φ

S         u SA, (C.7)
To compute the radial velocity field of a planet, we need to describe each component of the circulation model as follows:

V z = Z(λ, φ) u z (C.8) V m = M(λ, φ) u m (C.9) V v = V(λ, φ) u v (C.10) V SA = S(γ, θ) u SA (C.11)
where the functions Z, M, V, and S carry the information about the wind structure. For example, a solid-body circulation is Z(λ, φ) = V eq cos λ, where V eq is the velocity amplitude at equator.

Fig. 3 .

 3 Fig.3. Schematic view of the radial-velocity map of the Sun as seen from Venus. The solar rotation axis is tilted by the angle θ. The system of coordinates (x, y) is the frame in which we compute the integral in Eq. 1, while (x , y ) is the system in which the solar rotation is described.

Fig. 4 .

 4 Fig. 4. The Young effect computed along Venus's equator on the morning hemisphere for solar-zenith angles γ = [85π/180, π/2 + ρ ]. The black dashed line indicates the Young effect in the case the Sun is above the horizon (Eq. 7). The plain red line corresponds with the case of the rising Sun (Eq. 8). The vertical line indicates the first contact of the Sun with the horizon. The minimum value of ∆V Y is -2, 000 m s -1 corresponds with the Sun's equatorial velocity V eq used for the simulation.

Fig. 5 .

 5 Fig.5. Numerical computation of the Young effect for Venus as function of the incidence angle γ and the solar rotation axis inclination with respect to the horizon θ, both expressed in degrees. Note the x-axis is not linear to enhance the large values of the effect at large incidence angles. Angles θ > 0 correspond to the afternoon hemisphere of Venus, and θ < 0 to the morning hemisphere. Contours and color scale are in m s -1 .

PatrickFig. 6 .Fig. 7 .

 67 Fig. 6. Numerical computation of the Young effect for Venus as seen from the Earth with a phase angle (Sun-Venus-Earth) of 40 • and a apparent diameter of 11.6 arcsec. The morning hemisphere is mostly visible. The 0 • m s -1 isovalue corresponds with the meridian crossing the sub-solar point of Venus. Contours and color scale are in m s -1 .

Fig. 8 .

 8 Fig. 8. Simulation of a Lambert sphere of 15 arcsec of diameter seen with an atmospheric seeing of 0, 1 and 3 arcsec respectively (black, red, blue curves), and phase angles α = [0, 45, 90, 135] • . Lambert spheres with no seeing degradation are normalized with respect to the intensity at the sub-solar point. The length of the horizontal dashed lines are equal to the actual size of the planet without seeing degradation. Their ordinates correspond to the photometric threshold applied to match the observed planet with the expected size along the equator. The right frames represent the Lambert sphere without (top) and with a 3-arcsec seeing degradation (bottom). The red line indicate the edges of the sphere without seeing degradation.

Fig. 9 .

 9 Fig.9. Top panel: apparent diameter of Venus (arcsec) as function of phase angle (degrees), from June 2017 to June 2027 (black dots). Bottom panel: same for Jupiter, with axis reversed with respect to Venus, and a longer time span (2017-2117) to average enough orbits.

Fig. 10 .

 10 Fig. 10. Apparent shift of a Lambert sphere in the conditions Venus (top) and Jupiter (bottom) are seen from Earth. For Venus case, the shift along the equator relatively to the planet's radius δX/R is plotted as function of the phase angle (Sun-Venus-Earth). For Jupiter, δX/R is plotted as function of the apparent diameter. From top to bottom, the consecutive lines correspond to different seeing conditions, from 0.5 to 4 arcsec, as indicated near the minimum of each curve. The black indicate the planet and phases. Note that proportions are conserved in Venus and Jupiter panels.

Fig. 11 .

 11 Fig. 11. Radial velocities of a Lambert sphere of 15 arcsec of diameter, with a retrograde rotation of 1 m s -1 seen with an atmospheric seeing of 0, 1 and 3 arcsec respectively (black, red, blue curves), and phase angles α = [0, 45, 90, 135] • . We assume that radial velocities are measured on solar reflected lines (Fraunhofer), which explains why radial velocities can be larger than ±1 m s -1 . Left panels are cuts along the equator and right frames are maps. Top maps are the true radial velocity signal, and bottom are the same degraded by a 3-arcsec atmospheric seeing.

PatrickFig. 12 .

 12 Fig.12. Top: map of the Young effect on Venus as seen from the Earth when the atmospheric seeing's PSF is 2 arcsec, the phase angle (Sun-Venus-Earth) 40 • and the apparent diameter 11.6 arcsec. Bottom: cut of the apparent Young radial velocity V Y along the equator in the same conditions as above, with (red) and without (black) blurring by atmospheric seeing.

Fig. 13 .

 13 Fig.13. Left panels: cut of the radial velocity map of a Lambert sphere rotating in the retrograde direction (solid-body rotation of 1 m s -1 at equator), observed with a phase angle of ±40 • . Radial velocities are measured on solar reflected lines. The width of the PSF used to degrade the image and radial velocity map is 35 % of the planetary radius, which is equivalent to a 2-arcsec atmospheric seeing for Venus when it is observed at ±40 • (see Fig.9). The term δX refers to the apparent displacement of the planet on the sky caused by the atmospheric seeing.

Table 1 .

 1 Coefficient of the Young effect for Venus, Jupiter, Saturn, Uranus, and Neptune, for an optical wavelengths of 350, 550, and 800 nm from the Allen (1973) limb darkening laws. Orbital semi-major axes are used as Sun-to-planet distances.

	Planet	Y (350 nm) Y (550 nm) Y (800 nm)
	Mercury	5.07	5.37	5.59
	Venus	2.71	2.88	2.99
	Mars	1.29	1.36	1.42
	Jupiter	0.38	0.40	0.42
	Saturn/Titan	0.20	0.22	0.23
	Uranus	0.10	0.10	0.11
	Neptune	0.07	0.07	0.07

This comment about possible inaccurate pixel FOV refers to the case of Venus observations byGaulme et al. (in prep.) where the MultiRaies spectrometer of the THEMIS solar telescope was used in a specific configuration. The pixel FOV was know to about 5 % and no calibration was done to measure it.

Note that if high-resolution planetary images are obtained from the ground by amateur astronomers, it is only with lucky imaging techniques, i.e., with exposure time less than 0.1 second. Radial velocity measurements are the results of stacking hours of observations, composed of individual exposures of at least 30 sec, and are thus limited by seeing.

Acknowledgments

P. Gaulme was supported in part by the German space agency (Deutsches Zentrum für Luft-und Raumfahrt) under PLATO data grant 50OO1501. The JOVIAL project is supported by the Agence Nationale de la Recherche under the contract number ANR-15-CE31-0014-01. I. Gonçalves's PhD is granted by Observatoire de la Côte d'Azur and the JIVE in NM project (NASA EPSCoR grant #NNX14AN67A). THEMIS is UPS of the CNRS funded by INSU. The authors thank E. Lellouch for his very constructive review of the paper as the referee. P. Gaulme wishes to precise that he did most of the writing while listening to the 1969 Velvet Underground Live.

A&A proofs: manuscript no. Gaulme_Planet_Doppler_Spectro

Appendix A: Latitude and longitude of a 3D oblate spheroid

We consider the planet-centered system of coordinates (x, y, z), defined such as the y axis points at North (Fig. A.1). This frame is useful for avoiding confusion when dealing with data where an image of a planet is defined as function of the (x, y) axis of a detector (e.g., the simulated Venus in Fig. 6). Rotating planets or stars can be approximated with an oblate spheroid, i.e., axes along x and z are the same. For an oblate spheroid, the latitude λ and longitude φ on the planet are simply derived from the planetary-centered Cartesian frame:

In the case of a sphere ρ = R; for an oblate spheroid ρ expresses as:

where R eq is the equatorial radius and µ the oblateness factor defined as R pole = µR eq .

Appendix B: A Lambert sphere illuminated by an extended source

A Lambert sphere assumes that incoming light is isotropically scattered by its surface. A light coming from a given direction reaches the surface with an intensity that is proportional to the cosine of the incidence angle γ, defined as the angle between the incoming light direction and the normal to the surface. The intensity of a Lambert sphere is thus:

As illustrated in Fig. 2, the intensity received from the Sun in each point of a Lambert sphere is the solar photometry, which can be described by a limb-darkening law, modulated by the variation of incidence angle across the solar surface:

Appendix C: Projection factors

The radial velocity map of a planet in the case of reflected solarlines is the result of the relative motion at the surface of the planet with respect to both the Sun and the observer: ∆V = -V.O -V.S (C.1) where V is the velocity vector in the planet system of coordinates, O is the unitary vector pointing at the observer, and S is the unitary vector pointing at the Sun. The negative sign arises from the usual convention that a redshift corresponds with a positive radial velocity.

Let us consider a sphere, for which (λ O , φ O ) and (λ S , φ S ) are the coordinates of the sub-observer and sub-solar points in the planetary frame (Fig. A.1). The vectors O and S are expressed as: