# Atmospheric circulation of Venus measured with visible imaging spectroscopy at the THEMIS observatory

Patrick Gaulme, François-Xavier Schmider, Thomas Widemann, Ivan

Gonçalves, Arturo López Ariste, Bernard Gelly

► To cite this version:

Patrick Gaulme, François-Xavier Schmider, Thomas Widemann, Ivan Gonçalves, Arturo López Ariste, et al.. Atmospheric circulation of Venus measured with visible imaging spectroscopy at the THEMIS observatory. Astronomy and Astrophysics - A&A, 2019, 627, pp.A82. 10.1051/0004-6361/201833627. hal-02327296

HAL Id: hal-02327296

https://hal.science/hal-02327296

Submitted on 24 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Atmospheric circulation of Venus measured with visible imaging spectroscopy at the THEMIS observatory*

Patrick Gaulme<sup>1,2,3</sup>, François-Xavier Schmider<sup>4</sup>, Thomas Widemann<sup>5,6</sup>, Ivan Gonçalves<sup>4</sup>, Arturo López Ariste<sup>7</sup>, and Bernard Gelly<sup>8</sup>

Received 13 June 2018 / Accepted 24 May 2019

## ABSTRACT

Measuring the atmospheric circulation of Venus at different altitudes is important for understanding its complex dynamics, in particular the mechanisms driving super-rotation. Observationally, Doppler imaging spectroscopy is in principle the most reliable way to measure wind speeds of planetary atmospheres because it directly provides the projected speed of atmospheric particles. However, high-resolution imaging spectroscopy is challenging, especially in the visible domain, and most knowledge about atmospheric dynamics has been obtained with the cloud tracking technique. The objective of the present work is to measure the global properties of the atmospheric dynamics of Venus at the altitude of the uppermost clouds, which is probed by reflected solar lines in the visible domain. Our results are based on high-resolution spectroscopic observations with the long-slit spectrometer of the solar telescope THEMIS. We present the first instantaneous “radial-velocity snapshot” of any planet of the solar system in the visible domain, i.e., a complete radial-velocity map of the planet obtained by stacking data on less than 10% of its rotation period. From this, we measured the properties of the zonal and meridional winds, which we unambiguously detect. We identify a wind circulation pattern that significantly differs from previous knowledge about Venus. The zonal wind reveals a “hot spot” structure, featuring about 200 m s⁻<fnref n="1" /> at sunrise and 70 m s⁻<fnref n="1" /> at noon in the equatorial region. Regarding meridional winds, we detect an equator-to-pole meridional flow peaking at 45 m s⁻<fnref n="1" /> at mid-latitudes, i.e., about twice as large as what has been reported so far.

**Key words.** planets and satellites: individual: Venus – planets and satellites: atmospheres – methods: observational – techniques: imaging spectroscopy – techniques: radial velocities

## 1. Introduction

The atmosphere of Venus is well known for its super-rotation in a retrograde direction. The atmospheric zonal rotation period strongly varies with altitude, from a corotation with the surface of the planet of 243.02 Earth days at ground level, down to 4.4 days at cloud tops, where it peaks at about 100 m s⁻<fnref n="1" /> at equator. First evidenced from the ground (Boyer & Guérin 1969), the atmospheric super-rotation has been extensively studied both from space and ground-based telescopes (Gierasch et al. 1997; Limaye et al. 1988; Rossow et al. 1990). The cloud top region is important as it constrains the global mesospheric circulation in which zonal winds generally decrease with height while thermospheric subsolar-to-antisolar (SSAS) winds increase (Lellouch et al. 1997; Widemann et al. 2007, 2008). It also shows important variability at various spatial and temporal scales (Sánchez-Lavega et al. 2008; Hueso et al. 2012, 2015; Patsaeva et al. 2015; Khatuntsev et al. 2013; Machado et al. 2012, 2014).

Characterizing the meridional circulation is also important for understanding the maintenance of the super-rotation, by determining the global mean and eddy circulations and the associated meridional transport of angular momentum and energy (e.g., Limaye & Rengel 2013). The role of thermal tides to transport angular momentum vertically in low latitudes has been confirmed (Lebonnois et al. 2010; Takagi & Matsuda 2007). It has also been noticed that the latitudinal distribution of zonal wind at cloud tops may result from an equilibrium between the impact of thermal tides and the angular momentum transport by the meridional circulation (Lebonnois et al. 2010), providing grounds for systematic and simultaneous observations of both zonal and meridional regimes.

For Venus, as for dense atmospheres in the solar system, most atmospheric dynamics measurements come from the cloud tracking technique. The method consists of following cloud features at specific wavelengths taken on image pairs obtained at various times. Although the clouds are almost featureless in visible light, there are prominent features in UV and infrared wavelengths (Titov et al. 2008, 2012). Cloud motions are considered to be a good proxy for true atmospheric motions and are capable

<footnote n="1"><sup>1</sup> Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany
e-mail: gaulme@mps.mpg.de</footnote>

<footnote n="2"><sup>2</sup> Department of Astronomy, New Mexico State University, PO Box 30001, MSC 4500, Las Cruces, NM 88003-8001, USA</footnote>

<footnote n="3"><sup>3</sup> Physics Department, New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA</footnote>

<footnote n="4"><sup>4</sup> Laboratoire Lagrange, Université de Nice Sophia Antipolis, UMR 7293, Observatoire de la Côte d'Azur (OCA), Nice, France</footnote>

<footnote n="5"><sup>5</sup> LESIA, Observatoire de Paris, PSL, CNRS, Sorbonne Université, Université Paris-Diderot, France</footnote>

<footnote n="6"><sup>6</sup> DYPAC, Université Versailles-Saint-Quentin, Université Paris Saclay, France</footnote>

<footnote n="7"><sup>7</sup> IRAP, Université de Toulouse, CNRS UMR 5277, CNES, UPS, 14, Av. E. Belin. 31400 Toulouse, France</footnote>

<footnote n="8"><sup>8</sup> THEMIS Observatory, La Laguna, Tenerife, Spain</footnote>

<footnote>* Tables A.1–A.3 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/627/A82</footnote>
of providing a systematic long-term monitoring of the atmo-
spheric winds (e.g., Sánchez-Lavega et al. 2008; Peralta et al.
2008, 2012; Moissl et al. 2009). Outstanding measurements were
obtained by the ESA Venus Express mission (hereafter VEx;
Svedhem et al. 2007), whose main goal was a better understand-
ing of the atmospheric circulation, with a specific attention to the
origin of the super-rotation. Cloud tracking measurements were
provided by the Venus Monitoring Camera (VMC; Markiewicz
et al. 2007) and the Visible and InfraRed Thermal Imaging
Spectrometer (VIRTIS; Drossart et al. 2007).

Despite the exquisite quality of VEx measurements, cloud
tracking indicates the motion of large cloud structures (limited
by spatial resolution), which is an indication of the speed of
iso-pressure regions rather than the speed of the actual cloud
particles. In the case of Venus, cloud stuctures may represent the
phase speed of a condensation wave, possibly associated with
vertical mixing or chemical processes associated with the UV
absorber (Widemann et al. 2007; Hueso et al. 2015; Machado
et al. 2014, 2017). A complementary solution to access direct
wind speed measurement is Doppler spectrometry because it
measures the actual speed of cloud particles and has different
but complementary limitations in local time, latitudinal, and
temporal coverage than an orbiting spacecraft. The idea of com-
paring cloud tracking and Doppler spectroscopic measurements
has emerged with the idea of supporting the VEx, and eventu-
ally Venus Climate Orbiter (Akatsuki), which entered the Venus
orbit in 2015 (Nakamura et al. 2016).

Characterizing the atmospheric circulation of Venus based
on high-resolution spectroscopy is actually an old idea, but it
was logically considered to be challenging (e.g., Moreux 1928,
Chap. V, p. 101). The first reliable measurements with mod-
ern spectrographs started in the 1970s (Traub & Carleton 1975;
Young et al. 1979). Starting in 2007, the ground-based sup-
port to VEx kicked off many observational projects (Lellouch &
Witasse 2008), including the present project. Significant results
on the upper mesospheric dynamics were obtained using mid-
infrared heterodyne spectroscopy (Sornig et al. 2008, 2012)
and millimeter and submillimeter wave spectroscopy (Clancy
et al. 2008, 2012; Lellouch et al. 2008; Moullet et al. 2012),
but Doppler spectroscopy is more challenging at shorter wave-
lengths. Visible observations of solar Fraunhofer lines scattered
by Venus clouds were performed by Widemann et al. (2007,
2008), Gabsi et al. (2008), Gaulme et al. (2008), and Machado
et al. (2012, 2014, 2017). Most results regard average zonal
wind profiles as a function of latitude at the cloud-top level.
Widemann et al. (2007, 2008) and Machado et al. (2014, 2017)
also reported to have measured instantaneous zonal and merid-
ional wind circulation from data obtained with the ESPaDonS
échelle spectrometer of the Canada France Hawaii Telescope
(CFHT), which we discuss in this paper (Sect. 2).

Back in 2007, we proposed to use the THEMIS<fnref n="1" /> solar tele-
scope to get Doppler maps of Venus by scanning the planet
with the 100 arcsec long-slit spectrometer MulTiRaies (MTR),
whose resolution ranges from 100 000 to 1 000 000 in the visi-
ble (Mein & Rayrole 1985). Two factors motivated this choice.
First, the use of a long slit allows for reconstructing com-
plete radial-velocity (RV) maps of the planet, which is not
possible with single-fiber spectrographs, as ESPaDonS, or spec-
trographs with shorter slits, as that of the Ultraviolet and Visual
Echelle Spectrograph (UVES) of the Very Large Telescope

(VLT)<fnref n="2" />. Secondly, a solar telescope can easily observe during
the daytime, which is convenient for Venus. We led two test cam-
paigns in 2007 and 2008, and an “actual” campaign in 2009,
when the observing setup and data processing technique were
fully developed. The 2007 test did not meet good weather con-
ditions and the results were extensively discussed in a previous
publication (Gaulme et al. 2008, hereafter G08). Only a portion
of the visible phase of Venus was scanned in the 2007 obser-
vations. A rough estimate of the zonal wind was obtained, i.e.,
151 ± 16 m s⁻<fnref n="1" />, which was significantly larger than previously
measured. In 2008, the observing setup consisted of repeat-
edly scanning Venus from west to east. This observing protocol
showed that the mean RV would drift with time, for instrumen-
tal reasons, which introduced an uncontrollable bias in the zonal
wind estimate, as we could not disentangle RV variations due
to zonal winds from instrumental drifts. During the next run,
in September 2009, we opted to alternatively scan Venus from
north to south (NS) and south to north (SN) with the slit parallel
to Venus equator to circumvent this issue. That way, as described
later, we could ensure reliable RV values at each latitude and
erase possible drifts by combining NS and SN scans.

In this paper, we report the results of the 2009 campaign,
where we obtain for the first time an instantaneous and spa-
tially resolved RV map of any planet of the solar system. We
first review the previous works dealing with Venus RV mea-
surements, and we bring into question some aspects of what
was done so far in relation with our own observing approach
(Sect. 2). We then detail the data acquisition and processing tech-
niques that lead to photometric and RV image reconstruction
(Sect. 3). Next, we present the RV maps obtained on differ-
ent days and our empirical model, which involves zonal and
meridional winds as well as a local-time dependence of the
zonal component (Sect. 4). Finally, we discuss our results and
compare them to both cloud tracking and recent spectroscopic
observations (Sect. 5).

## 2. Two decades of efforts

2.1. Expected Doppler map

The Doppler technique has been extensively used from the
ground with a variety of observational techniques involving mil-
limeter, infrared, and visible spectroscopy. Spectroscopic tracers
of the atmospheric dynamics of Venus include: CO and isotopic
¹³CO rotational lines in the millimeter wave range to probe at
105 and 95 km, respectively; infrared emission of CO₂ in non-
local thermodynamical equilibrium (110 km level); and solar
Fraunhofer lines and Venus CO₂ molecular lines in the visible
probing cloud top level (70 km) and a few kilometers above,
respectively. Heterodyne spectroscopic measurements have been
interpreted as a combination of SSAS and zonal winds (e.g.,
Lellouch et al. 1994, 2008; Clancy et al. 2008; Sornig et al.
2008), with variable relative contributions. The cross-terminator
SSAS velocity is usually inferred to be on the order of 100 m s⁻<fnref n="1" />,
while equatorial zonal velocity varies between 0 and 200 m s⁻<fnref n="1" />.
Poleward meridional winds may also have been marginally
detected (Lellouch et al. 1994, 2008).

From visible observations, the atmospheric circulation of
Venus is mostly zonal at the top of cloud layers, peaking at about
100 m s⁻<fnref n="1" /> at the equator (e.g., Sánchez-Lavega et al. 2017). The
zonal wind was originally considered to follow either a solid

<footnote n="1">¹ Télescope Héliographique pour l'étude du Magnétisme et des Insta-
bilités Solaires.</footnote>

<footnote n="2">² In Machado et al. (2012), the slit length of UVES was 11 arcsec versus ≈20-arcsec for the diameter of Venus.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Schematic view of the Doppler shift of a prograde rotating planet as seen from the THEMIS observatory. The Doppler effect measured in the visible on the reflecting cloud deck is the sum of the motion relative to the Sun and the Earth. Radial velocities are zero along the bisector meridian, located halfway in between the subsolar and subterrestrial meridian.</figcaption></figure>

body rotation, i.e., the same rotation period at all latitudes, or cylindrical, i.e., the same velocity at all latitudes, except in the polar regions where it should drop to zero. Cloud tracking and RV measurements led to a slightly different picture; an average zonal wind profile as a function of latitude looks like the letter “M”: wind rising from 0 to 120 m s⁻¹ from poles to 45° latitude and decreasing to 100 m s⁻¹ around the equator (e.g., Machado et al. 2017). A meridional circulation corresponding to equatorial-to-pole Hadley cells was measured both from cloud tracking and RV measurements, peaking at about 20 m s⁻¹ at 45° latitude (e.g., Machado et al. 2017).

The Doppler shift of solar Fraunhofer lines reflected on a planet is the sum of the RV of the planet relatively to the Sun and to the observer. In case of a planet at exact opposition, the Doppler effect of reflected solar lines is then doubled. In the case of a non-zero phase angle, i.e., the angle Sun-planet-observer, the Doppler shifts cancel each other on the meridian located at the bisector of the subsolar and subterrestrial points (e.g., Gabsi et al. 2008). A retrogradely rotating zonal circulation therefore displays a blueshift in the morning and a redshift in the afternoon (Fig. 1). In our case, observations were performed during Earth morning elongation, which means that we were seeing the morning terminator of Venus. In Fig. 2, we represent theoretical RV maps of a solid-body rotator and a meridional circulation based on two Hadley cells at the phase angle corresponding to the 2009 campaign. On top of this, uniform RV offsets are expected, corresponding to the relative motion of Venus with respect to the Sun and mostly of Venus with respect to the observer on the rotating Earth. The motions are well documented in the ephemeris database and must be taken into account.

Beyond real RV fields, the rotation of the Sun as seen from Venus introduces a bias in RV measurements, as originally introduced by Young (1975) and subsequently completed by Gaulme et al. (2018). Rays from a different part of the Sun, which show different RVs, reach the planet with (slightly) different incidence

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Simulation of RV maps corresponding to a model of solid-body rotation, i.e., $V = V_{\text{equator}} \cos(\lambda)$ (top), and a meridional circulation, i.e., $V = V_{\text{1=45°}} \sin(2\lambda)$ (bottom), where $\lambda$ stands for the latitude. For the zonal wind $V_{\text{equator}} = 100 \text{ m s}^{-1}$, for the meridional wind $V_{\text{1=45°}} = 20 \text{ m s}^{-1}$. The Doppler shift assumes the use of reflected solar lines and a phase angle of 43°. The black boxes indicate the orientation of the entrance slit of the spectrometer and the arrows show the scanning direction. The slit width corresponds with the exact width (0.75 arcsec) relative to Venus, but the length was reduced from 100 to 18 arcsec for illustration purposes.</figcaption></figure>

angles. Regions of the Sun that are closer to the horizon contribute less to the reflected solar spectrum than regions closer to zenith. Thus, the RV integrated over the whole solar disk is not zero at a given point of Venus. In other words, even if Venus were not rotating, we would still measure a Doppler shift near Venus terminator, and that Doppler shift would mimic a retrograde rotation because the solar rotation is prograde. Gaulme et al. (2018) demonstrated that the RV field $\Delta V_Y$ associated with the so-called Young effect is expressed as

$$ \Delta V_Y(\gamma, \theta) = Y(\Lambda) \tan \gamma \sin \theta, \qquad (1) $$

where $\gamma$ is the solar-zenith angle, $\theta$ the inclination of the solar spin axis with respect to local horizon, and $Y(\Lambda)$ a coefficient that is about 2.9 m s⁻¹ at $\Lambda = 550$ nm. The analytical expression by Young (1975) of the artificial Doppler shift $\Delta V_Y$ on Venus was calculated for the equator and did not include the $\sin \theta$ term. In addition, his expression did not include the solar limb darkening and the Sun’s differential rotation, which leads the coefficient by Gaulme et al. (2018) to be smaller by about 10% than that of Young. We note that the new expression of the
Young effect extended to all latitudes by Gaulme et al. (2018) questions the wind measurements that have been done so far near Venus’s terminator, which made use of the Young (1975) analytical formulation (Widemann et al. 2007, 2008; Gabsi et al. 2008; Machado et al. 2012, 2014, 2017).

The observed RV map is not only the sum of all the above-listed contributions, it is also strongly biased by the atmospheric seeing. Firstly, as extensively studied by Gaulme et al. (2018), the atmospheric seeing modifies the apparent location of the planet in the sky whenever the planet is not observed at full phase (opposition). This leads to biases when inverting RV maps because it affects the position of longitudes and latitudes on the Venus image. Secondly, the seeing convolves regions of variable RV and photometry, which tends to reduce the apparent amplitude of atmospheric motions. As originally pointed out by Civeit et al. (2005), the resulting RV map is the convolution of the RV signal with the photometric map of the considered object, including its degradation by seeing. The mean Doppler $\Delta V_{\text{obs}}$ measured in a given pixel $(x, y)$ on the detector can be expressed as

$$ \Delta V_{\text{obs}}(x, y) = \frac{(\Delta V F * P)(x, y)}{(F * P)(x, y)}, \qquad (2) $$

where $\Delta V$ and $F$ are the real RV and photometric maps prior to seeing degradation, $P$ is the point spread function (PSF) of the atmospheric seeing, and the asterisk sign * indicates the convolution product. This effect tends to reduce the amplitude of RV variations from west to east. Besides, it cancels out most of the Young effect, making it almost negligible (Gaulme et al. 2018). In this paper, we take into account the data degradation by atmospheric seeing to extract zonal and meridional wind circulations, for the first time with RV measurements of Venus performed in the visible.

## 2.2. Review of previous works

In this section, we review the recent measurements of the atmospheric circulation of Venus with Doppler spectroscopy in the visible domain. This includes the works that started after 2000, slightly before or together with the VEx ground-based support, which were shortly described above. Techniques used for Doppler velocimetry in the visible solar spectrum on the dayside are mostly high-resolution échelle spectroscopy with single optical fiber feeding, i.e., single aperture measurements (Widemann et al. 2007, 2008; Gabsi et al. 2008; Machado et al. 2014, 2017), but also long-slit spectrometry (Gaulme et al. 2008; Machado et al. 2012). In the framework of coordinated campaigns to support VEx science investigations, the major breakthrough in terms of observational techniques for visible high-resolution spectroscopy of Venus has been led by Widemann and Machado (Widemann et al. 2007, 2008; Machado et al. 2012, 2014, 2017). We review all of these works, from the pioneering observations of Widemann et al. (2007), Gaulme et al. (2008), and Gabsi et al. (2008) to the robust observational protocols of Machado et al. (2017).

### 2.2.1. Bushwhacking: early works

Among pioneering works, we consider Widemann et al. (2007), Gabsi et al. (2008), and Gaulme et al. (2008) apart from other works, as they constitute the first attempts of the past two decades to measure the winds of Venus with high-resolution spectrometers in the visible wavelength. As for any new project, their initial observational protocols and early analysis assumptions made their results significantly discrepant with respect to

later studies. These three works nevertheless contributed toward renewing interest in the technique and kicking off the ground based-support to VEx and Akatsuki, which were fundamental steps toward reaching later successful measurements.

Widemann et al. (2007) reported measurements of the average global winds of Venus with the AURELIE high-resolution spectrometer at the 1.52 m telescope of Observatoire de Haute Provence, France. These authors measured the Doppler effect on Venus CO$_2$ absorption lines and a few reflected solar Fraunhofer lines in the range 8600–8800 Å. Their results confirmed the existence of a zonal retrograde flow, even though the measured mean equatorial velocity of 75 ± 15 m s$^{-1}$ was relatively low with respect to posterior measurements, and strong day-to-day variations (± 65 m s$^{-1}$) were identified. By combining the results from all data, these authors also reported the possible detection of a SSAS circulation component of amplitude of about 40 m s$^{-1}$ at the terminator. Given the consistency of observational results reported since then, it is possible that the daily variations of zonal winds were artifacts because most observations realized later on with improved protocols never displayed anything similar. Nevertheless, this work had the great merit of paving the way for the methods – especially the sequential pointing – used in Widemann et al. (2008) and Machado et al. (2014, 2017).

The Gabsi et al. (2008) observations were performed with the EMILIE high-resolution, cross-dispersed spectrograph and its associated calibrating instrument the Absolute Astronomical Accelerometer (AAA), at Observatoire de Haute-Provence, France. From their three best observing nights, they reported zonal wind values of 75 ± 6, 85 ± 3 and 91 ± 6 m s$^{-1}$ without considering the Young effect. These authors reported a better stability by introducing the Young effect according to the Young (1975) formula and retrieved a mean zonal circulation of 48, 47 and 51 ± 3 m s$^{-1}$, which is far from what has been measured otherwise. No other observations of Venus with this instrument have been published since then.

Gaulme et al. (2008) reported a test of the long-slit MTR spectrometer of the THEMIS solar telescope to scan Venus and retrieve a complete RV map of it, which had never been done before. The system planned for scanning Venus did not work and the weather was poor. The result is a partial RV map that was modeled with a global zonal circulation pattern (151 ± 16 m s$^{-1}$) without considering the Young effect. The analysis of these preliminary data was promising enough to justify further observations, which constitutes the base of the present work.

### 2.2.2. Stable radial-velocity measurements

The series of papers by Widemann et al. (2008) and Machado et al. (2014, 2017) shows a remarkable consistency both in terms of methods and results. All observations were done with the ESPaDonS spectro-polarimeter at the 3.6 m CFHT, in coordination with VEx/VIRTIS-M in 2011 and 2014 (Machado et al. 2014, 2017). The spectrometer ESPaDonS covers the whole visible range (3700–10 500 Å) at an average resolution of 80 000 and is fed by an optical fiber whose field of view (FOV) is 1.6 arcsec. The observation protocol was sequential, like in most other ground-based wind measurement techniques, either with CO rotational lines or infrared (IR) CO$_2$ non-LTE emission lines. Only Moullet et al. (2012) performed Doppler mapping with interferometric observations based on CO lines at the Plateau de Bure interferometer.

The sequential acquisition consists of selecting a set of positions on the Venus dayside and making one observation at a time per position to cover all positions. The whole scan of the dayside
hemisphere is then repeated once or twice, according to observ-
ing conditions. Centering and guiding is manually controlled
with the help of a Venus template that is taped on the display of
the guiding camera (Fig. 1 of Machado et al. 2017). To compen-
sate RV drifts of about 100 m s⁻¹, the spectrometer’s wavelength
calibrations were carried out using both Thorium-Argon (ThAr)
lamps and a set of telluric lines following the standard protocol
developed at CFHT (Donati et al. 1997). Final RV measurements
were corrected from the various motions (Venus–Sun, Earth–
Venus), as well as the Young effect according to the Young
(1975) expression.

In all three papers, the complete RV dataset is interpreted as
a horizontal zonal circulation, where two possible regimes are
considered: a solid-body or a cylindrical regime. Except for the
preliminary observations at CFHT by Widemann et al. (2008),
which suffered of a lack of stability (zonal winds in between 92
and 155 m s⁻¹), the results by Machado et al. (2014, 2017) pro-
vide very consistent values of the mean zonal wind at equator
from day to day and in between 2011 and 2014; i.e., 117.3 ± 18.0
and 117.5 ± 14.5 in Machado et al. (2014) and 119.6 ± 16.5,
122.6 ± 31.3, 119.6 ± 26.0, and 118.1 ± 19.5 in Machado et al.
(2017) with 2σ errors. Beyond mean zonal winds, Widemann
et al. (2008) introduced a meridional wind component as part of
their measurements, while Machado et al. (2014, 2017) looked
for the presence of equator-to-pole meridional circulation on the
spectra that were taken – on purpose – along the bisector merid-
ian, where zonal RV signal is canceled. Machado et al. (2014,
2017) reached a conclusion on the detection of a meridional
circulation peaking at about 20 m s⁻¹ at mid-latitudes, thus con-
firming the results obtained with Venus Express cloud tracking
measurements (Sánchez-Lavega et al. 2008; Hueso et al. 2012,
2015; Khatuntsev et al. 2013; Machado et al. 2014, 2017).

Machado et al. (2012) reported wind measurements per-
formed with the long-slit spectrometer UVES at the VLT in
Paranal, Chile. The UVES spectrometer covers the whole visi-
ble range (3000–10 000 Å) with an average spectral resolution
of 100 000 and a spatial resolution of 0.2 arcsec per pixel. The
slit size was set to 11 × 0.3 arcsec, while the diameter of the
planet was in between 20 and 22 arcsec with phase angles of
82.9° and 87.8° during two runs in May and June 2007, respec-
tively. The advantages of this instrument are its large band pass
(like ESPaDonS), large collecting area (8.2 m telescope), small
pixel FOV, and mapping capability. These authors focused on a
few localized positions on the planet: three that were parallel to
the rotation axis and six parallel to the equator. The positions
along the rotation axis aimed at putting into light any asym-
metry of zonal rotation in between both hemispheres. However,
such configuration does not allow for determining the absolute
speed of zonal circulation. To the contrary, the slit parallel to
the equator aimed at providing absolute measurements of the
amplitude of zonal circulation and possible longitudinal varia-
tions. Intentionally, no meridional circulation was considered at
the time.

For each spectrum, the Doppler shift was computed by cor-
relating the spectrum with the spectrum at the center of the slit,
as previously done by Luz et al. (2005, 2006) and Civeit et al.
(2005) with observations of Titan. In other words, the actual ref-
erence spectrum used for spectral calibration therefore changed
from position to position of the slit. To check for instrument
stability and to correct for optical slit curvature, exposures of
the built-in ThAr lamp were taken after the VLT/UVES science
exposures. To this set of corrections, the additional subtrac-
tion of the Young effect from the Young (1975) expression was
performed near the terminator.

From the spectra obtained with the slit parallel to the equa-
tor, Machado et al. (2012) reported an average zonal circulation
in between 106 and 127 ± 14 m s⁻¹. In addition, they investigated
possible variations of zonal wind amplitude as a function of local
time. They identified a slight increase of the wind speed near
the terminator, which corresponded to the evening side during
Spring 2007, at about 150 m s⁻¹ at 10° away from it. From the
spectra obtained with the slit parallel to the rotation axis, they
detected a slight asymmetry of zonal circulation by measuring
that winds are faster by 6 ± 5 m s⁻¹ in the southern hemisphere.
The M-shaped latitudinal profile of mean zonal circulation is
later compared with other published methods (e.g., Fig. 14 of
Machado et al. 2017).

2.2.3. A few concerns and the big picture

We base our concerns on a recent study by Gaulme et al. (2018), who in particular have shown that atmospheric seeing introduces biases regarding both the localization of the planet on the detector and the RV field. None of the cited works fully takes into account these effects. As regards observations by Machado et al. (2014, 2017) it is unlikely that it significantly influenced the final wind determination, thanks to the good seeing conditions. Nevertheless, we consider that neglecting the bias on RV by atmospheric seeing has likely led to underestimating the speed of zonal circulation next to the terminator. In the same way, making use of the Young (1975) expression of the Young effect is erroneous for data taken out of the equatorial region and may have biased part of the results. However, we note that the CFHT observations were done relatively far from the terminator, where the seeing effects on both RVs and the Young effect are small relative to the amplitude of zonal circulation. To the contrary, both UVES/VLT and MTR/THEMIS observations included the terminator region.

More specifically, because of this potential seeing effect, we
question an aspect of the Machado et al. (2012) analysis, where
the detection of small-scale longitudinal wave structure as a
function of local time on Venus is considered (RV variations
of about 10–30 m s⁻¹ over spatial scales of 216 km). Firstly,
a scale of 216 km corresponds to two pixels and is at the very
limit of the Shannon criterion on signal sampling. Secondly,
with a pixel FOV of 0.2 arcsec and an atmospheric seeing larger
than 1.2 arcsec, detecting such small-scale variations sounds
optimistic.

Regarding the CFHT observations (Machado et al. 2014,
2017), we identify two minor issues. Firstly, the fit of a pure
zonal circulation intentionally neglects the fraction of the RV
signal attributed to the meridional component. This choice could
alter the retrieved speed of the zonal wind. However, consider-
ing the relatively low amplitude of meridional wind reported so
far (~20 m s⁻¹ at mid latitudes), this approximation likely intro-
duced a marginal bias. Secondly, Fig. 4 of both papers show
local variations of the zonal wind as a function of longitude and
latitude, which were obtained from the sum of the model plus
the residuals (Machado, priv. comm.). This result is question-
able because local values of the zonal wind should be obtained
by dividing the residuals by the zonal-circulation projection fac-
tor, which is null at the bisector meridian. This would result in
divergent error bars on zonal wind speed along that meridian.
Still, we note a good agreement in between the Machado et al.
(2014, 2017) results with the simultaneous cloud tracking data
of VIRTIS-M/VEx. In the present paper, we follow the approach
proposed by Gaulme et al. (2018) in which we perform a global
fit of the observed RVs, including a zonal plus a meridional
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Example of spectral images obtained with THEMIS in 2008 while scanning Venus. Color scale is inverted, i.e., dark is bright and white is faint. The three spectral ranges are centered on Mg, Fe, and Na Fraunhofer lines in the mid-visible domain. The 2D image is a zoom of the actual image where the dark area corresponds to Venus. The Fraunhofer lines are the vertical white structures. Since this spectra were taken during daytime, the absorption lines are visible both on Venus and the Earth sky. For all of the deepest lines, we clearly distinguish the Doppler shift in between Venus and Earth atmosphere. In these figures, the slit was positioned parallel to Venus' equator close to the sub-Earth meridian. The black line, corresponding to the y-axis on the right side, is the Venus spectrum at about zero latitude.

component. This does not allow for a direct comparison of our
measurements with the local wind values of Machado et al.
(2014, 2017).

Overall, we retain that both observations conducted at CFHT with the fiber-fed spectrometer ESPaDonS and at VLT with the UVES long-slit spectrometer led to identifying a mean zonal circulation of about 120 m s⁻¹ in between ±45° parallels with an M-shaped profile. The UVES deprojected RV measurements per latitude indicated a slight increase of the zonal circulation toward terminator. The CFHT measurements along the bisector meridian show the presence of a double Hadley-cell equator-to-pole circulation with amplitude peaking at about 20 m s⁻¹ at mid-latitudes.

**3. Telescope, instrument, and methods**

**3.1. THEMIS telescope with the MTR spectrometer**

The THEMIS observatory is a solar telescope dedicated to accurate measurement of polarization of solar spectral lines with high spatial, spectral, and temporal resolutions (Mein & Rayrole 1985). It is a 90 cm diameter Ritchey–Chrétien telescope. As for G08, it has been operating in the MTRmode (Mein & Rayrole 1985) with no polarimetric analysis, which simultaneously permits spectral observations in up to four different spectral domains. Spectrometry with a slit produces 2D images, whose horizontal component is the optical spectrum and vertical component the spatial dimension (Fig. 3). The slit is 100 arcsec long, 0.75 arcsec wide, and the spectral resolution was set to R = 150 000. The guiding and positioning on the planet was controlled by a tip-tilt mechanism.

Using this solar observatory instead of a classical night-time telescope has two advantages. Firstly, a high-resolution spectrometer with such a long slit is not common at all. This allows us to get a whole cut of any planet at once; all planets have an apparent diameter lower than 60 arcsec – including a significant fraction of Earth skylight, which is helpful to monitor the stability of the spectrometer, as we see later. Secondly, it allows us to not be limited by the sunlight since it is designed to stare at the Sun. We could for instance observe Venus at phase angles that are not accessible with night-time telescopes, as in 2008, where the phase was about 13°, i.e., the elongation 9°.

As in G08, we worked with solar lines that are reflected
on the cloud decks of Venus. In the present observations, three
10 Å broad spectral domain were considered, all in the central

part of the visible spectrum and each centered around one or
several deep and sharp Fraunhofer lines: the magnesium line
at 5173 Å, the iron doublet at 5573 Å, and the sodium D1 at
5896 Å. The fourth detector was tested to get Doppler shifts on
CO₂ molecular lines, but the S/N of these data was not good
enough to go further. The theoretical velocity sensitivity of a
single line is roughly proportional to the slope of the line and
the total amount of photons. We refer G08 for the detailed esti-
mate. The maximum theoretical velocity sensitivities of each line
are 32, 17, and 42 m s⁻¹ arcsec⁻¹ min⁻¹ for Mg, Fe, and Na
respectively.

The observing protocol consisted of repeatedly scanning the
planet with the slit, either from NS or from SN. We would
place the slit at a given location on Venus and acquire several
exposures in a row. Then we would move to another position,
and so on until the planet was completely covered. Steps from
one position to the next were set to 0.8 arcsec, which roughly
corresponded to the slit width. Technically, it was actually the
opposite, as we positioned Venus with respect to the slit and not
the other way around. In practice, the position of the slit on the
planet was done by tilting the “M5” mirror, which sits at the
pupil plane of the telescope and whose angle is controlled by
a piezo-electric device. The rotation of the planet with respect
to the slit was obtained by rotating both the derotator and the
entrance slit by 180°. As we show in Sect. 4.2, we found out that
the steps increased while moving on the planet from SN, which
means that this occurs as a function of the elongation. This is
not fully understood but can be explained if one of the key com-
ponents was not perfectly aligned (e.g., M5 not exactly at pupil
plane, or derotator not perfectly aligned). In the following, we
indicate as a “scan” a set of observations that includes spectra
from one edge of the planet to the opposite (either from NS or
the reverse). Scans were typically composed of 16 consecutive
positions on sky – enough to account for the blurring caused by
atmospheric seeing – with 10 exposures per position.

**3.2. Calibration data and preprocessing steps**

Calibration files are fundamental to ensure high quality Doppler
measurements. In particular, errors on flat fields can induce
distortions of the shape of spectral lines along the spatial dimen-
sion. In other words, it can produce fake Doppler shift across
the planet. A quasi-absence of flat fields was one of the main
limitations of the data quality of G08, where a fixed Doppler pat-
tern on the detector was identified and removed a posteriori from
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Example of mean solar spectrum (**left panel**) and resulting flat field (**right panel**) on the detector dedicated to the Fe line during the 2009 campaign. A mask was applied to the redressed flat field. Vertically, it corresponds to regions that Venus never crosses, while horizontally, the image quality was affected by the rectification/interpolation.</figcaption></figure>

the RV map. We considerably improved our protocol during the 2008 and 2009 campaigns.

On the one hand, dark fields were obtained before each scan, and a master dark field was obtained by averaging them together. Using individual dark fields or the master dark field did actually not make any significant difference on RV measurements. On the other hand, obtaining a flat field is somewhat delicate with a slit spectrometer, given that daylight presents spectral lines. Using incandescent lights into the dome was not possible because of the specific configuration of this solar telescope. We chose to acquire about 400 spectra by directly pointing at the Sun to get a very high S/N spectral image on the whole detector by averaging them all (Fig. 4). We note that the 400 spectra were taken at random locations on the solar disk to average out irregular Doppler shifts in the solar photosphere. Indeed, spatial and spectral resolutions of THEMIS is high enough to clearly distinguish by eye distortions of the Fraunhofer lines due to the granulation and presence of p-waves close to the surface.

These three spectral images were then used to measure and correct the geometrical distortions on the detector: instead of being vertical, absorption lines appeared to be slightly bent (Fig. 4). We determined row by row the center of each absorption line from the image derivative with respect to the spectral axis (x-axis). A third-order polynomial fitting was used to fit the measured position of the bottom of the lines and to recify the images with a cubic spline interpolation algorithm. We note that distortion is different for each detector and was estimated and corrected independently. The optical distortions were stable enough to not repeat this calibration in between observations. However, for cautiousness, we repeated the process every day and did not notice any significant variation. An average line profile could be computed by vertically collapsing the rectified image of the average solar spectrum. Each row of the solar spectral image was finally divided by the average line profile to obtain a flat-field image (see Fig. 4).

The dark and flat fields therefore obtained were subtracted and divided, respectively, from each spectral image of Venus. One last important item to preprocess the data taken during daytime, i.e., most of the data, consists of removing spectral lines that are scattered by the Earth's sky. For this, we selected the area of the detector where Venus is not present and collapsed all rows to get an average background line profile, as we did for the flat field. This contribution was then subtracted from the spectral image. This process had to be repeated for each spectral image because the solar lines scattered by Earth's atmosphere drift with

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Example of data preprocessing on an image taken during the 2008 campaign on May 5. Color scale is inverted, i.e., dark is bright and white is faint. The image is centered on the Fe doublet at 557 nm. *Left panel:* raw image, which is the result of a 1 min exposure on Venus during daytime. The dark area corresponds to the region occupied by Venus spectrum, while the Fraunhofer lines scattered by Earth's sky are clearly visible out of it. The sky lines are tilted because of the geometric distortions mentioned in the text. The areas delimited by dotted lines and indicated by two vertical arrows are used to compute the average sky spectrum, which is then subtracted from the whole image. *Right panel:* “cleaned” image after dark field removal, rectification of geometric distortion, flat-field division, and mean sky lines subtraction.</figcaption></figure>

Earth's rotation. The whole process was demonstrated to work adequately as illustrated in Fig. 5.

To ensure a good data quality at each position of the slit, several exposures were taken every time, six in 2008 and ten in 2009. Thanks to the excellent pointing stability due to the tip-tilt mechanism, and after inspecting the images, we considered that Venus was not moving during the six or ten exposures at each slit position. In other words, shifts along the spatial dimension could happen while translating the slit on Venus from one position to another, but not otherwise, except if a cloud hid Venus in the meantime. The six or ten images taken at the same position were then averaged. The final product, i.e., ready to be used for measuring Doppler shifts, consisted of rebinning each of these mean spectral images along the spatial dimension. Indeed, the pixel FOV (0.20 arcsec) of the MTR is well below the seeing – from 2 to 4 arcsec in 2009 – so that we could gain a factor two in S/N by rebinning four pixels along the y-axis, leading to a rebinned pixel resolution of 0.8 arcsec.

### 3.3. Measuring Doppler shifts

Measuring absolute RVs on any astronomical object is very challenging, and we seek these RVs for tens of m s⁻¹. This has been discussed in the case of Venus wind measurements in Young et al. (1979) and Widemann et al. (2007, 2008) with attempts of making absolute RV measurements using visible lines. These authors concluded the need of a reference point on Venus that would serve as a relative velocity reference, and they used such a point to perform differential velocity measurements on the disk.

In our case, we measured Doppler shifts by comparing the position of a given spectrum on Venus with a high S/N reference spectrum, which was not one of Venus spectra. Instead, the reference spectrum is the mean solar spectrum used to build a reference flat field. Indeed, the solar spectrum scattered by Venus atmosphere is the result of integrating all upwelling and downwelling flows on the solar surface. The average of 400 spectra randomly located on the solar disk makes it representative of what Venus receives and reflects, and has the advantage of displaying a much larger S/N (≈500) than any individual spectrum obtained on Venus. With that method, the estimated shifts
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Two-step process used to measure Doppler shifts. Top panel: the reference spectrum is the average of 400 spectra obtained on the Sun, while spectrum of Venus corresponds to a single location on Venus disk. Bottom panel: same spectra after shifting the reference spectrum at the same location of Venus and after normalizing them. The wavelength ranges indicated by two arrows and vertical dotted lines indicate the selected areas where cross-correlation is applied (see text). Spectra were obtained on September 14, 2009.</figcaption></figure>

need to be corrected from the motion of Venus with respect to observer and the motion of the observer with respect to the Sun. All these components are well known and can be subtracted with the help of ephemeris data retrieved from the website of the Institut de Mécanique Céleste et de Calcul des Ephémerides<fnref n="3" />. The RV differences between Venus and the Sun with respect to the Earth were about 3 and 9 km s⁻¹ in 2008 and 2009, respectively, i.e., 0.06 and 0.18 Å in the mid-visible. We note that a typical Fraunhofer line is about 0.3 Å wide at half maximum.

We employed a two-step process (Fig. 6) to estimate the Doppler shift. Let us consider a spectrum obtained on Venus. We first computed the cross-correlation of the Venus spectrum with the reference solar spectrum. To give more weight to the spectral ranges that are sensitive to Doppler shifts, i.e., those with the steepest spectral slope, we multiplied each spectrum by the absolute value of its derivative before computing the cross-correlation. Then, we interpolated the reference spectrum on a grid shifted by this first estimate of the Doppler shift. We note that we interpolated the reference spectrum and not the Venus spectrum because it shows a good enough S/N to not be altered by interpolation. Interpolation was performed with a spline algorithm.

A second step is necessary because even though we give more weight to the sensitive parts of the spectral range, the contribution of the noisy continuum still alters the measurement.

Now that the two spectra are almost overlapping, we computed the cross-correlation of the two spectra on the cores of the Fraunhofer lines, instead of the complete spectral range. In the case of the Fe doublet, we selected two regions on the detector (Fig. 6), then cross-correlated the spectra on each region, and took the average.

We note that we considered using the Connes method to measure the Doppler shifts, which is commonly used in exoplanetary science (Connes 1985). However, this method is best suited when working with thousands of lines, while we are working with a maximum of nine lines. We still tested it instead of performing a second cross-correlation on the core of the lines but the results were clearly noisier, as seen by the presence of many outliers.

Last but not least, we estimated the error on Doppler velocity values from a mix of measurements and simulated data. For a given spectral range (Mg, Fe, or Na), we first built a model spectrum from the solar spectrum used for making the flat fields, which we smoothed over three spectral bins. Secondly, we simulated 100 000 of simulated spectra with a given S/N (either 10, 50, 100, or 1000) based on our model spectrum and a normally distributed noise, with a given Doppler shift (either 10, 50, 100, 200, or 1000 m s⁻¹). For each simulated spectrum, we ran our Doppler shift estimation routine. The histograms of the estimated Doppler shifts with respect to their exact values gave us an estimate of the measurement error. We observed that the error was independent from the absolute value of the Doppler shift and is – as expected – a linear function of the inverse of the S/N. For a S/N of 100, the average error on Doppler velocity estimates was written as

$$ \sigma_{v,Mg,ref} = 17 \text{ m s}^{-1}, \qquad (3) $$

$$ \sigma_{v,Fe,ref} = 18 \text{ m s}^{-1}, \qquad (4) $$

$$ \sigma_{v,Na,ref} = 32 \text{ m s}^{-1}, \qquad (5) $$

for the three spectral ranges we consider in this paper. Then, for each actual spectrum, we measured its S/N from the standard deviation of the data minus the model spectrum in the continuum regions, and we retrieved the velocity error by comparing the S/N to the reference simulated values, i.e.,

$$ \sigma_v = \sigma_{v,ref} \frac{100}{S/N}. \qquad (6) $$

For the data taken on September 17, 2009, which include Mg, Fe, and Na spectra, the velocity error $\sigma_{v,avg}$ on an average map is written as

$$ \sigma_{v,avg} = \frac{1}{3} \sqrt{\sigma_{v,Mg}^2 + \sigma_{v,Fe}^2 + \sigma_{v,Na}^2}, \qquad (7) $$

where $\sigma_{v,Mg}$, $\sigma_{v,Fe}$, and $\sigma_{v,Na}$ are the error maps in each spectral band.

# 4. Observations

## 4.1. Observation setting

As mentioned in the introduction, we were awarded observing time from May 2 to May 8, 2008 to extend the test reported in G08. The planet diameter was about 9.8 arcsec and the phase angle $\phi \approx 13^\circ$. This campaign ended up being a test campaign as it put into light some difficulties in reaching our goal, as we identified a spurious Doppler shift in the measurements that was not from astrophysical origin. For this run, we opted for an exclusive west-east (WE) scanning of the planet in the Venus coordinate

<footnote n="3">³ http://www.imcce.fr</footnote>
<table><caption>**Table 1.** Observation properties of Venus during the 2009 campaign at Teide observatory.</caption><thead><tr><th>Date</th><th>Start (MJD)</th><th>End (MJD)</th><th>RA (h m s)</th><th>Dec (°′′)</th><th>Dist. (AU)</th><th>V (km s<sup>-1</sup>)</th><th>Phase (°)</th><th>Elong. (°)</th><th>R (arcsec)</th><th>Scan #</th><th>Seeing (arcsec)</th></tr></thead><tbody><tr><td>2009.09.14</td><td>55 088.300350</td><td>55 088.633799</td><td>9 39 05</td><td>14 44 35</td><td>1.403</td><td>8.688</td><td>43.10</td><td>29.12</td><td>5.95</td><td>12</td><td>2.2–3.5</td></tr><tr><td>2009.09.16</td><td>55 090.359586</td><td>55 090.445014</td><td>9 48 40</td><td>14 01 20</td><td>1.413</td><td>8.529</td><td>42.29</td><td>28.64</td><td>5.91</td><td>2</td><td>3.0–3.7</td></tr><tr><td>2009.09.17</td><td>55 091.357399</td><td>55 091.502237</td><td>9 53 26</td><td>13 39 05</td><td>1.418</td><td>8.449</td><td>41.89</td><td>28.40</td><td>5.88</td><td>6</td><td>3.1–5.2</td></tr></tbody></table>

**Notes.** Columns 1–3: observation days, start, and end dates. Start and end times are in modified Julian date (MJD), i.e., Julian date JD −2 400 000.5 days. Columns 4–10: ephemerides from the IMCCE database. Ephemerides are calculated every day of both campaign at 06:00:00 UTC. The last two columns indicate the number of scans done each day and the atmospheric seeing.

frame, i.e., slit parallel to rotation axis, which was repeated over and over as long as the weather permitted. The reason to choose a WE scanning was the possibility to construct a map easily from individual spectral images. By building a map, we actually mean two maps: a “photometric” image and an RV image of Venus. We intend by photometric maps an image of Venus in the visible, but we do not aim to quantify the exact photon flux.

With a WE scan, the photometric map can be obtained in a straightforward manner. Firstly, at each position of the slit, we considered the mean spectral image (result of averaging ten images), which we projected along the *x*-axis to get a 1D vector that contains the north-south (NS) photometric profile at the slit position along the equator. Knowing that slit positions are spaced by 0.8 arcsec, we just placed each vector one after the other into a 2D table and produced a map. However, this is not enough because the slit sometimes moves along the *y*-axis in between two positions, introducing vertical shifts on the detector. The fact of having a WE scan makes the vertical adjustment easy because all photometric profiles must be symmetrically centered on the apparent equator of Venus. Both the photometric and RV maps were then interpolated according to the shift measured from the photometric profiles. However, this choice, which was motivated by its simplicity for reconstructing the map, ended up being terrible. Indeed, the RV measurements appeared to be dominated by a signal that was not related to Venus but to the instrument configuration. The spectrometer had never been tested for measuring such small RV shifts, so the drift was a surprise. We had no way to a posteriori disentangle the instrumental drift from the Venus wind circulation because both the instrumental drift and the zonal circulation were oriented along the WE direction. This is why we rotated the slit by 90° and scanned Venus along the NS direction during the following campaign.

The second campaign occurred from September 11 to 17, 2009, when Venus was displaying a phase of about 43° and an apparent diameter of 11.9 arcsec. The weather allowed us to scan Venus 12, 2, and 6 times on September 14, 16, and 17, respectively (Table 1). On September 14 and 16, only spectra of the Fe doublet were taken, while spectra in the three ranges were taken the other day. In addition to scanning the planet along the NS direction we also scanned it along the reverse direction (SN) to help characterize the RV bias induced by the spectrometer. If truly deterministic, averaging the NS and SN maps would have allowed us to get rid of the bias without needing to model it.

## 4.2. Long path to assemble radial velocity maps

The drawback of the NS/SN scanning approach is a challenging map construction. With a WE scan, the middle of the photometric profile at each position of the slit corresponds to the equator. With an NS scan, the middle of the photometric profile has no specific meaning for a planet that is not observed at opposition.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Mean photometric profile from September 14, 2009 data. A so-called photometric profile is the sum of each photometric map projected along the *y*-axis on the *x*-axis. Intensity was normalized such as maximum value is 1. The “NS” and “SN” photometric profiles are the average profiles for each scanning direction (NS or SN). The symmetric profile is the bisector of both profiles.</figcaption></figure>

It is even impossible to connect the location of the maximum of intensity to a given longitude, as its position is a function of the atmospheric seeing, which we do not know (Gaulme et al. 2018). To be able to associate a latitude and a longitude on Venus to a given position along the slit, we must consider that the photometric profile is the result of the actual photometric profile, the latitude on the planet, the pixel FOV, and the atmospheric seeing. We note that for our specific instrumental configuration, the pixel FOV was known with an accuracy of about 5%, which we aimed at refining for accurate wind measurements. If we had a complete photometric map of the planet, we could estimate both the atmospheric seeing and the apparent diameter (i.e., pixel FOV), by fitting a photometric profile degraded by a PSF, assuming the true photometric limb-darkening profile to be known.

Unfortunately, we do not have yet a photometric map and we need to know the seeing and the pixel FOV to build it. Where to start then? The first observational input consists of computing the photometric profile along the NS and the SN directions, directly from the data. For each position of the slit, we already have a photometric profile along the EW direction. If we collapse these photometric profiles, we get a total intensity as a function of the projected latitude. We can retrieve a proxy of the pixel FOV and the atmospheric seeing from this simple profile. In principle, modeling this profile should be enough to extract the parameters.

However, our planned processing strategy needed to be tuned because of an unexpected issue. In Fig. 7, we plot the NS and SN photometric profiles of Venus, which are the average of the
12 scans done on September 14, 2009. They appear to be both asymmetric and in opposite directions, which means that the interval between each position of the slit increase as we move it. The farther the slit is from its original position, the larger is the step. This stretching along the NS direction thus prevents us from simply modeling the photometric profile and recombining NS and SN maps. We first need to symmetrize the photometric profile and then interpolate the maps on an evenly sampled grid of slit position. We compute the median photometric profile as the bisector of the NS and SN profiles (plain black line on Fig. 7), which is symmetrical and we can model as the result of a limb-darkening law modulated by atmospheric seeing and pixel FOV. We note that the atmospheric seeing was not monitored by any device at the observatory, especially because most observing time was during daytime.

Fitting a photometric profile that is a function of the pixel FOV and atmospheric seeing (convolution by a Gaussian function) is not straightforward. We developed a dedicated routine based on a Monte Carlo Markov Chain (MCMC) optimization method, including the Metropolis-Hasting algorithm and parallel tempering, and Bayesian inference (Gaulme et al. 2018). At each step of the iterative process, a 2D photometric profile is generated on a grid of pixels with free FOV, then the image is convolved with a Gaussian function of free full width at half maximum (FWMH). The center of the planet is also a free parameter. The image is then projected along the NS direction and compared with the actual profile. We minimize the square of the difference (least-squares fitting) until convergence is reached. The error bars are retrieved from the posterior density function of each parameter. We note that we make use of a Lambert law to describe the limb darkening because the smooth aspect of Venus in the visible, coupled with a small apparent diameter (~5 times the seeing value), does not justify models that are undistinguishable from a Lambert sphere in our observing conditions. In Fig. 8, we present the fitting of the symmetrized photometric profile along the NS direction, done with the MCMC routine for the September 14 data. Seeing is estimated to be $3.21_{-0.78}^{+0.67}$ arcsec and apparent radius$^4$ on the detector $5.54_{-0.22}^{+0.17}$ arcsec if the steps were 0.80 arcsec. It actually means that the steps were of 0.86 arcsec.

This step is however not sufficient to estimate accurately the seeing and pixel FOV because the atmospheric seeing varies with time and the symmetric photometric profile is an average of the scans. Besides, the fact of having symmetrized the profile could insert a small bias in the apparent dimension of the planet upon the detector. The second step thus consists of attributing a latitude to each position of the slit based on the symmetric NS profile. For each step, the photometric cut along the EW direction is fitted in the same way as the NS profile with an MCMC routine. This process allows us to know the location of the center of the planets at each step of the scan and to monitor the seeing along the day (Fig. 9). On September 14, where we had the best conditions (no clouds), the seeing ranged from about 2 to 3.5 around noon. On September 16 and 17, the seeing conditions were definitely worse with a mean value above 3 arcsec.

Once the relative position along the EW direction are determined, we place each photometric vector in a 2D matrix and interpolate (linear interpolation) it on the regular grid along the NS direction, which we got by symmetrizing the profile.

<figure><img src="image_9.png" /></figure>

**Fig. 8.** Modeling Venus photometric profile with a Lambert law. Top panel: the black line is the mean photometric profile of Venus along the NS direction, and the red line is the best-fit profile. 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 subpanel 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. The PDFs are plotted at the bottom and right of the panels. The dotted lines indicate the maximum probability value of the PDF of each parameter. After 500 000 iterations, the fitted values are $3.21_{-0.78}^{+0.67}$ arcsec for the seeing, $5.54_{-0.22}^{+0.17}$ arcsec for the apparent radius, and $0.03_{-0.11}^{+0.10}$ arcsec for the planet decentering with respect to the center of the detector. The red curve on the top panel indicates the result of the fitted values.

The same is done with RV maps, except that the interpolation is performed in a way to not alter the velocity value at the edges. Indeed, by definition the RV field is zero out of Venus (we do not compute it), so that simply interpolating a column on a shifted vertical axis would bias velocities at the edges. Therefore, we extend the RV vector out of Venus with the last value met on the west for the northern space above Venus, and vice versa for the eastern part of the detector. After interpolation, all regions out of Venus are set back to a default zero RV.

### 4.3. Biases in the radial-velocity measurements

Most of our data were acquired during daytime. Therefore, we can monitor the possible drift of the spectrometer thanks to the solar spectrum that is scattered by the Earth's blue sky. We note that the MTR spectrometer of the THEMIS telescope was not designed to perform high precision RV measurements – spectro-polarimetric observations of the Sun instead – and is not thermally stabilized. Drifts of hundreds of meters per second are expected.

<footnote n="4">⁴ 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.</footnote>
<figure><img src="image_10.png" /></figure>

**Fig. 9.** Atmospheric seeing measured by fitting a Lambert limb- darkening law degraded by atmospheric seeing, as function of time. Time origin is set at the first exposure of the day. The red curve indicates a fit of it with a third order polynomial.

Figure 10 shows the mean RV that is measured on the sky. The variation of RV is much larger than that expected from the Earth's rotation (gray lines). Even though we do not have measurements of the temperature of the spectrometer, it is very likely that the variations we observe are dominated by temperature, which is typical of spectrometers. The amplitude of the observed variations is about 3000 m s⁻¹ (≈0.06 Å) both for observations performed on September 14 and 17, 2009. We can correct from these variations by subtracting the mean RV that is measured on sky, or at least a smooth of it to reduce the measurements uncertainties on the relatively low S/N sky RVs. For the only scan that was done during night time (the first of September 14th), we extrapolated the drift of the sky from a fitting of the six following scans, where the trend is well fitted by a second order polynomial. Again, such a drift was expected and correcting it is straightforward.

What was not expected and surprised us is an RV drift in between the Earth's sky and Venus. In Fig. 11, we plot the

<figure><img src="image_11.png" /></figure>

**Fig. 10.** Mean RV of the Fraunhofer lines scattered by the Earth's sky as functions of time. Each chunk (separated by empty spaces) corresponds to a scan. The first scan taken on September 14, 2009 is not shown as it was taken before sunrise. In both plots, the thin gray line represents the theoretical variation of the sky-scattered solar lines according to the ephemeris (Earth's rotation essentially). The abbreviation “cam” refers to the CCD camera's number.

<figure><img src="image_12.png" /></figure>

**Fig. 11.** Difference of RVs in between Venus and the Earth's sky. The horizontal line on the *bottom panel* corresponds to a spatial coverage of 10 arcsec. The abbreviation “cam” refers to the CCD camera's number.

difference of the Earth's sky RVs with the mean RV that is measured on Venus at each position. Since we are scanning the planet along the north-south axis (slit parallel to equator), we do not expect large variations of the average RV along the slit. It should display the meridional circulation on top of relative motions of the planet with respect to the observatory and the Sun. Actually, it arises that whatever the direction of the scan (NS or SN), there is a spurious RV that increases monotonically at each step of the scan. For Fe spectra, the RV drift per scan reaches 2000 m s⁻¹ per scan on September 14, while it is of 400 m s⁻¹ three days later (Fig. 11).

In 2008, we scanned the planet along the WE direction and found out this bias. This was particularly terrible as it overlapped and overwhelmed the zonal circulation of Venus. That made the observation taken in 2008 impossible to use scientifically. We thus decided to scan Venus in the perpendicular direction during
the 2009 observations. That way, the bias does not compete with
the zonal circulation of Venus, and we can combine the RV maps
obtained by scanning the planet from SN with those done from
NS to hopefully cancel this systematic effect.

The origin of such a bias is still unclear and we have no
definitive explanation. A possible explanation is related to the
way we scan the planet. As mentioned earlier, positioning the
slit on the planet actually consists of positioning the planet with
respect to the slit thanks to a titled mirror (M5). If the M5 mir-
ror is not perfectly placed at the pupil focus of the telescope, the
small rotation of the mirror applied to position the planet on the
slit can cause some astigmatism. In other words, the optical dis-
tortions of the field, which we measure on the solar flat field (see
Fig. 4), could slightly vary while the mirror is moving. In such
a case, an imperfect rectification of geometrical distortion could
lead to a differential RV drift in between Venus and the sky.

To correct the RV maps from this bias, we first flatten each
individual RV map by fitting a tilted plane (2D linear polyno-
mial) oriented along the NS axis (Fig. 12). Indeed, if the RV bias
is pretty much stable from one scan to the next, it actually seems
to diminish slightly with time. Therefore, by subtracting a tilted
plane to the RV map is a way to include the amplitude varia-
tion of the bias. Subtracting a tilted plane affects any uniform
pole-to-pole meridional RV field, but not equator-to-pole merid-
ional fields. When we then combine the RV maps, we assume
the residuals to cancel each other. This assumption is the main
weakness of the whole process as we cannot be fully sure that
the bias completely cancels out when averaging scans done in
both directions. We discuss that aspect in Sect. 6 regarding the
strong variations of zonal wind speed as a function of local time
that we observe.

In Fig. 12, we show the individual flattened RV maps and they appear to be very consistent from day to day and spectral line to line, although the data taken on September 14 are the only data that benefited from good observing conditions. From the photometric images, it clearly appears that the data taken on September 17 suffered from worse atmospheric seeing and transparency (clouds). We note RV discontinuity during severe photometric drops caused by clouds, which is especially visible in the second half of the data taken on September 17. Conditions on September 16 were poor and allowed only for two scans (NS and SN). The error bars on Doppler velocity are about twice as bad as for the other two nights. For each day, all flattened RV maps were then averaged to get a final RV map (Fig. 13). For obvious observing condition differences, we only consider the data taken on September 14 to model the atmospheric circulation. However, it is worth noticing that despite significantly worse conditions and strong biases in some latitude ranges, the overall aspects of the September 16 and 17 RV maps are very consistent with the September 14 map.

**5. Models of atmospheric circulation**

In this section, we aim to interprete the RV map obtained
on September 14, 2009 (Fig. 13, upper panel). We first note
that this map was computed by stacking data over about eight
hours, which represents a little less than 10% of Venus’ rota-
tion at that altitude. We can therefore consider this map as
an “instantaneous” Doppler snapshot of Venus. As far as we
know, it is the first instantaneous optically resolved RV map
of any planet of the solar system in the visible domain. The
RV image of Jupiter obtained by Gonçalves et al. (2019) with
the JOVIAL/JIVE instrument are averaged over several complete
planetary rotations.

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Photometric and RV maps for September 14, 16, and 17, 2009. Data were taken on Fe lines only on September 14 and 16, and on Fe, Na, and Mg lines on September 17. Each RV map was “flattened” by subtracting a fitted plane along the NS direction, i.e., x-axis. Photometric intensity is expressed in ADU and RV in m s⁻¹.</figcaption></figure>

Because we expect, or at least we cannot exclude, the pres-
ence of both zonal and meridional winds and the variations of
winds with local time, we must extract the atmospheric circu-
lation pattern with forward modeling. It consists of fitting the
RV map with a simulated RV map that is composed of a sum of
different terms (zonal, meridional, etc), which we alter accord-
ing to the atmospheric seeing. The model is optimized with an
MCMC routine that we developed to interpret these data. Such
an approach, with respect to classical least-squares fitting, is a
straightforward and reliable way to compute error bars on the
estimated parameters. Two examples of outputs of the MCMC
code are shown in Fig. 14. These examples correspond with two
models that we comment on the next paragraphs.

We consider several types of types of models, from a sim-
ple solid-body rotator to circulation patterns involving variations
of wind speed with local time and hemispheric asymmetries.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** *From left to right:* photometric, Doppler, and Doppler error maps obtained during the 2009 campaign. **Top:** September 14, middle: September 16, and bottom: September 17. **Middle panels:** colors indicate redshift/blueshifts and the color bar indicates the Doppler shift in m s⁻¹. **Right panels:** color indicates the error on Doppler measurements σ<sub>V</sub>. Spatial coordinates are expressed in arcsec. The white and black contours indicates the theoretical size of the planet on the detector and the maps are oriented such as north is on top. With respect to Fig. 12, the images and maps were averaged, then rotated by 90° to align the vertical axis with Venus polar axis. Pixel size is ≈0.8 × 0.8 arcsec. We note the change in scale for Doppler errors: from 8 to 21 m s⁻¹ on Sept. 14, 17 to 41 m s⁻¹ on Sept. 16, and 11 to 21 m s⁻¹ on Sept. 17. The data corresponding to the 9 maps displayed in this figure are available in Tables A.1–A.3.</figcaption></figure>

We show that if a simple solid-body rotator is not sufficient to properly fit our observed map, models that are too sophisticated cannot be discriminated as the S/N level of the RV map does not allow us to look for small-scale structures. We note that the atmospheric seeing is about 3 arcsec in average on the map and the apparent diameter 11.9 arcsec.

The first impression we get by looking at the RV map is that it looks similar to the simulated solid-body rotator from Fig. 2 (top panel), which is expressed as

$$ \Delta V_{\text{model 1}} = V_{z,0} \cos \lambda + \kappa, \qquad (8) $$

where $V_{z,0}$ is the zonal wind amplitude at the equator, $\lambda$ is the latitude (from 0 to $\pi/2$ in each hemisphere), and $\kappa$ is an offset. We note that the offset term $\kappa$ is meant to take into account the fact that our measurements are not absolute. We introduce the subscript “0” in $V_{z,0}$ to be consistent with other models (3–7), where more than one term is employed to describe the zonal circulation. However, fitting the map with a pure solid-body zonal circulation pattern plus a global RV offset (Model 1, see Fig. 15 and Table 2) appears to be rather unsatisfactory given the structure of the residuals (Fig. 15, top row, middle panel), and the large value of the reduced $\chi^2$ at 6.14. To help evaluate the fitting quality, we also plot both data and model across the planet’s equator in the

third column of Fig. 15. The disagreement with a solid rotator is obvious. Furthermore, the residual map shows a symmetrical feature in both hemispheres at mid-high latitudes, which recalls the Doppler signature of an equator-to-pole circulation peaking at mid-latitudes (Fig. 2, bottom panel).

We then tested a second model composed of a solid-body zonal and an equator-to-pole meridional circulation pattern, i.e.,

$$ \Delta V_{\text{model 2}} = V_{z,0} \cos \lambda + V_{m,0} |\sin 2\lambda| + \kappa, \qquad (9) $$

where $V_{m,0}$ is the speed of the meridional circulation at $\lambda = 45^\circ$. The results of model 2 are shown on the second row of Fig. 15. The residual map and reduced $\chi^2$ (3.31) indicate a better agreement of the model with respect to the data, but a redshift bulb located in the equatorial region indicates that there are more terms to include in the model. Model 3 includes a variation of the zonal circulation as a function of longitude with either a maximum or a minimum located at the subsolar longitude (local noon meridian), i.e.,

$$ \Delta V_{\text{model 3}} = (V_{z,0} - V_{z,1} \sin \phi) \cos \lambda + V_{m,0} |\sin 2\lambda| + \kappa, \qquad (10) $$

where $\phi$ is the longitude, ranging from 0 to $\pi$ from the evening to the morning terminators, and $V_{z,1}$ is the amplitude of the noon-centered departure to a uniform solid-body rotation. The fit does
<figure><img src="image_15.png" /><figcaption>Fig. 14. Fitting the global circulation parameters with the MCMC routine. Top panel: models 4; bottom: model 7, both obtained from 250 000 iterations. Color scales and subpanels are the same as in Fig. 8.</figcaption></figure>

not significantly improve with respect to the previous fit (red. $\chi^2 = 2.90$). To make the model agree better, we add the option that the longitudinal variation is concentrated at low latitudes, i.e., within $\pm 20/30^\circ$ around the equator. Model 4 is built upon this assumption, by including a departure to zonal wind that is modulated by a Gaussian function centered around the equator

$$ \Delta V_{\text{model 4}} = V_{z,0} \cos \lambda - V_{z,1} \sin \phi e^{-\lambda^2/(2\sigma_{z,1})^2} + V_{m,0} |\sin 2\lambda| + \kappa, \quad (11) $$

where $V_{z,0}$ is the solid-body component at equator (i.e., the equatorial velocity at dawn for this specific model), $V_{z,1}$ is the amplitude of the “equatorial noon-centered” departure from a solid-body zonal circulation, and $\sigma_{z,1}$ is the standard deviation of the Gaussian function. Residuals are still present but are much less significant; i.e., $\approx 90\%$ of the points within $2\sigma$ with respect to the model. This is illustrated by a much lower reduced $\chi^2$ (1.42) and the cut along the equator, which shows a satisfactory agreement between the data and model. In Fig. 16 we represent

the map of the zonal wind<fnref n="5" /> deduced from the model 4, where it appears that zonal winds are larger than 200 m s⁻¹ in the morning – and evening, by extrapolating – and get slower in the subsolar area, down to 70 m s⁻¹. These values are not considered as accurate values of the wind speed at noon or in the morning, but as the result of a global fit of the atmospheric circulation of Venus in which global trends are retrieved but not local wind speed. The main fact here is that zonal winds are faster by a factor 2 to 3 in the morning and evening with respect to the subsolar region. In the following, we refer to model 4 as a circulation pattern in which global zonal winds display a hot spot structure – recalling some hot Jupiter’s nomenclature – in which winds are slower in the region where solar heating is maximum.

That being said, the difficulty is to not overinterpret data that likely suffer from biases and were obtained in nonoptimal conditions. Is it possible to go further without interpreting noise? We make other three attempts to refine the atmospheric circulation model (Fig. 17). The first (model 5) lets the hot spot free to shift along the equator, instead of being centered around noon, i.e.,

$$ \Delta V_{\text{model 5}} = V_{z,0} \cos \lambda - (V_{z,1,\sin} \sin \phi + V_{z,1,\cos} \cos \phi) e^{-\lambda^2/(2\sigma_{z,1})^2} + V_{m,0} |\sin 2\lambda| + \kappa, \quad (12) $$

where $V_{z,0}$ is the solid-body component at equator, $V_{z,1,\sin}$ and $V_{z,1,\cos}$ are the amplitude of the sin and cos component of the departure to pure solid-body zonal wind. The result is that the spot shifts of $\delta\phi \sim 7^\circ$ away for the subsolar meridian (noon), which is not significant given the noise<fnref n="6" />. Model 6 is identical to model 4 (hot spot is centered around noon) with the exception that both zonal and meridional flows may be periodic functions of the longitude, i.e.,

$$ \Delta V_{\text{model 6}} = V_{z,0} \cos \lambda - V_{z,1} \sin \phi e^{-\lambda^2/(2\sigma_{z,1})^2} + V_{m,0} |\sin 2\lambda| + V_{m,1} \sin \phi + \kappa, \quad (13) $$

where $V_{m,1}$ is the amplitude of the noon-meridional departure to uniform Hadley-cell meridional circulation. Unfortunately, in both cases the residual maps and reduced $\chi^2$ (1.41 and 1.51, respectively) do not show any significant difference with respect to model 4; this tends to indicate that model 4 reached the best possible model of this RV map and that going beyond is not reliable.

An option to interpret these data a little further consists of identifying large-scale patterns among the residuals of model 4. We may note, by eye, that to the west of the bisector meridian, the residuals are mostly red in the northern hemisphere and blue in the southern, and vice versa to the east of the meridian. Interpreting the “quadrupolar” appearance of the residuals out of the equatorial region could be caused by an asymmetric zonal wind circulation, faster in the southern hemisphere. Model 7 includes the possibility of a faster zonal circulation in either hemisphere as follows:

$$ \Delta V_{\text{model 7}} = V_{z,0} \cos \lambda - V_{z,1} \sin \phi e^{-\lambda^2/(2\sigma_{z,1})^2} - V_{z,2} \cos \lambda \sin (2\lambda) + V_{m,0} |\sin 2\lambda| + \kappa, \quad (14) $$

where the term $V_{z,2} \cos \lambda \sin (2\lambda)$ is a departure to symmetrical solid-body rotation. The term $\sin (2\lambda)$ involves that the NS zonal wind asymmetry is smooth (0 at equator and maximum

<footnote n="5">⁵ We do not make a similar plot for meridional circulation as it is a simple equator to pole regime (sine curve from equator to pole), which is uniform as a function of the longitude.</footnote>

<footnote n="6">⁶ $\delta\phi = \tan^{-1}\left(\frac{V_{z,1,\cos}}{V_{z,1,\sin}}\right)$.</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Model vs. data for models 1 to 4 for data taken on September 14, 2009. **Left panels:** best-fit model. **Second column of panels:** observation minus model, i.e., residual maps. **Third column of panels:** residual maps in terms of measurement errors. White color indicates agreement within $1\sigma$, light gray in between 1 and $2\sigma$, dark gray in between 2 and $3\sigma$, and black above $3\sigma$. In the four corners of each error map are indicated the fractions of points in each group ($[0-1]\sigma$, $[1-2]\sigma$, $[2-3]\sigma$, or $> 3\sigma$). **Right panels:** best fit models overplotted with data along the equator. Error bars are $1\sigma$ errors.</figcaption></figure>

at 45°). Contrarily to models 5 and 6, model 7 improves the match between data and fit (red. $\chi^2 = 1.08$, Fig. 18, Table 2). According to the results shown in Table 2, the zonal wind would be slightly larger in the northern hemisphere by about $18 \times \cos(\pi/4) \approx 13 \text{ m s}^{-1}$ at mid latitudes.

## 6. Discussion and prospects

This paper presents the first complete Doppler snapshots of a planet in the visible domain. The map obtained on September 14, 2009 is the result of integrating eight hours of data, which represents about 10% of the rotation period at the cloud-top altitude<fnref n="7" />. Despite poor observing conditions – even on September 14 the average seeing was 3 arcsec – we note a strong coherence of the RV maps (Fig. 13). Beyond biases on RVs, Doppler shifts vary

on the same range from side to side (color scale are the same in Fig. 13), indicating a clear retrograde rotation.

From a technical point of view, this paper makes use of an innovative method to measure and analyze the RVs of planetary atmosphere, following the recipes introduced by Gaulme et al. (2018) and applied in parallel to this work to observations of Jupiter (Gonçalves et al. 2019). The main improvements consist of using the correct expression of the Young effect, and including the atmospheric seeing both to locate (a posteriori) the spectrometer on the planet and to model RVs. Besides, we perform forward modeling to extract the wind components, by simultaneously fitting different components of the atmospheric circulation. The development of a dedicated MCMC routine ensures a good exploration of the parameter space and a proper estimate of error bars.

The first main result confirms what was expected from both cloud tracking and recent spectroscopic observations: solid body rotation alone is not sufficient to model observations and

<footnote n="7">⁷ The maps obtained on September 16 and 17, 2009 are the result of integrating 2 and 3.5 h of data, respectively, which represent about 2.5 and 4.5% of the rotation period at the cloud-top altitude.</footnote>
<table><caption>**Table 2.** Global circulation models.</caption><thead><tr><th></th><th colspan="3">Model 1</th><th colspan="2">Red. χ<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="2">Pure zonal (solid-body)</td><td>V<sub>z,0</sub> cos λ</td><td>Zonal amplitude</td><td>V<sub>z,0</sub></td><td>157.6<sup>+0.3</sup><sub>-2.3</sub></td><td rowspan="2">6.14</td></tr><tr><td>+κ</td><td>Offset</td><td>κ</td><td>5.5<sup>+0.7</sup><sub>-1.0</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 2</td></tr><tr><td rowspan="3">Zonal (solid-body)<br>and Meridional</td><td>V<sub>z,0</sub> cos λ</td><td>Zonal amplitude</td><td>V<sub>z,0</sub></td><td>156.9<sup>+1.2</sup><sub>-1.4</sub></td><td rowspan="3">3.31</td></tr><tr><td>+V<sub>m,0</sub> |sin λ|</td><td>Meridional</td><td>V<sub>m,0</sub></td><td>57.9<sup>+2.4</sup><sub>-2.7</sub></td></tr><tr><td>+κ</td><td>Offset</td><td>κ</td><td>-25.3<sup>+1.4</sup><sub>-1.8</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 3</td></tr><tr><td rowspan="4">Zonal (solid-body),<br>Noon zonal variation<br>and Meridional</td><td>V<sub>z,0</sub> cos λ</td><td>Zonal amplitude</td><td>V<sub>z,0</sub></td><td>211.8<sup>+7.8</sup><sub>-6.4</sub></td><td rowspan="4">2.90</td></tr><tr><td>-V<sub>z,1</sub> sin φ</td><td>Departure to zonal</td><td>V<sub>z,1</sub></td><td>79.9<sup>+12.0</sup><sub>-8.4</sub></td></tr><tr><td>+V<sub>m,0</sub> |sin λ|</td><td>Meridional</td><td>V<sub>m,0</sub></td><td>47.9<sup>+2.8</sup><sub>-2.7</sub></td></tr><tr><td>+κ</td><td>Offset</td><td>κ</td><td>-10.4<sup>+2.3</sup><sub>-2.6</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 4</td></tr><tr><td rowspan="4">Zonal (solid-body),<br>Equatorial noon zonal variation<br>and Meridional</td><td>V<sub>z,0</sub> cos λ</td><td>SB zonal amplitude</td><td>V<sub>z,0</sub></td><td>214.4<sup>+5.1</sup><sub>-5.8</sub></td><td rowspan="4">1.42</td></tr><tr><td>-V<sub>z,1</sub> sin φ e<sup>-λ<sup>2</sup>/(2σ<sub>z,1</sub>)<sup>2</sup></sup></td><td>Departure to zonal</td><td>V<sub>z,1</sub></td><td>143.2<sup>+6.8</sup><sub>-8.3</sub></td></tr><tr><td>+V<sub>m,0</sub> |sin λ|</td><td>Sigma of the Gaussian</td><td>σ<sub>z,1</sub></td><td>16.7<sup>+1.3</sup><sub>-1.7</sub></td></tr><tr><td>+κ</td><td>Offset</td><td>V<sub>m,0</sub></td><td>44.9<sup>+2.5</sup><sub>-3.0</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 5</td></tr><tr><td rowspan="4">Zonal (solid-body),<br>Equatorial decentered zonal variation<br>and Meridional</td><td>V<sub>z,0</sub> cos λ</td><td>SB zonal amplitude</td><td>V<sub>z,0</sub></td><td>235.4<sup>+16.3</sup><sub>-12.9</sub></td><td rowspan="4">1.41</td></tr><tr><td>-V<sub>z,1,sin</sub> sin φ e<sup>-λ<sup>2</sup>/(2σ<sub>z,1</sub>)<sup>2</sup></sup></td><td>Departure to zonal</td><td>V<sub>z,1,sin</sub></td><td>157.4<sup>+15.5</sup><sub>-13.2</sub></td></tr><tr><td>-V<sub>z,1,cos</sub> cos φ e<sup>-λ<sup>2</sup>/(2σ<sub>z,1</sub>)<sup>2</sup></sup></td><td>Departure to zonal</td><td>V<sub>z,1,cos</sub></td><td>20.2<sup>+2.7</sup><sub>-2.4</sub></td></tr><tr><td>+V<sub>m,0</sub> |sin λ|</td><td>Sigma of the Gaussian</td><td>σ<sub>z,1</sub></td><td>20.2<sup>+2.7</sup><sub>-2.4</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 6</td></tr><tr><td rowspan="4">Zonal (solid-body),<br>Equatorial noon zonal variation,<br>and Meridional including<br>longitudinal variation (noon centered)</td><td>V<sub>z,0</sub> cos λ</td><td>SB zonal amplitude</td><td>V<sub>z,0</sub></td><td>218.6<sup>+6.7</sup><sub>-11.5</sub></td><td rowspan="4">1.51</td></tr><tr><td>-V<sub>z,1</sub> sin φ e<sup>-λ<sup>2</sup>/(2σ<sub>z,1</sub>)<sup>2</sup></sup></td><td>Departure to zonal</td><td>V<sub>z,1</sub></td><td>149.9<sup>+9.2</sup><sub>-17.2</sub></td></tr><tr><td>+V<sub>m,0</sub> |sin λ|</td><td>Sigma of the Gaussian</td><td>σ<sub>z,1</sub></td><td>16.8<sup>+1.2</sup><sub>-1.8</sub></td></tr><tr><td>+V<sub>m,1</sub> sin φ</td><td>Meridional</td><td>V<sub>m,0</sub></td><td>49.2<sup>+22.9</sup><sub>-19.6</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 7</td></tr><tr><td rowspan="4">Zonal (solid-body),<br>Equatorial noon zonal variation<br>N-S zonal asymmetry<br>and Meridional</td><td>V<sub>z,0</sub> cos λ</td><td>SB zonal amplitude</td><td>V<sub>z,0</sub></td><td>215.4<sup>+6.0</sup><sub>-5.5</sub></td><td rowspan="4">1.08</td></tr><tr><td>-V<sub>z,1</sub> sin φ e<sup>-λ<sup>2</sup>/(2σ<sub>ENZ</sub>)<sup>2</sup></sup></td><td>Departure to zonal</td><td>V<sub>z,1</sub></td><td>144.4<sup>+6.9</sup><sub>-8.5</sub></td></tr><tr><td>-V<sub>z,2</sub> cos λ sin (2λ)</td><td>Sigma of the Gaussian</td><td>σ<sub>z,1</sub></td><td>17.1<sup>+1.4</sup><sub>-1.6</sub></td></tr><tr><td>+V<sub>m,0</sub> |sin λ|</td><td>Amplitude NS asymmetry</td><td>V<sub>z,2</sub></td><td>18.0<sup>+2.1</sup><sub>-2.7</sub></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Model 8</td></tr><tr><td></td><td>+κ</td><td>Offset</td><td>κ</td><td>-7.6<sup>+1.9</sup><sub>-3.0</sub></td><td></td></tr></tbody></table>

**Notes.** The value $V_{z,0}$ is the amplitude of zonal winds along the equator, $V_{z,1}$ is the amplitude of the solid-body (SB) zonal component for models that have longitudinal variations of the zonal circulation (models 4 to 7), $V_{m,0}$ is the amplitude of meridional winds at latitude $λ = 45°$, and $κ$ is a global RV offset. Model 3: $V_{z,1}$ (“noon zonal” component) is the amplitude of departure to solid body rotation which is minimum at noon. Models 4, 6, and 7: $V_{z,1}$ is the amplitude of the equatorial noon-centered zonal component. Model 5: $V_{z,1,sin}$ and $V_{z,1,cos}$ are the amplitudes of the equatorial decentered zonal components for the sine and cosine, respectively. Model 6: $V_{z,1}$ is the amplitude of the noon-centered meridional component. All terms are in m s$^{-1}$ except $σ_{z,1}$, which is in degrees.
<figure><img src="image_17.png" /><figcaption>Fig. 16. Zonal wind map that corresponds to Model 4.</figcaption></figure>

equator-to-pole meridional circulation is needed. However, we find the speed of zonal and meridional winds to be larger than previously measured. It is hard to compare the zonal wind values because we identify a strong longitudinal and latitudinal variation, however, no observations have indicated winds as large as 200 m s⁻¹ at the morning terminator in the equatorial region, so far. As regards meridional winds, our model is identical to that used to interpret previous observations (e.g., Machado et al. 2017) and we find an amplitude about twice as large (about 45 instead of 22 m s⁻¹ at mid latitudes).

The speeds of the winds that we report are significantly higher than those that have been reported so far from both cloud tracking and spectroscopic measurements. We first note that no study has ever produced a complete RV map of Venus in the visible, which makes a direct comparison hard to perform. However, we now explore how compatible our results are with respect to those from the literature, in particular the most recent results by Machado et al. (2012, 2014, 2017). For this, we fit the RV map obtained on September 14, 2009 with a simple zonal plus meridional wind without taking account the biases on RVs caused by the atmospheric seeing. The model consists of an equator-to-pole meridional circulation (∝ |sin 2λ|) and a zonal wind whose amplitude is fitted for each band of latitude. The zonal wind profile as a function latitude obtained that way is shown in Fig. 19. We observe a (noisy) M-shaped profile where zonal winds peak at mid-latitudes with an amplitude of about 120 m s⁻¹ and a local minimum around the equator at about 100 m s⁻¹, which

is very similar to what is reported in the literature. This suggests that if atmospheric seeing and instrumental PSF were included in the previous high-resolution spectroscopic works, a larger zonal wind would have been measured. We are aware that reported seeing and pointing accuracy by Machado et al. (2014, 2017) are below one arcsec each. However, the diameter of the spectrometer's fiber is 1.6 arcsec, which implies that combined instrumental plus atmospheric PSF is larger than 2 arcsec and that their RV measurements are likely to be biased too, especially toward the planetary edges. Interestingly, we note that wind measurements of Venus by Machado et al. (2012, 2014, 2017) provide very consistent results, with average zonal winds estimated in between 117 and 123 m s⁻¹, whereas Widemann et al. (2007) and Gabsi et al. (2008) reported much lower values (≈75 m s⁻¹). The earlier 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) reported rather poor observing conditions at the Observatoire de Haute Provence, which could explain in part why they obtain such low wind estimates.

The second main result is the hot spot structure of the atmospheric circulation, at least for the zonal component, which had never been suggested that clearly so far. Cloud tracking measurements and Doppler spectroscopic measurements indicate possible longitudinal variations of the wind as functions of local time with faster circulation toward the terminator (e.g., Khatuntsev et al. 2013; Hueso et al. 2015; Machado et al. 2017). In particular, Khatuntsev et al. (2013) found diurnal variations of the zonal winds but at the ~10 m s⁻¹ level only for local times 6–18 h. Since then, based on VIRTIS VEx observations, Hueso et al. (2015) reported zonal winds variations as a function of longitude, from 90 m s⁻¹ slightly before noon (local time) up to 130 m s⁻¹ at 17 h for latitudes of about -30°. Machado et al. (2017) also reported zonal wind variations at [20, 40]° latitudes, ranging from 125 at ≈10 h to 167 m s⁻¹ at ~7 h. We also note that their local values of zonal winds are a little larger in the northern hemisphere, as we find in this work. However, no such hot spot pattern has been identified. In contrast, from clouds tracking measurements performed on images from the VEx VMC, Bertaux et al. (2016) did not find zonal wind speed variation as a function of local time, but as a function of geographical features with values ranging from 101 to 83 m s⁻¹ at latitudes in between 5 and 15°S.

We note that a hot spot regime of zonal winds is compatible with an M-shaped mean zonal wind profile as a function of latitude. If we average the zonal winds at a given latitude, we obtain that type of profile, with no wind at poles, a maximum amplitude at mid-latitudes and a local minimum in the equatorial region. Figure 19 also shows the zonal profile obtained by averaging the zonal wind along longitudes on the visible dayside of Venus: it shows that the observed M-shape can be the result of a hot spot structure. Such a regime recalls global circulation models that are proposed to explain light curves of stars hosting hot Jupiters. In particular, Showman & Guillot (2002) studied the atmospheric circulation of hot Jupiters that are locked to their host stars. Venus can be seen to be somewhat similar as its atmosphere is very dense (90 bars on the ground) and its rotation period (243 days) is very close to its orbital period (225 days). A significant difference from the case studied by Showman & Guillot (2002) regards the global circulation regime, which is assumed to be in geostrophic regime for hot Jupiters, and which is considered to be cyclostrophic in the case of planets with little planetary rotation such as Venus or Titan. Detailed numerical
<figure><img src="image_18.png" /><figcaption>Fig. 17. Same as Fig. 15 for models 5–7.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Zonal wind map that corresponds to Model 7.</figcaption></figure>

<caption>Fig. 19. Plain black line: fitting of the RV map obtained from the observations done on September 14, 2009 without taking account of the atmospheric seeing on RV profiles. Gray lines: mean zonal wind on the visible part of the dayside of Venus as computed from model 7 based on the observations done on September 14, 2009. In other words, it is the mean of the map shown in Fig 18, bottom panel, with 100 executions to take into account the errors on the parameters of the model. The black dashed line indicates the average mean profile of model 7.</caption>

simulations are needed to investigate this question further, which is beyond the scope of our paper.

Our two main results, i.e., both the zonal and meridonal large wind speeds as well as the hot spot pattern of the zonal flow, are significantly different with respect to what has been measured so far. They will likely remain controversial until further new independent observations, either by another team or with another instrument. Indeed, the presence of uncontrolled
and not fully understood biases (see Sect 4.3) makes the conclusion partially questionable. This works points out the difficulty of constructing RV maps from individual spectra, even with a long-slit spectrometer that covers the whole planetary diameter at once. The constantly changing observation conditions and pointing stability of the telescope, as well as the wavelength stability of the spectrometer, are major challenges to conduct clean measurements.

There are two options for go further. The first consists of using the same type of instrumentation (long-slit spectrometry) with adaptive optics (AO) to reduce the seeing issue, which has a huge impact on map reconstruction. At the same time, even with AO, these kinds of observations are dependent on thermal drifts of the spectrometer and on clouds during observations. The THEMIS telescope plans to get an AO system soon, however it will be designed for working on the Sun, which is likely not adapted for observing planetary targets. In the short term, the most promising alternative is the JOVIAL/JIVE project (PI: F.-X. Schmider), which consists of three identical Doppler imaging spectrometers that will be installed in France, New Mexico, and Japan (Gonçalves et al. 2016, 2019). The original purpose of the project is to get RV maps of Jupiter and Saturn at short cadence (30 s) to conduct seismic observations of these planets similar to techniques of helioseismology. This second objective of this instrument is to measure the atmospheric dynamics of the dense atmospheres of the solar system, i.e., mostly Venus, Jupiter, and Saturn, given that Uranus, Neptune, and Titan are faint and barely resolved from the ground. The JOVIAL/JIVE instrument is able to provide a complete RV map in only one exposure and does not have any of the issues we encountered in this work in building maps from 1D segments. A paper dedicated to the first Doppler measurements of the wind of Jupiter was recently published by Gonçalves et al. (2019). Preliminary observations of Venus were led in spring 2018, and more observations could be done with AO with the C2PU 1 m telescope at the Calern observatory (France) in 2019 or 2020. We will first look at whether the hot spot structure is real with the help of high spatial resolution. The second exciting opportunity will be to make use of the three instruments mentioned above in order to observe Venus continuously for several days to study time variations of the atmospheric circulation of the planet.

*Acknowledgements*. P.G. was supported in part by the German space agency (Deutsches Zentrum für Luft- und Raumfahrt) under PLATO data grant 500O1501. I.G.'s PhD was granted by Observatoire de la Côte d'Azur and the JIVE in NM project (NASA EPSCoR grant #NNX14AN67A). THEMIS is a UPS of the CNRS funded by INSU. The authors are grateful to the referee, E. Lellouch, for his intensive review of the paper, which contributed to improving its quality.

## References

Bertaux, J.-L., Khatuntsev, I. V., Hauchecorne, A., et al. 2016, J. Geophys. Res. Planets, 121, 1087

Boyer, C., & Guérin, P. 1969, Icarus, 11, 338

Civeit, T., Appourchaux, T., Lebreton, J.-P., et al. 2005, A&A, 431, 1157

Clancy, R. T., Sandor, B. J., & Moriarty-Schieven, G. H. 2008, Planet. Space Sci., 56, 1344

Clancy, R. T., Sandor, B. J., & Moriarty-Schieven, G. 2012, Icarus, 217, 794

Connes, P. 1985, Ap&SS, 110, 211

Donati, J.-F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658

Drossart, P., Piccioni, G., Gérard, J. C., et al. 2007, Nature, 450, 641

Gabsi, Y., Bertaux, J. L., Hauchecorne, A., Schmitt, J., & Guibert, S. 2008, Planet. Space Sci., 56, 1454

Gaulme, P., Schmider, F.-X., Grec, C., et al. 2008, Planet. Space Sci., 56, 1335

Gaulme, P., Schmider, F.-X., & Gonçalves, I. 2018, A&A, 617, A41

Gierasch, P. J., Goody, R. M., Young, R. E., et al. 1997, in Venus II: Geology, Geophysics, Atmosphere, and Solar Wind Environment, eds. S. W. Bougher, D. M. Hunten, & R. J. Phillips (Tucson, AZ: University of Arizona Press), 459

Gonçalves, I., Schmider, F.-X., Bresson, Y., et al. 2016, in Ground-based and Airborne Instrumentation for Astronomy VI, Proc. SPIE, 9908, 99083M

Gonçalves, I., Schmider, F. X., Gaulme, P., et al. 2019, Icarus, 319, 795

Hueso, R., Peralta, J., & Sánchez-Lavega, A. 2012, Icarus, 217, 585

Hueso, R., Peralta, J., Garate-Lopez, I., Bandos, T. V., & Sánchez-Lavega, A. 2015, Planet. Space Sci., 113, 78

Khatuntsev, I. V., Patsaeva, M. V., Titov, D. V., et al. 2013, Icarus, 226, 140

Lebonnois, S., Hourdin, F., Eymet, V., et al. 2010, J. Geophys. Res. Planets, 115, E06006

Lellouch, E., & Witasse, O. 2008, Planet. Space Sci., 56, 1317

Lellouch, E., Goldstein, J. J., Rosenqvist, J., Bougher, S. W., & Paubert, G. 1994, Icarus, 110, 315

Lellouch, E., Clancy, T., Crisp, D., et al. 1997, in Venus II: Geology, Geophysics, Atmosphere, and Solar Wind Environment, eds. S. W. Bougher, D. M. Hunten, & R. J. Phillips (Tucson, AZ: University of Arizona Press), 295

Lellouch, E., Paubert, G., Moreno, R., & Moullet, A. 2008, Planet. Space Sci., 56, 1355

Limaye, S. S., & Rengel, M. 2013, Atmospheric Circulation and Dynamics, eds. L. Bengtsson, R.-M. Bonnet, D. Grinspoon, et al. (Heidelberg: Springer), 55

Limaye, S. S., Grassotti, C., & Kuetemeyer, M. J. 1988, Icarus, 73, 193

Luz, D., Civeit, T., Courtin, R., et al. 2005, Icarus, 179, 497

Luz, D., Civeit, T., Courtin, R., et al. 2006, J. Geophys. Res. Planets, 111, E08S90

Machado, P., Luz, D., Widemann, T., Lellouch, E., & Witasse, O. 2012, Icarus, 221, 248

Machado, P., Widemann, T., Luz, D., & Peralta, J. 2014, Icarus, 243, 249

Machado, P., Widemann, T., Peralta, J., et al. 2017, Icarus, 285, 8

Markiewicz, W. J., Titov, D. V., Limaye, S. S., et al. 2007, Nature, 450, 633

Mein, P., & Rayrole, J. 1985, Vistas in Astronomy (Amsterdam: Elsevier), 28, 567

Moissl, R., Khatuntsev, I., Limaye, S. S., et al. 2009, J. Geophys. Res. Planets, 114, 0

Moreux, T. A. 1928, Le Ciel et l'Univers, ed. G. Doin & Cie (Paris: Librairie Octave Doin, 101

Moullet, A., Lellouch, E., Moreno, R., Gurwell, M., & Sagawa, H. 2012, A&A, 546, A102

Nakamura, M., Imamura, T., Ishii, N., et al. 2016, Earth, Planets, and Space, 68, 75

Patsaeva, M. V., Khatuntsev, I. V., Patsaev, D. V., et al. 2015, Planet. Space Sci., 113, 100

Peralta, J., Hueso, R., Sánchez-Lavega, A., et al. 2008, J. Geophys. Res. Planets, 113, 0

Peralta, J., Luz, D., Berry, D. L., et al. 2012, Icarus, 220, 958

Rossow, W. B., del Genio, A. D., & Eichler, T. 1990, J. Atmos. Sci., 47, 2053

Sánchez-Lavega, A., Hueso, R., Piccioni, G., et al. 2008, Geophys. Res. Lett., 35, 13204

Sánchez-Lavega, A., Lebonnois, S., Imamura, T., Read, P., & Luz, D. 2017, Space Sci. Rev., 212, 1541

Showman, A. P., & Guillot, T. 2002, A&A, 385, 166

Sornig, M., Livengood, T., Sonnabend, G., et al. 2008, Planet. Space Sci., 56, 1399

Sornig, M., Livengood, T. A., Sonnabend, G., Stupar, D., & Kroetz, P. 2012, Icarus, 217, 863

Svedhem, H., Titov, D. V., Taylor, F. W., & Witasse, O. 2007, Nature, 450, 629

Takagi, M., & Matsuda, Y. 2007, J. Geophys. Res. Atmos., 112, D09112

Titov, D. V., Taylor, F. W., Svedhem, H., et al. 2008, Nature, 456, 620

Titov, D. V., Markiewicz, W. J., Ignatiev, N. I., et al. 2012, Icarus, 217, 682

Traub, W. A., & Carleton, N. P. 1975, J. Atmos. Sci., 32, 1045

Widemann, T., Lellouch, E., & Campargue, A. 2007, Planet. Space Sci., 55, 1741

Widemann, T., Lellouch, E., & Donati, J.-F. 2008, Planet. Space Sci., 56, 1320

Young, A. T. 1975, Icarus, 24, 1

Young, A. T., Schorn, R. A., Young, L. D. G., & Crisp, D. 1979, Icarus, 38, 435