# Amateur Observers Witness the Return of Venus’ Cloud Discontinuity

E. Kardasis Corresponding author: E. Kardasis (astromanos2002@yahoo.gr).
Hellenic Amateur Astronomy Association, Alopekis 42, 10676,Athens, Greece.  J.
Peralta  Departamento de FAMN, Facultad de Física, Universidad de Sevilla,
Sevilla 41012, Spain.  G. Maravelias  Hellenic Amateur Astronomy Association,
Alopekis 42, 10676,Athens, Greece.  IAASARS, National Observatory of
Athens,Lofos Nymfon, Thission, 118 10,Athens, Greece.  Institute of
Astrophysics, FORTH, Heraklion, Greece.  M. Imai  Faculty of Science, Kyoto
Sangyo University, Japan.  A. Wesley  Astronomical Society of Australia,
Rubyvale, QLD 4702, Australia.  T. Olivetti  Union of Italian Amateur
Astronomers, Via Lazio, 14-00040 Rocca di Papa, RM, Italy.  Y. Naryzhniy
Kagarlyk, Kiev Region, Ukraine.  L. Morrone  AstroCampania
Association,Agerola, Italy.  British Astronomical Association, Burlington
House, Piccadilly, London W1J 0DU, UK.  A. Gallardo  Asociación Astronómica
del Campo de Gibraltar “Luz Cero”, Spain.  G. Calapai  Private Astronomical
Observatory in Massa S. Giorgio, Messina, Italy.  J. Camarena  Agrupación
Astronómica de la Safor, Valencia, Spain.  P. Casquinha  Portuguese
Association of Amateur Astronomers, Portugal.  D. Kananovich  Faculty of
Chemistry,Tallinn University of Technology, Tallinn, Estonia.  N. MacNeill
British Astronomical Association, Burlington House, Piccadilly, London W1J
0DU, UK.  Astronomy Society of NSW, PO Box 870 Epping NSW 1710, Australia.  C.
Viladrich  Astroqueyras, France.  A. Takoudi  Hellenic Amateur Astronomy
Association, Alopekis 42, 10676,Athens, Greece.

###### Abstract

Firstly identified in images from JAXA’s orbiter Akatsuki , the cloud
discontinuity of Venus is a planetary-scale phenomenon known to be recurrent
since, at least, the 1980s. Interpreted as a new type of Kelvin wave, this
disruption is associated to dramatic changes in the clouds’ opacity and
distribution of aerosols, and it may constitute a critical piece for our
understanding of the thermal balance and atmospheric circulation of Venus.
Here, we report its reappearance on the dayside middle clouds four years after
its last detection with Akatsuki /IR1, and for the first time, we characterize
its main properties using exclusively near-infrared images from amateur
observations. In agreement with previous reports, the discontinuity exhibited
temporal variations in its zonal speed, orientation, length, and its effect
over the clouds’ albedo during the 2019/2020 eastern elongation. Finally, a
comparison with simultaneous observations by Akatsuki UVI and LIR confirmed
that the discontinuity is not visible on the upper clouds’ albedo or thermal
emission, while zonal speeds are slower than winds at the clouds’ top and
faster than at the middle clouds, evidencing that this Kelvin wave might be
transporting momentum up to upper clouds.

_K_ eywords terrestrial planets $\cdot$ Venus $\cdot$ atmosphere $\cdot$
atmospheric dynamics $\cdot$ atmospheric waves

##  1 Introduction

Venus is permanently covered by a thick stratified atmosphere. Within
altitudes $\sim$48–70 km above the surface, we find clouds of sulfuric acid
(Titov et al., 2018). These are divided into upper, middle, and lower clouds.
At the top (56.5–70 km above the surface), the clouds are dominated by a
retrograde superrotation that reaches speeds of
$\sim$100$\,\textrm{m}\cdot\textrm{s}^{-1}\,$at the clouds’ top of lower
latitudes, which is 60 times faster than the planetary rotation (Sánchez-
Lavega et al., 2017). This is the layer that has been observed continuously
for almost a century with violet and ultraviolet (UV) filters (Wright, 1927;
Ross, 1928), whose dark/bright features are caused by known (SO2) and unknown
absorbers (Titov et al., 2018; Pérez-Hoyos et al., 2018). Under this layer,
the middle clouds (50.5–56.5 km above the surface) are observed on the dayside
at visible and near-infrared (NIR) wavelengths up to 1 $\mu$m (Peralta et al.,
2019; Titov et al., 2018). The lower clouds (47.5–50.5 km) can be observed as
silhouettes partially blocking the deeper thermal emission on the nightside
mainly at the spectral windows 1.74 $\mu$m and 2.32 $\mu$m (Limaye et al.,
2018; McGouldrick et al., 2012; Titov et al., 2018). The middle and lower
clouds contribute to the greenhouse effect and the radiative energy balance
(Peralta et al., 2020). The morphology and motions of the clouds differ at
each of the three layers (Horinouchi et al., 2018; Limaye et al., 2018;
Peralta et al., 2019; Sánchez-Lavega et al., 2017; Titov et al., 2018).

Multiple observations have provided evidence of planetary-scale waves in this
stratified atmosphere, and these waves are suspected to play a critical role
in powering of the superrotation (Sánchez-Lavega et al., 2017). These have
been mainly observed at the upper clouds, such as the Y feature (Peralta et
al., 2015), the stationary bow-shaped wave generated at the surface (Fukuhara
et al., 2017), or the solar tides, which has been shown to maintain the
superrotation near the equator (Horinouchi et al., 2020). Recently, Peralta et
al. (2020) reported the discovery of a giant discontinuity/disruption at the
deeper clouds of Venus, which was shown to be a cyclical phenomenon limited
between 30∘N and 40∘S, which importantly alters clouds’ properties and
aerosols (Peralta et al., 2020; McGouldrick et al., 2021), propagates faster
than the superrotating winds, and is able to keep a coherent shape after
several revolutions. Based on numerical simulations, Peralta et al. (2020)
interpreted this phenomenon as a new type of planetary-scale wave that may be
feeding the superrotation at the upper clouds by bringing momentum from the
deeper atmosphere, which is precisely where the Venus atmosphere accumulates
most of its angular momentum (Schubert, 1983). Its absence at the top of the
clouds and range of phase speeds also support that this wave should dissipate
within the upper clouds (Peralta et al., 2020). Peralta et al. (2020) also
showed that the cloud discontinuity (CD henceforth) has been a recurrent
feature for decades at the nightside lower clouds, although it is less
frequent at the dayside middle clouds where it was observed for the last time
on 27 November 2016 by the camera Akatsuki /IR1, before this instrument and
IR2 stopped working indefinitely (Iwagami et al., 2018).

Here, we report the reappearance of the CD on the dayside middle clouds during
the 2019/2020 eastern elongation of Venus, using observations from small
private telescopes around the world coordinated with JAXA’s orbiter Akatsuki
(Nakamura et al., 2016). In Section 2, the outline of the observations and the
data reduction method are presented. The results of our analysis are provided
and discussed in Section 3. Finally, Section 4, summarizes our work.

##  2 Observations and Methods

###  2.1 Venus Images Taken with Ground-Based Telescopes

To study the CD events, we inspected the amateur observations of Venus hosted
in the database of ALPO-Japan http://alpo-j.sakura.ne.jp/indexE.htm (accessed
on 2021) (Sato, 2018). We selected and examined 458 NIR and 552 UV
observations of the eastern elongation 2019/2020. The NIR observations were
obtained with filters in the range 685–1050 nm and the UV observations with
filters in the $\sim$320–410 nm range. Table 1 provides details about these
NIR and UV observations as well as the Venus phase angle and the nominal best
resolution. The resolution of the amateur images ranged from 187 to 663 km
depending on the apparent angular size of Venus and the best resolution
allowed by the diameter of the telescope (Sánchez-Lavega et al., 2016). One of
us (EK), on 11 March 2020 (with Venus having a solar elongation of 45.6∘ East,
$20.6^{\prime\prime}$ in diameter, and 58% of illuminated disk), identified a
strong manifestation of a CD and alerted Venus observers and researchers. That
event stimulated intense monitoring and led to the current work. We inspected
a period covering from 13 October 2019 to 25 April 2020 (Solar Elongation
16.2∘–40.8∘ East, angular size 10.3′′–35.8′′, illuminated fraction
0.962–0.295).

Table 1: Summary of NIR and UV images of Venus used in this work. We present the date and time of each observation (column 1), the Venus phase angle (column 2), the telescope diameter and type (column 3), the best spatial resolution (column 4), which was calculated theoretically from telescope’s properties and apparent disk diameter, the effective transmission of the equipment used (column 5), the integration time (column 6), which corresponds to the total exposure time for each image, the observer (column 7), and country (column 8). Date and Time (UT) | Venus Phase | Telescope | Best Spatial | Effective | Integration | Observer | Country
---|---|---|---|---|---|---|---
(YYYY-MM-DDThh:mm.m) | Angle (deg) | Diameter (m) | Resolution (km) | Transmission (nm) | Time (s) |  |
2019-11-14T07:22.7 | 33.3 | 0.41 (N) | 589 | 750–1020 | 180 | Anthony Wesley | Australia
2019-11-17T07:34.9 | 34.3 | 0.41 (N) | 589 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-18T07:53.2 | 34.6 | 0.41 (N) | 583 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-19T07:50.3 | 35.0 | 0.41 (N) | 583 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-20T07:50.5 | 35.3 | 0.41 (N) | 583 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-22T07:53.6 | 36.0 | 0.41 (N) | 578 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-23T08:13.9 | 36.3 | 0.41 (N) | 578 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-23T15:00.6 | 36.3 | 0.35 (S) | 663 | 884–900 | 300 | Emmanuel Kardasis | Greece
2019-11-24T08:07.9 | 36.7 | 0.41 (N) | 573 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-27T14:19.0 | 37.7 | 0.35 (S) | 652 | 884–900 | 300 | Emmanuel Kardasis | Greece
2019-11-29T07:58.2 | 38.4 | 0.41 (N) | 568 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-11-30T08:06.5 | 38.7 | 0.41 (N) | 563 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-12-01T14:41.0 | 39.1 | 0.35 (S) | 646 | 884–900 | 300 | Emmanuel Kardasis | Greece
2019-12-02T08:09.8 | 39.4 | 0.41 (N) | 558 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-12-03T15:20.0 | 39.7 | 0.35 (S) | 641 | 884–900 | 300 | Emmanuel Kardasis | Greece
2019-12-04T08:13.4 | 40.1 | 0.41 (N) | 558 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-12-05T08:02.4 | 40.4 | 0.41 (N) | 554 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-12-06T08:48.0 | 40.8 | 0.35 (S) | 635 | 850–1020 | 180 | Niall MacNeill | Australia
2019-12-07T14:55.0 | 41.1 | 0.35 (S) | 630 | 884–900 | 300 | Emmanuel Kardasis | Greece
2019-12-08T15:29.0 | 41.5 | 0.35 (S) | 630 | 884–900 | 300 | Emmanuel Kardasis | Greece
2019-12-18T08:09.0 | 45.0 | 0.41 (N) | 527 | 1000–1020 | 180 | Anthony Wesley | Australia
2019-12-18T15:23.0 | 45.0 | 0.35 (S) | 604 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-01-11T11:34.0 | 53.7 | 0.50 (D) | 393 | 807–1050 | 249 | Tiziano Olivetti | Thailand
2020-01-11T14:51.0 | 53.7 | 0.35 (S) | 547 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-01-29T11:42.2 | 60.8 | 0.50 (D) | 356 | 807–1050 | 315 | Tiziano Olivetti | Thailand
2020-02-01T11:27.6 | 62.0 | 0.50 (D) | 352 | 807–1050 | 415 | Tiziano Olivetti | Thailand

Table 2: _Cont_. Date and Time (UT) | Venus Phase | Telescope | Best Spatial | Effective | Integration | Observer | Country
---|---|---|---|---|---|---|---
(YYYY-MM-DDThh:mm.m) | Angle (deg) | Diameter (m) | Resolution (km) | Transmission (nm) | Time (s) |  |
2020-02-01T15:51.0 | 62.0 | 0.35 (S) | 490 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-02-05T11:38.2 | 63.7 | 0.50 (D) | 343 | 807–1050 | 415 | Tiziano Olivetti | Thailand
2020-02-16T07:55.0 | 68.6 | 0.41 (N) | 387 | 1000–1020 | 180 | Anthony Wesley | Australia
2020-02-24T17:27.0 | 72.3 | 0.35 (S) | 419 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-02-24T18:48.5 | 72.3 | 0.35 (S) | 419 | 807–1050 | 380 | Paulo Casquinha | Portugal
2020-02-28T17:41.7 | 74.2 | 0.28 (S) | 509 | 742–1050 | 240 | Antonio Gallardo | Spain
2020-03-07T15:59.0 | 78.3 | 0.35 (S) | 377 | 850–1020 | 160 | Raimondo Sedrani | Italy
2020-03-08T11:37.0 | 78.9 | 0.50 (D) | 269 | 807–1050 | na | Tiziano Olivetti | Thailand
2020-03-09T11:32.0 | 79.4 | 0.50 (D) | 266 | 807–1050 | na | Tiziano Olivetti | Thailand
2020-03-10T16:57.0 | 79.9 | 0.35 (S) | 367 | 850–1020 | 160 | Raimondo Sedrani | Italy
2020-03-11T16:44.8 | 80.5 | 0.35 (S) | 366 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-03-11T17:09.1 | 80.5 | 0.35 (S) | 366 | 1000–1050 | 300 | Luigi Morrone | Italy
2020-03-11T16:30.0 | 80.5 | 0.35 (S) | 366 | 807–1050 | 60 | Joaquin Camarena | Spain
2020-03-12T15:59.0 | 81.0 | 0.35 (S) | 360 | 807–1050 | 60 | Joaquin Camarena | Spain
2020-03-13T11:22.0 | 81.6 | 0.50 (D) | 256 | 807–1050 | 482 | Tiziano Olivetti | Thailand
2020-03-14T16:25.0 | 82.1 | 0.35 (S) | 354 | 807–1050 | 60 | Joaquin Camarena | Spain
2020-03-15T15:10.0 | 82.7 | 0.40 (D) | 314 | 950–1050 | 180 | Yaroslav Naryzhniy | Ukraine
2020-03-16T15:00.0 | 83.3 | 0.40 (D) | 313 | 950–1050 | 180 | Yaroslav Naryzhniy | Ukraine
2020-03-16T16:01.0 | 83.3 | 0.28 (S) | 438 | 685–1050 | 1000 | Giovanni Calapai | Italy
2020-03-16T16:54.4 | 83.3 | 0.35 (S) | 349 | 1000–1050 | 300 | Luigi Morrone | Italy
2020-03-21T11:34.4 | 86.2 | 0.50 (D) | 237 | 807–1050 | 301 | Tiziano Olivetti | Thailand
2020-03-21T17:35.0 | 86.2 | 0.25 (K) | 330 | 830–1020 | 205 | Dzmitry Kananovich | Estonia
2020-03-21T17:10.2 | 86.2 | 0.35 (S) | 330 | 850–1020 | 160 | Raimondo Sedrani | Italy
2020-03-21T17:35.0 | 86.2 | 0.35 (S) | 330 | 884–900 | 300 | Emmanuel Kardasis | Greece

Table 3: _Cont_.

Date and Time (UT) | Venus Phase | Telescope | Best Spatial | Effective | Integration | Observer | Country
---|---|---|---|---|---|---|---
(YYYY-MM-DDThh:mm.m) | Angle (deg) | Diameter (m) | Resolution (km) | Transmission (nm) | Time (s) |  |
2020-03-26T15:20.0 | 89.3 | 0.40 (D) | 281 | 950–1050 | 180 | Yaroslav Naryzhniy | Ukraine
2020-03-26T17:01.6 | 89.3 | 0.28 (S) | 394 | 742–1050 | 240 | Antonio Gallardo | Spain
2020-03-31T16:15.0 | 92.6 | 0.40 (D) | 266 | 950–1050 | 180 | Yaroslav Naryzhniy | Ukraine
2020-03-31T17:04.0 | 92.6 | 0.35 (S) | 296 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-04-05T16:40.0 | 96.1 | 0.40 (D) | 250 | 950–1050 | 180 | Yaroslav Naryzhniy | Ukraine
2020-04-10T16:18.0 | 99.9 | 0.35 (S) | 261 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-04-10T17:00.1 | 99.9 | 0.35 (S) | 261 | 1000–1050 | 300 | Luigi Morrone | Italy
2020-04-20T11:31.7 | 108.5 | 0.50 (D) | 187 | 807–1050 | 247 | Tiziano Olivetti | Thailand
2020-04-25T17:27.0 | 113.4 | 0.35 (S) | 211 | 884–900 | 300 | Emmanuel Kardasis | Greece
2020-03-08T11:37.6 | 78.9 | 0.50 (D) | – | 320–380 | 1029 | Tiziano Olivetti | Thailand
2020-03-09T11:39.7 | 79.4 | 0.50 (D) | – | 320–380 | 804 | Tiziano Olivetti | Thailand
2020-03-10T17:04.0 | 79.9 | 0.35 (S) | – | 320–380 | 300 | Emmanuel Kardasis | Greece
2020-03-11T16:55.0 | 80.5 | 0.35 (S) | – | 320–380 | 300 | Emmanuel Kardasis | Greece
2020-03-13T11:31.3 | 81.6 | 0.50 (D) | – | 320–380 | 569 | Tiziano Olivetti | Thailand
2020-03-14T16:32.2 | 82.1 | 0.25 (M) | – | 320–380 | 500 | Christian Viladrich | France
2020-03-15T15:10.0 | 82.7 | 0.40 (D) | – | 350–410 | 180 | Yaroslav Naryzhniy | Ukraine
2020-03-16T15:00.0 | 83.3 | 0.40 (D) | – | 350–410 | 180 | Yaroslav Naryzhniy | Ukraine

Notes:
(N) Newtonian Telescope,
(S) Schmidt Cassegrain Telescope,
(D) Dall-Kirkham Telescope,
(K) Klevtzov Cassegrain Telescope.

The images from amateur observations were obtained thanks to the technique of
lucky imaging, where thousands of images are captured in a few minutes with
fast planetary cameras and then stacked to obtain a single one more nitid and
with better resolution (by up to a factor of four in poor seeing conditions)
Farsiu et al. (2004); Law et al. (2006). Software such as Autostakkert
(http://www.autostakkert.com , accessed on 2021) or Registax
(http://astronomie.be/registax, accessed on 2021) (Berrevoets et al., 2012)
were used following this strategy: (a) selection of the highest quality
images, (b) alignment, (c) stacking, and (d) processing of the final image to
get a sharper version Kardasis et al. (2016). To obtain these observations,
sensitive CMOS cameras (e.g., ZWO290 mm) were equipped on $\sim$0.2–0.5 m
aperture telescopes. In the case of NIR images, Venus’s clouds exhibit lower
contrast than in UV images, and a different image processing is required
(Peralta et al., 2019).

The first step in the analysis was to properly define the navigation of each
image, i.e., the correct orientation (north–south), fitting of the limb, and
the phase of the planet to a synthetic outline (calculated by software’s
ephemeris). For this, and further processing of the images, we used the free
software WinJupos (http://www.grischa-hahn.homepage.t-online.de , accessed on
2021) (Hahn and Jacquesson, 2012). We also used this software to measure the
motions of clouds patterns (latitude–longitude position, size, and wind
speed). The phase and the ephemeris of the planet are automatically calculated
by the software. Two rotation periods are implemented in WinJupos for the
determination of the planet longitudes: System 1 (S1) accounting for the
surface rotation with a period of $\sim$243 days, and System 2 (S2) with a
period of 4.2 days, corresponding to the superrotating upper clouds in UV
images (Sánchez-Lavega et al., 2017). For creating maps of the dayside middle
clouds with NIR images, we used a slower rotation rate with a period of almost
5 days (calculated from the average equatorial speed of CD events, see Table
4). CD speeds were measured using feature tracking of the following wave-front
edge (FWFE henceforth) in pairs of images (with time separation ranging from a
few hours to some days).

Most of the cloud motions were measured near the equator, where the CD is
customarily found and also because the spatial resolution approaches the best
that can be attained with ground-based telescopes (see Table 1). Using the
standard deviations of the measurements as a metric, we compared those derived
from the optimal fitting of the images with those obtained by slightly
changing the parameters that affect the navigation (orientation and limb
fitting), by adjusting to the largest and lowest acceptable values (by visual
inspection). We found that limb fitting (with an error of 2–3 pixels) was the
one resulting in systematically larger standard deviations than the
orientation. However, even in those cases, the standard deviations were always
smaller than those produced by the image resolution.

We obtained an estimate of the speed error by two methods. In the first
approach, we used the theoretical resolution of each telescope to calculate
the possible error ($\delta\upsilon=\delta x/\Delta t$, where $\delta x$
corresponds to the theoretical expected image resolution according to the
telescope’s aperture expressed in meters, and $\Delta t$ is the time
separation between image pairs). Image pairs separated by a temporal interval
of $\sim$2–8 h allow us to measure wind speeds with errors in the range
$\sim$10–40$\,\textrm{m}\cdot\textrm{s}^{-1}\,$. Image pairs separated by
$\Delta t\sim$ 5–10 days provide smaller errors in the range of $\sim$1
$\,\textrm{m}\cdot\textrm{s}^{-1}\,$, because of the large time difference. We
note here that these error estimates are lower limits, as they consider
optimal seeing conditions. However, given the lucky imaging approach and
deconvolution techniques, this can be lowered ever more. A second approach was
to measure wind speed from various latitude bins of 5–12.5∘. The standard
deviation of these measurements within each bin results in errors of $\sim$1
$\,\textrm{m}\cdot\textrm{s}^{-1}\,$for $\Delta t\sim$ 5–10 days and
$\sim$5–25 $\,\textrm{m}\cdot\textrm{s}^{-1}\,$for $\Delta t\sim$ 2–8 h,
respectively. Given the similarity between the two approaches, we opted to
quote the error estimates as derived from the second method (i.e., the
statistical approach). Although this approach depends on the determination of
the CD by “visual inspection”, it takes into account the high quality of the
selected images.

Table 4: Summary of cloud discontinuity events’ properties as derived from
this work. The first two columns present the pair of images that were used to
extract the wind measurements near the equator. For each detection, we show
the orientation (column 3), the maximum latitude range (column 4), the maximum
length (column 5), the width at the equator (column 6), and the derived
rotation period (column 7) and zonal speed (column 8).

Date and Time (UT) 1 | Date and Time (UT) 1 | Orientation 2 | Maximum | Maximum | Width at | Rotation | Zonal
---|---|---|---|---|---|---|---
(1st obs.) | (2nd obs.) |  | Latitude Range | Length3 | Equator 4 | Period | Speed 5
(YYYY-MM-DDThh:mm) | (YYYY-MM-DDThh:mm) | (degrees) |  (degrees) | ($10^{3}\,\textrm{km}$) | ($10^{3}\,\textrm{km}$) | (Earth days) | (m$\cdot$ s-1)
2019-11-23T08:14* | 2019-11-23T15:51 | 124 | 20∘ N–9∘ S | 3.3 | 0.6 | 4.1 $\pm$ 0.9 | 105 $\pm$ 19
2019-11-23T15:51 | 2019-12-03T15:20* | 97 | 26∘ N–28∘ S | 5.7 | 0.6 | 5.06 $\pm$ 0.01 | 84.9 $\pm$ 0.1
2019-12-03T15:20 | 2019-12-08T15:29* | 39 | 14∘ N– 17∘ S | 5.1 | 1.1 | 5.06 $\pm$ 0.01 | 84.9 $\pm$ 0.1
2019-12-18T08:09* | 2019-12-18T15:23 | 77 | 27∘ N–28∘ S | 5.9 | 0.6 | 4.2 $\pm$ 0.2 | 104 $\pm$ 5
2020-01-11T11:34* | 2020-01-11T14:51 | 137 | 12∘ N–14∘ S | 3.7 | 0.5 | 4.3 $\pm$ 0.8 | 100 $\pm$ 16
2020-02-01T11:27 | 2020-02-01T15:51* | 39 | 19∘ N–25∘ S | 6.5 | 0.6 | 4.0 $\pm$ 0.4 | 108 $\pm$ 9
2020-02-16T07:55* | – | 84 | 13∘ N–25∘ S | 4.0 | 0.3 | – | –
2020-02-24T17:27 | 2020-02-24T18:36* | 85 | 28∘ N–28∘ S | 6.1 | 0.5 | 3.5 $\pm$ 0.9 | 125 $\pm$ 26
2020-02-24T18:36 | 2020-02-28T17:41* | 66 | 15∘ N–13∘ S | 3.3 | 0.6 | 4.18 $\pm$ 0.02 | 103.4 $\pm$ 0.5
2020-02-28T17:41 | 2020-03-11T16:45* | 100 | 25∘ N–30∘ S | 5.9 | 0.5 | 6.00 $\pm$ 0.01 | 71.6 $\pm$ 0.1
2020-03-11T16:45 | 2020-03-16T15:00* | 95 | 23∘ N–36∘ S | 6.3 | 0.7 | 5.14 $\pm$ 0.02 | 83.9 $\pm$ 0.3
2020-03-16T15:00 | 2020-03-21T17:35* | 88 | 27∘ N–32∘ S | 6.3 | 0.5 | 5.10 $\pm$ 0.01 | 84.5 $\pm$ 0.2
2020-03-21T11:34 | 2020-03-21T17:35 | – | – | – | – | 5.09 $\pm$ 0.08 | 84.7 $\pm$ 1.3
2020-03-21T17:35 | 2020-03-26T15:20* | 94 | 13∘ N–33∘ S | 5.0 | 0.4 | 5.10 $\pm$ 0.02 | 84.6 $\pm$ 0.4
2020-03-26T15:20 | 2020-03-31T17:04* | 99 | 20∘ N–29∘ S | 5.2 | 1.0 | 5.15 $\pm$ 0.03 | 83.6 $\pm$ 0.5
2020-03-31T17:04 | 2020-04-05T16:40* | 90 | 31∘ N–25∘ S | 6.0 | 0.5 | 4.91 $\pm$ 0.03 | 87.8 $\pm$ 0.5
2020-04-05T16:40 | 2020-04-10T16:12* | 93 | 31∘ N–29∘ S | 6.5 | 0.6 | 5.09 $\pm$ 0.01 | 84.7 $\pm$ 0.2
2020-04-10T16:12 | 2020-04-20T11:32* | 80 | 10∘ N–20∘ S | 3.3 | 0.3 | 5.05 $\pm$ 0.01 | 85.4 $\pm$ 0.1
AVERAGE | 88 | 20.8∘ N–21.5∘ S | 5.2 | 0.6 | 4.8 $\pm$ 0.2 | 92.1 $\pm$ 4.6

Notes:
1 We note with an * the observations from which the properties of CD event
were measured.
2 Degrees relative to equator (scale 0–180∘) right to left.
3 North–south direction, error of $\pm 0.2\times 10^{3}\,\textrm{km}$.
4 East–west direction, error of $\pm 0.2\times 10^{3}\,km$.
5 Near equator, i.e., within 5∘ N–5∘ S.

The position and the inclination of the CD were also measured using the FWFE.
Inclinations were measured with respect to the equator (0–180∘). Maps were
created with equirectangular projection that simply displays the relative
position of large formations. The gaps correspond to areas where no
observations were performed or where the spatial resolution was very poor.
Latitudes higher then 60∘ were also excluded because of the low resolution and
severe distortion near the polar region expected in equirectangular
projections.

###  2.2 Venus Images Taken with JAXA’s Akatsuki

Amateur observations during March 2020 were compared with nearly simultaneous
observations from Akatsuki (Murakami et al., 2019; Nakamura et al., 2016).
This gave us the chance to test whether the CD events can be visualized at
other vertical levels and if there is any interaction with other planetary
scale waves at the upper clouds. Such phenomena include the Y-feature, a dark
cloud structure observed for decades in UV images that looks like a
horizontally oriented Y letter visible every $\sim$4 days (Peralta et al.,
2015; Imai et al., 2019), and the stationary bow-shape gravity waves (Kouyama
et al., 2017; Peralta et al., 2017). These phenomena at the upper clouds can
be observed in images taken by the Longwave Infrared Camera (LIR) (Fukuhara et
al., 2011) and the Ultraviolet Imager (UVI) (Yamazaki et al., 2018)) onboard
Akatsuki . The bolometer LIR maps the brightness temperature of the upper
clouds within $\sim$60–75 km (Taguchi et al., 2007) in the spectral band 8–12
$\mu$m, both on the dayside and nightside, while the camera UVI senses the
albedo of the cloud tops at $\sim$70 km (Yamazaki et al., 2018; Horinouchi et
al., 2018) using filters centered at 283 and 365 nm. Five events of the CD
during March (11 March 2020, 16 March 2020, 21 March 2020, 26 March 2020, and
31 March 2020) were studied in combination with Akatsuki images, which covered
orbits 142–144 of the mission (Murakami et al., 2017, 2018, 2017, 2018).
During March 2020, the phase angle of Venus as observed from Akatsuki varied
between 30∘ and 160∘, while the spatial resolution of UVI and LIR images at
lower latitudes ranged 12–90 and 40–360 km/pixel, respectively. The UVI images
were photometrically corrected, i.e., correcting for the limb-darkening effect
(Lee et al., 2017). Both UVI and LIR images were processed and projected as
described by Gonçalves et al. (2020). A devoted software written in IDL was
used to process and perform measurements in the Akatsuki images (Peralta et
al., 2018).

##  3 Results and Discussion

Three years after its last detection on the dayside middle clouds in December
2016 (Peralta et al., 2020), the CD was captured again in November–December of
2019 (see Figures 1 and 2). Unfortunately, from mid-December 2019 until early
March 2020, we have gaps in the observations that prevented a continuous
monitoring of these new events of the CD. From March 2020, the presence of a
CD became more evident, manifesting as a long dark feature sometimes followed
by a brighter streak (see Figure 1). The larger number of observations enabled
us to monitor with more detail this phenomenon until the middle of April, when
the CD seemed to vanish.

Figure 1: Observations of cloud discontinuities, observed in the 2019/2020
eastern elongation of Venus, showing different morphologies. Details of the
observations are presented in Table 1. Figure 2: Time composites of the of the
middle clouds (NIR $\sim$900–1000 nm) of Venus showing large atmospheric
features when the cloud discontinuity is observed. These large-scale features
are marked with red ovals. Equirectangular projections (60∘ N to 60∘ S) placed
from left to right with increasing dates: (A) 14–18 November 2019 (A. Wesley,
E. Kardasis), (B) 19–23 November 2019 (A. Wesley, E. Kardasis), (C) 24–28
November 2019 (A. Wesley, E. Kardasis), (D) 29 November–3 December 2019 (A.
Wesley, E. Kardasis, N. MacNeill), (E) 4–8 December 2019 (A. Wesley,
E.Kardasis). Details of the observations are presented in Table 1.

###  3.1 CD during November–December 2019

The first clear manifestation of CD occurred in late November 2019 (Figures 1
and 2). On 18 November 2019, Venus displayed no features on the area where we
would expect to observe the CD. There are three dark stripes or bands
(Northern, Equatorial, Southern) in the proceeding area. On 23 November 2019,
the CD seems apparent as an extension of the Northern Dark Stripe, which was
previously reported as a cloud pattern frequently preceding the arrival of the
discontinuity on the dayside middle clouds during the year 2016 (Peralta et
al., 2020) (see Figure 1b). On 3 December 2019, the CD is observed between 26∘
N and 28∘ S, bending at its northern end and eastward extending $\sim$5700 km
as an horizontal sharp dark stripe similar to those reported by Peralta et al.
(2019) on the nightside lower clouds. On the west of CD, an equatorial stripe
was formed also. In the central image of Figure 2D, a bright equatorial area
is observed, surrounded north and south with wavy borders.

Taking advantage of the smaller phase angle for these earlier observations, we
studied the morphology and development of CD along several rotations through
maps of time composites (see Figure 2) similar to the ones created by Peralta
et al. (2020) for the nightside lower clouds. Due to the difficulty in
accurately measuring the middle clouds’ mean flow with amateur images, the
equirectangular projections of the NIR images were combined after shifting
them according to the average CD speed ($\sim$5 day period, see Peralta et al.
(2020)) at the equator rather the mean wind speed from passive tracers. During
the last rotation (Figure 2E), the CD changed its inclination (starting with a
tilt of 124∘ relative to the parallels, then 97∘ and finally 39∘) and a bright
area follows the CD in the southward of the equator. The CD drifted to the
west with an average equatorial zonal speed of $-94.7\pm
6$$\,\textrm{m}\cdot\textrm{s}^{-1}\,$, with a derived period of $4.6\pm 0.3$
days.

###  3.2 CD during January–April 2020

From November–December 2019 to March 2020, the increasing phase angle and more
sparse data coverage hindered the monitoring of the CD, though some of the
images show its presence (see Figure 1). During this period, the CD suffered
noticeable morphological changes after each full revolution, which is
consistent with the observations by Akatsuki /IR1 during the year 2016
(Peralta et al., 2020). Between January 2020 and April 2020, the CD was always
observed between 31∘ N and 36∘ S, with a length within 3300–6500 km and tilt
relative to the parallels ranging 39–137∘. Its average speed was $-91\pm
4$$\,\textrm{m}\cdot\textrm{s}^{-1}\,$with a period of $4.8\pm 0.2$ days.

On 11 March 2020, our observations exhibited an intensification of the CD,
being apparent as a long dark sharp discontinuity followed by a brighter
cloudy area, and barely suffering distortions on its shape and size along 11,
16, and 21 March 2020. A similar event was also reported in August 2016
(Peralta et al., 2020). Combining images from several ground-based observers,
we were able to track the CD until 25 April 2020, and we estimated for the CD
a rotation period of about 5 days. From 26 March 2020 onwards, it appears to
have shrunk, partly dissipated, and displayed much lower contrast during the
next weeks. In general, it spanned between 31∘ N and 36∘ S, ranging 3300–6500
km in length. The dark streak was sometimes followed by a bright one. The CD
average width at the equator was $\sim$600 km, with a lower value of $\sim$300
km (similar to the spatial resolution of the images) and up to $\sim$ 1100 km.
The March–April event presented some constant properties such as inclination
(almost vertical relative to the equator, 80–100∘) and speed ($-$83.6 to
$-$87.8$\,\textrm{m}\cdot\textrm{s}^{-1}\,$, with an average of $84.9\pm 0.4$
$\,\textrm{m}\cdot\textrm{s}^{-1}\,$, and a period $5.08\pm 0.03$ days.).

Figure 3 exhibits the zonal speeds of the discontinuity during the CD events
of 11, 16, 21, 26, and 31 March 2020 apparent in images of Venus taken by
ground-based observers (see Figure 1). These zonal speeds are compared with
the zonal winds measured with 365 nm images taken by Akatsuki /UVI during the
same dates as the CD events occurring during 11, 16, 21, and 26 March 2020. In
the case of the winds at the middle clouds, it was not possible to identify
reliable passive tracers in the NIR ground-based images due to their low
contrast, and the profile of winds obtained by Akatsuki /IR1 in 2016 (Peralta
et al., 2019) is shown instead. The CD speeds are consistent with the
manifestation of a Kelvin wave (Peralta et al., 2020), with phase speeds
faster than the background winds at the middle clouds observed in the 900 nm
images from Akatsuki /IR1. On the other hand, the phase speeds of the CD are
slower than the winds measured in the UVI 365 nm images sensing the top of the
upper clouds, supporting the argument that the CD cannot be observed at the
clouds’ top because it becomes dissipated during its vertical propagation
(Peralta et al., 2020).

Figure 3: Zonal speeds of the clouds’ discontinuity (CD) during March 2020
(dots) compared with the zonally-averaged profiles of the zonal wind for the
top of the upper clouds at $\sim$70 km (violet line) and the middle clouds
within 50–55 km (Peralta et al., 2019) (red line). The speeds of the
discontinuity and the winds at the cloud tops were obtained from ground-based
near-infrared images and Akatsuki /UVI 365 nm images respectively, considering
observations matching the CD events: 11, 16, 21, 26, and 31 March 2020. For
the winds at the middle clouds, we considered the profile obtained by Peralta
et al. (2019) during 2016. As a result of the worse spatial resolution at
higher latitudes, speed vectors had larger error bars (errors of about $\pm
15$$\,\textrm{m}\cdot\textrm{s}^{-1}\,$for $\Delta$t$\sim$6 hours and $\pm
5$$\,\textrm{m}\cdot\textrm{s}^{-1}\,$for $\Delta$t$\sim$5 terrestrial days).

###  3.3 Long-Term Evolution

Figure 4 exhibits the rotation period, orientation, and latitude coverage of
the discontinuity during the years 2016–2020 from Venus images taken by the
orbiter Akatsuki (Peralta et al., 2020) and ground-based observations by
NASA’s Infrared Telescope Facility (Peralta et al., 2020) and small telescopes
(this work). Concerning the rotation period of the CD, we obtain a period of
$4.8\pm 0.2$ days considering the averaged zonal speed ($92.1\pm
4.6$$\,\textrm{m}\cdot\textrm{s}^{-1}\,$) from the 17 CD events documented in
this work (see Table 4), while the mean rotation period considering the full
data from 2016 to 2020 can be updated to $4.9\pm 0.5$ days, the same as
estimated by Peralta et al. (2020). The discontinuity exhibited noticeable
variations in its orientation between November 2019 and February 2020, with a
clear tendency to keep an orientation perpendicular to the equator from March
until April 2020. Given that gravity waves are subject to refraction as they
propagate, due to the varying background atmospheric temperature structure and
winds (Heale and Snively, 2015), these changes in the orientation of the
discontinuity might be tentatively attributed to a slight atmospheric
refraction whose evidence should be investigated in future works. Finally, the
CD events reported in this work confirm the hemispherical asymmetry for the
discontinuity previously reported by Peralta et al. (2020).

Figure 4: Long-term evolution of CD properties. Combination of past
(2016–2018) with current (2019–2020) measurements (marked with the purple
outline). The rotation period, orientation, and latitude coverage of the
disruption along 2016–2020 are exhibited. The periods were measured from the
position of the disruption at the equator in images separated by hours and
several days. When the CD did not intersect the equator, we considered its
longitude closest to the equator. The mean period is shown with a blue dotted
line. Day and night-side data are shown in green and red, respectively.

###  3.4 The Cloud Tops during the CD Events of March 2020

The Y-feature is the most noticeable planetary-scale pattern that can be
identified in the ultraviolet images that sense the top of the upper clouds at
the dayside of Venus (Boyer and Camichel, 1961; Belton et al., 1976; Rossow et
al., 1980; Kouyama et al., 2012). The Y-feature shares some of the
characteristics of the lower clouds’ discontinuity: the Y-feature is a
recurrent phenomenon subject to cycles of generation and destruction (Rossow
et al., 1980), it manifests about the equator, it propagates faster than the
mean zonal flow (Kouyama et al., 2012), and it has been long interpreted as
the manifestation of a kelvin-type wave (Peralta et al., 2015). On the other
hand, the Y-feature exhibits tends to have a rotation period of about 4 days
(Imai et al., 2019) (generally faster than the discontinuity, see Figure 4),
it can extend to higher latitudes (Rossow et al., 1980; Peralta et al., 2015),
and it has never been reported on observations of altitudes below the top of
the clouds (Belton et al., 1991; Peralta et al., 2019), which is in agreement
with the idea that the Y-feature may be vertically trapped (Peralta et al.,
2015).

Figures 5 and 6 display a comparison between the morphology of the cloud tops
(UV images) and the middle clouds (NIR images) during some of the CD events
occurring in March 2020. The cycle of the discontinuity at the dayside middle
clouds is confronted by that of the Y-feature at the cloud tops in Figure 5,
using NIR and UV images obtained by ground-based observers and time composites
similar to those published by Sánchez-Lavega et al. (2016). When comparing the
cycles of 7–11 March and 12–16 March 2020, we observe that the position of the
sharp vertical streak-form of the CD (evident in NIR images) does not seem
linked to the propagation of the Y-feature in UV images or to any specific
cloud pattern in its cycle. The time composites also display an evident phase-
lag between both phenomena, supporting that we may be observing independent
waves that propagate at different vertical levels.

Figure 5: Comparison among simultaneous NIR amateur maps showing the CD
position in the days observed in relation to Y-feature. NIR map between 7–11
March 2020 (A) and 12–16 March 2020 (C) from amateur observations compared
with UV amateur maps presenting the Y-feature on 8 March 2020 to 11 March 2020
(B) and 13 March 2020 to 16 March 2020 (D). NIR and UV maps were constructed
by taking as reference images the simultaneous observations in both
wavelengths when the CD was observed.

Akatsuki /UVI imaging of the top of the clouds provide higher resolution to
study the details of individual cloud features and patterns that cannot be
resolved by ground-based telescopes. In fact, UVI allows observing Venus with
two filters: 283 nm to sense the absorption of the SO2 and 365 nm to observe
the absorption due to the unknown absorber (Yamazaki et al., 2018). Figure 6
displays a representative sampling of upper clouds in images obtained very
close to the time of the NIR amateur data and covers the exact area of CD
observations on March 2020. The sharp vertical form of CD (Figure 6A) due to
variations of the middle-cloud optical thickness is also absent in high-
quality UV data (Figures 6B,C), confirming previous findings by Peralta et al.
(2020) and supporting the idea that the critical level where the CD dissipates
should be located between the middle clouds and upper clouds’ top (see Figure
3).

Figure 6: Contemporaneous views of NIR amateur set ($\sim$700-1000 nm) where
the cloud discontinuity is observed (A), along with Akatsuki /UVI at 365 nm
(B), UVI at 285 nm (C), and LIR at 10 $\mu$m (D) where the cloud discontinuity
seems absent.

###  3.5 Thermal Emission from the Upper Clouds during the CD Events of March
2020

Large stationary gravity waves appear mostly above four specific highland
regions in the low latitudes, specially when these regions are in the local
afternoon (Kouyama et al., 2017). The exact mechanism for generating these
stationary waves at the cloud top is as yet unclear, although there is a
general agreement that we must be visualizing Lee waves triggered at the
surface of Venus (stationary waves in the atmosphere, happening where winds
encounter obstacles such as high mountains; Peralta et al. (2017); Fukuhara et
al. (2017); Navarro et al. (2018); Kitahara et al. (2019); Lefèvre et al.
(2020)). These planetary scale waves tend to be aligned almost north–south
with a slight curvature, and they present silhouette similarities with the
cloud discontinuity.

During March 2020, Akatsuki /LIR exhibits clear stationary waves in the
brightness temperature images taken during 11–16 and 21 of March 2020. The
wavefronts were located above Alta Regio and almost aligned in a north–south
direction spanning several thousand kilometers. They were not visible during
26 and 31 March 2020. In Figure 6D, the discontinuity is not apparent in LIR
images acquired in almost the same time as the CD events, while the LIR image
taken on 26 March 2020 does not cover the region where the CD was observed at
the middle clouds with ground-based observations (Figure 6A). Since the
weighting function of LIR extends within about 50–60 km above the surface
(Taguchi et al., 2007), LIR might be (partially) sensitive to the thermal
emission from middle clouds (50.5–56.5 km) as well. In this sense, the absence
of the CD in the LIR images may provide new hints to better constrain the
altitude range for the critical level of the Kelvin wave interpreted to be
responsible for the CD.

##  4 Summary

In this work, we report the reappearance of the equatorial discontinuity at
the dayside middle clouds of Venus (Peralta et al., 2020), which is a global-
scale atmospheric wave comparable with other planetary phenomena such as the
Y-feature (Peralta et al., 2015) and the giant stationary bow-shape wave
(Fukuhara et al., 2017). We studied its morphology, the evolution, and the
properties of the most prominent events in 2019/2020 eastern elongation of
Venus, combining both NIR and UV data from ground-based amateur observations
and Akatsuki imaging.

The most intense event occurred in March 2020, which was visible in NIR images
as a dark vertical streak followed by a bright streak. In general, the CD
events extended approximately 20∘ above and below the equator, displaying a
variable length (3300–6500 km) with widths in the range of $\sim$300–1100 km
(with an average of 600 km). Although the orientation displayed noticeable
variations, it is generally perpendicular to the equatorial plane. In total,
we have measured the mean rotational period from all CD events to be $\sim$4.9
days with a mean wind speed of $\sim$92.1
$\,\textrm{m}\cdot\textrm{s}^{-1}\,$, which is consistent with Peralta et al.
(2020). We also confirm that the morphology of the CD is often hemispherically
asymmetric (Peralta et al., 2020).

In general, the speeds we obtained are higher/lower when compared with the
average zonal speeds of the middle/upper clouds found by Akatsuki in
2016/2020. In combination with imaging, we conclude that the CD is a
phenomenon limited to the middle clouds of the atmosphere of Venus. When
combining ground-based amateur images with Akatsuki UV and LIR imaging, we
find no direct relation between the CD in the middle and the phenomena in
upper clouds (such as Y-feature or the bow-shape waves). Moreover, the evident
phase lags between these phenomena point to a scenario with independent waves,
which propagate at different vertical levels. Our results further support that
the CD dissipates during its vertical propagation, with the critical level
located between the middle and upper clouds.

With the current work, we aim to both contribute to the additional
characterization of this phenomenon and to highlight the importance of the
amateur contributions. Amateur observations can play an important role
especially when used complementary to professional ones. Moreover, systematic
monitoring by amateurs can trigger, and of course, support coordinated
observing campaigns, such as the next campaigns of observations coordinated
with Akatsuki (http://pvol2.ehu.eus/bc/Venus/) of Venus flybys of NASA’s
PARKER (Wood et al., 2022) on 21 August 2023 and 6 November 2024. The
professional–amateur collaboration offers an invaluable resource in the
advance of our understanding of the atmosphere of Venus.

## Funding

J.P. thanks EMERGIA funding from Junta de Andalucía Spain (code:
EMERGIA20_00414). G.M. acknowledges funding support from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 772086). M.I. is supported by the Japan Society
for the Promotion of Science (JSPS) (Grant-in-Aid for JSPS Research Fellow:
JP21J00351).

## Acknowledgments

We are grateful for support from all the members of the Akatsuki mission. We
are also grateful to many amateur observers that observed Venus intensively
during the investigated period providing data for this research. Special
thanks to the observer Raimondo Sedrani for using some of his observations.
E.K. would like to offer his special thanks to Grischa Hahn developer of
WinJupos software for providing help and advice for the proper use. Moreover,
E.K. would like to acknowledge the continuous support of his wife Dimitra.

## Data availability statement

This study was based on ground-based images of Venus provided by worldwide
amateur observers and Akatsuki observational data provided by ISAS. Ground-
based observations are publicly available in the Venus section of ALPO-Japan
at http://alpo-j.sakura.ne.jp/Latest/Venus.htm (accessed on 2021). All
original Akatsuki/UVI and LIR data are available via Data ARchives and
Transmission System (DARTS), provided by the Center for Science-satellite
Operation and Data Archive (C-SODA) at ISAS/JAXA, and NASA PDS Atmospheres
Node. The analyzed data are the L2b and L3bx products, and these data are
provided in the UVI datasets (accessed on 2020)
https://data.darts.isas.jaxa.jp/pub/pds3/vco-v-uvi-3-cdr-v1.0/ and at
https://data.darts.isas.jaxa.jp/pub/pds3/extras/vco_uvi_l3_v1.0/, and LIR
datasets at https://data.darts.isas.jaxa.jp/pub/pds3/vco-v-lir-3-cdr-v1.0/ and
https://data.darts.isas.jaxa.jp/pub/pds3/extras/vco_lir_l3_v1.0/.

## References

  * Belton et al. (1976) Belton, M.J.S.; Smith, G.R.; Schubert, G.; del Genio, A.D.  Cloud patterns, waves and convection in the Venus atmosphere.  J. Atmos. Sci. 1976, 33, 1394–1417.
  * Belton et al. (1991) Belton, M.J.S.; Gierasch, P.J.; Smith, M.D.; Helfenstein, P.; Schinder, P.J.; Pollack, J.B.; Rages, K.A.; Morrison, D.; Klaasen, K.P.; Pilcher, C.B.  Images from Galileo of the Venus cloud deck.  Science 1991, 253, 1531–1536.
  * Berrevoets et al. (2012) Berrevoets, C.; DeClerq, B.; George, T.; Makolkin, D.; Maxson, P.; Pilz, B.; Presnyakov, P.; Roel, E.; Weiller, S.  RegiStax: Alignment, stacking and processing of images. _arXiv_ 2012, arXiv:1206.001.
  * Boyer and Camichel (1961) Boyer, C.; Camichel, H.  Observations photographiques de la planète Vénus.  Ann. D’Astrophysique 1961, 24, 531.
  * Farsiu et al. (2004) Farsiu, S.; Robinson, M.D.; Elad, M.; Milanfar, P.  Fast and Robust Multiframe Super Resolution.  IEEE Trans. Image Process. 2004, 13, 1327–1344, doi:10.1109/TIP.2004.834669.
  * Fukuhara et al. (2011) Fukuhara, T.; Taguchi, M.; Imamura, T.; Nakamura, M.; Ueno, M.; Suzuki, M.; Iwagami, N.; Sato, M.; Mitsuyama, K.; Hashimoto, G.L.; et al.  LIR: Longwave Infrared Camera onboard the Venus orbiter Akatsuki.  Earth Planets Space 2011, 63, 1009–1018, doi:10.5047/eps.2011.06.019.
  * Fukuhara et al. (2017) Fukuhara, T.; Futaguchi, M.; Hashimoto, G.L.; Horinouchi, T.; Imamura, T.; Iwagaimi, N.; Kouyama, T.; Murakami, S.Y.; Nakamura, M.; Ogohara, K.; et al.  Large stationary gravity wave in the atmosphere of Venus.  Nat. Geosci. 2017, 10, 85–88, doi:10.1038/ngeo2873.
  * Gonçalves et al. (2020) Gonçalves, R.; Machado, P.; Widemann, T.; Peralta, J.; Watanabe, S.; Yamazaki, A.; Satoh, T.; Takagi, M.; Ogohara, K.; Lee, Y.J.; et al.  Venus’ cloud top wind study: Coordinated Akatsuki/UVI with cloud tracking and TNG/HARPS-N with Doppler velocimetry observations.  Icarus 2020, 335, 113418, doi:10.1016/j.icarus.2019.113418.
  * Hahn and Jacquesson (2012) Hahn, G.; Jacquesson, M.  WinJUPOS-Database for Object Positions on Planets and the Sun.  Available online: http://jupos.privat.t-online.de (accessed in 2021).
  * Heale and Snively (2015) Heale, C.J.; Snively, J.B.  Gravity wave propagation through a vertically and horizontally inhomogeneous background wind.  J. Geophys. Res.  2015, 120, 5931–5950, doi:10.1002/2015JD023505.
  * Horinouchi et al. (2018) Horinouchi, T.; Kouyama, T.; Lee, Y.J.; Murakami, S.y.; Ogohara, K.; Takagi, M.; Imamura, T.; Nakajima, K.; Peralta, J.; Yamazaki, A.; Yamada, M.; Watanabe, S.  Mean winds at the cloud top of Venus obtained from two-wavelength UV imaging by Akatsuki.  Earth Planets Space 2018, 70, 10, doi:10.1186/s40623-017-0775-3.
  * Horinouchi et al. (2020) Horinouchi, T.; Hayashi, Y.Y.; Watanabe, S.; Yamada, M.; Yamazaki, A.; Kouyama, T.; Taguchi, M.; Fukuhara, T.; Takagi, M.; Ogohara, K.; et al.  How waves and turbulence maintain the super-rotation of Venus’ atmosphere.  Science 2020, 368, 405–409, doi:10.1126/science.aaz4439.
  * Imai et al. (2019) Imai, M.; Kouyama, T.; Takahashi, Y.; Yamazaki, A.; Watanabe, S.; Yamada, M.; Imamura, T.; Satoh, T.; Nakamura, M.; Murakami, S.y.; Ogohara, K.; Horinouchi, T.  Planetary-Scale Variations in Winds and UV Brightness at the Venusian Cloud Top: Periodicity and Temporal Evolution.  J. Geophys. Res. Planets 2019, 124, 2635–2659, doi:10.1029/2019JE006065. Available online: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JE006065 (accessed in 2021).
  * Iwagami et al. (2018) Iwagami, N.; Sakanoi, T.; Hashimoto, G.L.; Sawai, K.; Ohtsuki, S.; Takagi, S.; Uemizu, K.; Ueno, M.; Kameda, S.; Murakami, S.Y.; et al. Initial products of Akatsuki 1-$\mu$m camera.  Earth Planets Space 2018, 70, 6, doi:10.1186/s40623-017-0773-5.
  * Kardasis et al. (2016) Kardasis, E.; Rogers, J.H.; Orton, G.; Delcroix, M.; Christou, A.; Foulkes, M.; Yanamandra-Fisher, P.; Jacquesson, M.; Maravelias, G.  The need for professional-amateur collaboration in studies of Jupiter and Saturn.  J. Br. Astron. Assoc. 2016, 126, 29–39.
  * Kitahara et al. (2019) Kitahara, T.; Imamura, T.; Sato, T.M.; Yamazaki, A.; Lee, Y.J.; Yamada, M.; Watanabe, S.; Taguchi, M.; Fukuhara, T.; Kouyama, T.; Murakami, S.y.; Hashimoto, G.L.; Ogohara, K.; Kashimura, H.; Horinouchi, T.; Takagi, M.  Stationary Features at the Cloud Top of Venus Observed by Ultraviolet Imager Onboard Akatsuki.  J. Geophys. Res. (Planets) 2019, 124, 1266–1281, doi:10.1029/2018JE005842.
  * Kouyama et al. (2012) Kouyama, T.; Imamura, T.; Nakamura, M.; Satoh, T.; Futaana, Y.  Horizontal structure of planetary-scale waves at the cloud top of Venus deduced from Galileo SSI images with an improved cloud-tracking technique.  Planet. Space Sci. 2012, 60, 207–216, doi:10.1016/j.pss.2011.08.008.
  * Kouyama et al. (2017) Kouyama, T.; Imamura, T.; Taguchi, M.; Fukuhara, T.; Sato, T.M.; Yamazaki, A.; Futaguchi, M.; Murakami, S.; Hashimoto, G.L.; Ueno, M.; et al.  Topographical and Local Time Dependence of Large Stationary Gravity Waves Observed at the Cloud Top of Venus.  Geophys. Res. Lett. 2017, 44, 12,098–12,105.  doi: 10.1002/2017GL075792.
  * Law et al. (2006) Law, N.M.; Mackay, C.D.; Baldwin, J.E.  Lucky imaging: High angular resolution imaging in the visible from the ground.  Astron. Astrophys. 2006, 446, 739–745, doi:10.1051/0004-6361:20053695.
  * Lee et al. (2017) Lee, Y.J.; Yamazaki, A.; Imamura, T.; Yamada, M.; Watanabe, S.; Sato, T.M.; Ogohara, K.; Hashimoto, G.L.; Murakami, S.  Scattering Properties of the Venusian Clouds Observed by the UV Imager on board Akatsuki.  Astron. J. 2017, 154, 44, doi:10.3847/1538-3881/aa78a5.
  * Lefèvre et al. (2020) Lefèvre, M.; Spiga, A.; Lebonnois, S.  Mesoscale modeling of Venus’ bow-shape waves.  Icarus 2020, 335, 113376, doi:10.1016/j.icarus.2019.07.010.
  * Limaye et al. (2018) Limaye, S.S.; Watanabe, S.; Yamazaki, A.; Yamada, M.; Satoh, T.; Sato, T.M.; Nakamura, M.; Taguchi, M.; Fukuhara, T.; Imamura, T.; et al. Venus looks different at different wavelengths: Morphology from Akatsuki multispectral images.  Earth Planets Space 2018, 70, 38, doi:10.1186/s40623-018-0789-5.
  * McGouldrick et al. (2012) McGouldrick, K.; Momary, T.W.; Baines, K.H.; Grinspoon, D.H.  Quantification of middle and lower cloud variability and mesoscale dynamics from Venus Express/VIRTIS observations at 1.74 $\mu$m.  Icarus 2012, 217, 615–628, doi:10.1016/j.icarus.2011.07.009.
  * McGouldrick et al. (2021) McGouldrick, K.; Peralta, J.; Barstow, J.K.; Tsang, C.C.C.  Using VIRTIS on Venus Express to Constrain the Properties of the Giant Dark Cloud Observed in Images of Venus by IR2 on Akatsuki.  Planet. Sci. J. 2021, 2, 153, doi:10.3847/PSJ/ac0e39.
  * Murakami et al. (2017) Murakami, S.; Yamada, M.; Yamazaki, A.; McGouldrick, K.; Yamamoto, Y.  _Venus Climate Orbiter Akatsuki UVI Calibrated Data v1.0_ ; VCO-V-UVI-3-CDR-V1.0, NASA Planetary Data System; 2017; doi:10.17597/isas.darts/vco-00003..
  * Murakami et al. (2017) Murakami, S.; Yamada, M.; Yamazaki, A.; McGouldrick, K.; Yamamoto, Y.  _Venus Climate Orbiter Akatsuki LIR Calibrated Data v1.0_ ; VCO-V-LIR-3-CDR-V1.0, NASA Planetary Data System; 2017; doi:10.17597/isas.darts/vco-00012.
  * Murakami et al. (2018) Murakami, S.; Ogohara, K.; Takagi, M.; Kashimura, H.; Yamada, M.; Kouyama, T.; Horinouchi, T.; Imamura, T.  _Venus Climate Orbiter Akatsuki UVI Longitude-Latitude Map Data v1.0_ ; JAXA Data Archives and Transmission System; 2018; doi:10.17597/isas.darts/vco-00016
  * Murakami et al. (2018) Murakami, S.; Ogohara, K.; Takagi, M.; Kashimura, H.; Yamada, M.; Kouyama, T.; Horinouchi, T.; Imamura, T.  _Venus Climate Orbiter Akatsuki LIR Longitude-Latitude Map Data v1.0_ ; JAXA Data Archives and Transmission System; 2018; doi:10.17597/isas.darts/vco-00019.
  * Murakami et al. (2019) Murakami, S.; Yamamoto, Y. McGouldrick, K.; Hashimoto, G.L.; Yamada, M.; Yamazaki, A.; Sato, T.M.; Hirata, N.; Ohtsuki, S.; Kameda, S.; Satoh, T.; Kouyama, T.; Fukuhara, T.; Taguchi, M.; Ando, H.; Imamura, T.; Ogohara, K.; Takagi, M.; Kashimura, H.; Horinouchi, T.  Venus Climate Orbiter Akatsuki Dataset Collection. In JAXA Data Archives and Transmission System; Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, 2019, doi:10.17597/ISAS.DARTS/VCO-00000.
  * Nakamura et al. (2016) Nakamura, M.; Imamura, T.; Ishii, N.; Abe, T.; Kawakatsu, Y.; Hirose, C.; Satoh, T.; Suzuki, M.; Ueno, M.; Yamazaki, A.; et al. AKATSUKI returns to Venus.  Earth Planets Space 2016, 68, 1–10, doi:10.1186/s40623-016-0457-6.
  * Navarro et al. (2018) Navarro, T.; Schubert, G.; Lebonnois, S.  Atmospheric mountain wave generation on Venus and its influence on the solid planet’s rotation rate.  Nat. Geosci. 2018, 11, 487–491, doi:10.1038/s41561-018-0157-x.
  * Peralta et al. (2015) Peralta, J.; Sánchez-Lavega, A.; López-Valverde, M.A.; Luz, D.; Machado, P.  Venus’s major cloud feature as an equatorially trapped wave distorted by the wind.  Geophys. Res. Lett. 2015, 42, 705–711.  2014GL062280, doi:10.1002/2014GL062280.
  * Peralta et al. (2017) Peralta, J.; Hueso, R.; Sánchez-Lavega, A.; Lee, Y.J.; Muñoz, A.G.; Kouyama, T.; Sagawa, H.; Sato, T.M.; Piccioni, G.; Tellmann, S.; Imamura, T.; Satoh, T.  Stationary waves and slowly moving features in the night upper clouds of Venus.  Nat. Astron. 2017, 1, 0187, doi:10.1038/s41550-017-0187.
  * Peralta et al. (2018) Peralta, J.; Muto, K.; Hueso, R.; Horinouchi, T.; Sánchez-Lavega, A.; Murakami, S.Y.; Machado, P.; Young, E.F.; Lee, Y.J.; Kouyama, T.; et al.  Nightside Winds at the Lower Clouds of Venus with Akatsuki/IR2: Longitudinal, local time and decadal variations from comparison with previous measurements.  Astrophys. J. Suppl. Ser. 2018, 239, 17, doi:10.3847/1538-4365/aae844.
  * Peralta et al. (2019) Peralta, J.; Iwagami, N.; Sánchez-Lavega, A.; Lee, Y.J.; Hueso, R.; Narita, M.; Imamura, T.; Miles, P.; Wesley, A.; Kardasis, E.; Takagi, S.  Morphology and Dynamics of Venus’s Middle Clouds With Akatsuki/IR1.  Geophys. Res. Lett. 2019, 46, 2399–2407, doi:10.1029/2018GL081670.
  * Peralta et al. (2019) Peralta, J.; Sánchez-Lavega, A.; Horinouchi, T.; McGouldrick, K.; Garate-Lopez, I.; Young, E.F.; Bullock, M.A.; Lee, Y.J.; Imamura, T.; Satoh, T.; et al.  New cloud morphologies discovered on the Venus’s night during Akatsuki.  Icarus 2019, 333, 177–182, doi:10.1016/j.icarus.2019.05.026.
  * Peralta et al. (2020) Peralta, J.; Navarro, T.; Vun, C.W.; Sánchez-Lavega, A.; McGouldrick, K.; Horinouchi, T.; Imamura, T.; Hueso, R.; Boyd, J.P.; Schubert, G.; et al.  A Long-Lived Sharp Disruption on the Lower Clouds of Venus.  Geophys. Res. Lett. 2020, 47, e87221, doi:10.1029/2020GL087221.
  * Pérez-Hoyos et al. (2018) Pérez-Hoyos, S.; Sánchez-Lavega, A.; García-Muñoz, A.; Irwin, P.G.J.; Peralta, J.; Holsclaw, G.; McClintock, W.M.; Sanz-Requena, J.F.  Venus Upper Clouds and the UV Absorber From MESSENGER/MASCS Observations.  J. Geophys. Res. (Planets) 2018, 123, 145–162, doi:10.1002/2017JE005406.
  * Ross (1928) Ross, F.E.  Photographs of Venus.  Astrophys. J. 1928, 68, 57, doi:10.1086/143130.
  * Rossow et al. (1980) Rossow, W.B.; del Genio, A.D.; Limaye, S.S.; Travis, L.D.  Cloud morphology and motions from Pioneer Venus images.  J. Geophys. Res. 1980, 85, 8107–8128.
  * Sánchez-Lavega et al. (2016) Sánchez-Lavega, A.; Peralta, J.; Gomez-Forrellad, J.M.; Hueso, R.; Pérez-Hoyos, S.; Mendikoa, I.; Rojas, J.F.; Horinouchi, T.; Lee, Y.J.; Watanabe, S.  Venus Cloud Morphology and Motions from Ground-based Images at the Time of the Akatsuki Orbit Insertion.  Astrophys. J. Lett. 2016, 833, L7, doi:10.3847/2041-8205/833/1/L7.
  * Sánchez-Lavega et al. (2017) Sánchez-Lavega, A.; Lebonnois, S.; Imamura, T.; Read, P.; Luz, D.  The Atmospheric Dynamics of Venus.  Space Sci. Rev. 2017, 212, 1541–1616, doi:10.1007/s11214-017-0389-x.
  * Sato (2018) Sato, T.  A Brief History of ALPO-Japan.  J. Assoc. Lunar Planet. Obs. Stroll. Astron. 2018, 60, 20–21.
  * Schubert (1983) Schubert, G., General circulation and the dynamical state of the Venus atmosphere.  In Venus; Hunten, D.M., Colin, L., Donahue, T.M., Moroz, V.I., Eds.; University of Arizona Press: Tucson, AZ, USA, 1983; pp. 681–765.
  * Taguchi et al. (2007) Taguchi, M.; Fukuhara, T.; Imamura, T.; Nakamura, M.; Iwagami, N.; Ueno, M.; Suzuki, M.; Hashimoto, G.L.; Mitsuyama, K.  Longwave Infrared Camera onboard the Venus Climate Orbiter.  Adv. Space Res. 2007, 40, 861–868, doi:10.1016/j.asr.2007.05.085.
  * Titov et al. (2018) Titov, D.V.; Ignatiev, N.I.; McGouldrick, K.; Wilquet, V.; Wilson, C.F.  Clouds and Hazes of Venus.  Space Sci. Rev. 2018, 214, 126, doi:10.1007/s11214-018-0552-z.
  * Wood et al. (2022) Wood, B.E.; Hess, P.; Lustig-Yaeger, J.; Gallagher, B.; Korwan, D.; Rich, N.; Stenborg, G.; Thernisien, A.; Qadri, S.N.; Santiago, F.; Peralta, J.; et al.  Parker Solar Probe Imaging of the Night Side of Venus.  Geophys. Res. Lett. 2022, 49, e96302, doi:10.1029/2021GL096302.
  * Wright (1927) Wright, W.H.  Photographs of Venus made by Infra-red and by Violet Light.  Publ. Astron. Soc. Pac. 1927, 39, 220, doi:10.1086/123718.
  * Yamazaki et al. (2018) Yamazaki, A.; Yamada, M.; Lee, Y.J.; Watanabe, S.; Horinouchi, T.; Murakami, S.; Kouyama, T.; Ogohara, K.; Imamura, T.; Sato, T.M.; et al.  Ultraviolet imager on Venus orbiter Akatsuki and its initial results.  Earth Planets Space 2018, 70, 23, doi:10.1186/s40623-017-0772-6.


# Influence of atomic relaxations on the moiré flat band wavefunctions in
antiparallel twisted bilayer WS${}_{\text{2}}$

Laurent Molino  Department of Physics, University of Ottawa, Ottawa, Canada
Leena Aggarwal  Department of Physics, University of Ottawa, Ottawa, Canada
Indrajit Maity  Department of Materials, Imperial College London, and the
Thomas Young Centre for Theory and Simulation of Materials, London, UK  Ryan
Plumadore  Department of Physics, University of Ottawa, Ottawa, Canada
Johannes Lischner  Department of Materials, Imperial College London, and the
Thomas Young Centre for Theory and Simulation of Materials, London, UK  Adina
Luican-Mayer  Department of Physics, University of Ottawa, Ottawa, Canada
<EMAIL_ADDRESS>
###### Abstract

Twisting bilayers of transition metal dichalcogenides (TMDs) gives rise to a
periodic moiré potential resulting in flat electronic bands with localized
wavefunctions and enhanced correlation effects. In this work, scanning
tunneling microscopy is used to image a WS2 bilayer twisted approximately
$3^{\circ}$ off the antiparallel alignment. Scanning tunneling spectroscopy
reveals the presence of localized electronic states in the vicinity of the
valence band onset. In particular, the onset of the valence band is observed
to occur first in regions with a Bernal stacking in which S atoms are located
on top of each other. In contrast, density-functional theory calculations on
twisted bilayers which have been relaxed in vacuum predict the highest lying
flat valence band to be localized in regions of AA’ stacking. However,
agreement with the experiment is recovered when the calculations are carried
out on bilayers in which the atomic displacements from the unrelaxed positions
have been reduced reflecting the influence of the substrate and finite
temperature. This demonstrates the delicate interplay of atomic relaxations
and the electronic structure of twisted bilayer materials.

*These authors contributed equally

When a pair of two-dimensional (2D) materials are stacked vertically and then
rotated relative to each other, a moiré superlattice is created. Electrons in
such a superlattice experience a moiré potential which favors electron
localization in specific stacking regions, reduces the electrons’ kinetic
energy, and results in the formation of flat bands with enhanced electron-
electron interactions, promoting phenomena such as superconductivity, Mott
insulating phases, topological phases, coupled spin-valley states, and
interlayer excitons.1, 2, 3, 4, 5, 6, 7

In contrast to graphene bilayers which exhibit flat bands at a specific
“magic” twist angle, bilayers of transition metal dichalcogenides (TMDs) offer
a much larger range of twist angles at which flat bands are realized8, 9, 10,
11, 12, 13, 14, opening broader possibilities for creating previously
inaccessible correlated systems. Both TMD homobilayers (in which two layers of
the same material are combined) and heterobilayers (in which two different TMD
monolayers are combined) have been studied. In heterobilayers, the length
scale of the moiré pattern is less sensitive to the twist-angle than in
homobilayers, leading to reduced twist angle disorder in large samples. As a
result, to date, several heterobilayer TMDs such as WSe2/WS2 11, 15, MoS2/WS2
16, 17, 18, 19, MoS2/WSe220, 13 and WSe2/MoSe2 21 have been investigated
experimentally, with reports of correlated phenomena such as
superconductivity, Mott insulating phases, generalized Wigner crystals,
etc.22, 9, 23, 11. Moiré patterns of TMD homobilayers are highly sensitive to
the twist angle, and thus expected to exhibit more twist-angle disorder. As a
consequence, these systems are less studied despite having great potential for
inducing strong moiré potentials 8, 24, 3, 30, 26. For example, recent
scanning probe microscopy experiments on a WSe2/WSe2 homobilayer have
demonstrated twist-angle-dependent lattice reconstructions near a $60^{\circ}$
twist angle resulting in multiple energy-separated ultra-flat valence bands
for long moiré wavelengths (6 nm-12 nm).10

Semiconducting TMDs have metal and chalcogen atoms occupying the lattice sites
in their unit cell breaking sub-lattice symmetry. For this reason, twisted TMD
bilayers fall into two classes: parallel (P) aligned bilayers with a zero
twist angle close to $0^{\circ}$, and antiparallel (AP) aligned bilayers with
twist angles near $60^{\circ}$. These two classes exhibit different stacking
configurations in the moiré unit cell which results in different atomic
relaxations in the form of both out-of-plane buckling and in-plane
reconstructions. In-plane reconstructions can be understood as the result of
the layers deforming to increase the area of the energetically preferred
stacking order. 27, 24, 28, 8 For very small twist angles (or twist angles
very close to $60^{\circ}$), this reconstruction produces networks of uniform
stacking order domains separated by highly strained domain walls. In the case
of P-aligned TMD homobilayers, the system tends towards a network of
triangular stacking order domains. Due to the lack of inversion symmetry,
these domains host interlayer charge transfer leading to an out-of-plane
ferroelectric dipole moment.24, 21, 27, 29 In the case of AP-aligned TMD
homobilayers, no such charge transfer occurs. While AP systems do form domain
wall networks at angles very close to $60^{\circ}$, they generally exhibit
less pronounced in-plane reconstructions than P-aligned bilayers at a given
twist-angle offset, leading to a more smoothly varying periodic lattice
distortion over the moiré unit cell. 25 In all cases, both out-of-plane
buckling and in-plane reconstruction have a profound impact on the electronic
properties of TMD homobilayers. In-plane and out-of-plane deformations cause
the simple band-folding description of a moiré pattern band-structure to break
down, for example, causing a modification in the real-space position of the
valence band edge wavefunctions. Strain is furthermore critical to
wavefunction localization responsible for strong correlated effects observed
in these systems. 8, 24

As the field of moiré-induced flat band physics progresses, it is important to
explore a larger range of homobilayer systems, especially using scanning probe
techniques that can directly reveal the local twist angle. Here, we report
room temperature scanning tunnelling microscopy/spectroscopy (STM/STS) studies
of moiré superlattices of $\approx$ 6 nm wavelength in the WS2/WS2 bilayer
system near $60^{\circ}$ twist angle or antiparallel (AP) alignment, together
with extensive large-scale density-functional theory calculations.
Importantly, we demonstrate that the localization of the flat-band
wavefunctions depends sensitively on the degree of lattice relaxations.

Figure 1:  (a) Schematic of the STM experiment: a twisted bilayer WS2 with
$\theta\approx 60^{\circ}$ is prepared on a graphite/Au substrate. (b) STM
topograph showing a region with 6.0 nm ± 0.2 nm moiré pattern, corresponding
to a twist angle of $57.1^{\circ}\pm 0.1^{\circ}$ (scanning parameters: bias
voltage $V_{b}=-1.3V$, tunneling current $I_{t}=100$ pA). Inset: Fourier
transform of the moiré pattern showing hexagonal symmetry. (c) STM topograph
of 6.0 nm ± 0.2 nm moiré pattern (scanning parameters: $V_{b}=-1.3$ V,
$I_{t}=150~{}$pA). The dashed rhombus corresponds to the moiré unit cell, with
high-symmetry stacking regions ($\mathrm{B^{S/S}}$, $\mathrm{AA}$’,
$\mathrm{B^{W/W}}$) indicated by circles. (d) Atomic structure of the three
high-symmetry stackings present in an antiparallel aligned twisted WS2
homobilayer. (e) Simulated STM topography for a $56.86^{\circ}$ twisted AP-WS2
homobilayer obtained from first-principles density-functional theory
calculations on a relaxed atomic structure. (f) Side view of twisted bilayer
WS2 showing out-of-plane buckling, exaggerated by a factor of ten for clarity.
Locations of maximum and minimum interlayer spacing are indicated. The maximum
interlayer spacing, corresponding to B${}^{\text{S/S}}$ stacking, is 8% larger
than the minimum interlayer spacing, corresponding to AA’ stacking.

We fabricated the AP-WS2 twisted bilayer using a variation of the “tear-and-
stack” dry transfer method (see Methods section) and performed room
temperature STM/STS experiments, as indicated schematically in Figure 1(a).
Across the entire sample, we identified moiré patterns with a range of periods
(6 nm-14 nm) (Supplementary Figure S1), corresponding to a twist-angle range
of $57^{\circ}$-$58.8^{\circ}$, indicating the presence of twist-angle
disorder within the AP-WS2 sample. A large area of uniform hexagonal moiré
pattern with a period of 6.0 nm $\pm$ 0.2 nm (corresponding to a twist angle
of $57.1^{\circ}\pm 0.1^{\circ}$) is shown in Figure 1(b). A zoomed-in
scanning tunnelling micrograph of this pattern in Figure 1(c) highlights the
unit cell and the different high-symmetry stacking regions within it
($\mathrm{B^{S/S}}$, $\mathrm{AA}$’, $\mathrm{B^{W/W}}$), with corresponding
atomic configurations described in Figure 1(d).

A simulated STM topograph of a $56.86^{\circ}$ AP-twisted WS2 bilayer obtained
from density functional theory (DFT) calculations based on the Tersoff-Hamann
approximation 31 is shown in Figure 1(e). Before the DFT calculations, the
relaxed atomic positions of the twisted bilayer, shown in Fig. 1(f), were
obtained using classical force fields (see Methods). The simulated topograph
is in good agreement with the experimental micrograph of Figure 1(c). In
particular, the bright regions correspond to $\mathrm{B^{S/S}}$ stacking. This
stacking exhibits the largest buckling and therefore also the largest
tunnelling matrix element.

Figure 2:  (a) Average of scanning tunneling spectra acquired over the moiré
region indicated in Supplementary Figure S3.(b) Electronic band structure of
$56.86^{\circ}$ twisted bilayer WS2 including atomic relaxations. The valence
band maximum is set to 0. (c) Calculated density of states for the same
relaxed twisted bilayer WS2 as in (b).

To explore the local electronic properties of the moiré superlattice in Figure
1(b), we performed scanning tunnelling spectroscopy (STS) measurements. In
Figure 2(a) we present the average tunnelling conductance (dI/dV) acquired
over five moiré periods (Supplementary Figures S2 and S3). Two notable
features are present in the spectrum: a first one (denoted V1) near the
valence band edge (VBE), with an energy of $\approx-1.0$ eV, and a second
feature V2 deeper in the valence band at $\approx-1.7$ eV.

These two features are also found in the density of states (DOS) calculated
using DFT for the relaxed atomic structure, see Fig. 2(c). The energy
separation between the two features is found to be 0.7 eV, in excellent
agreement with the experimental result. We note that the valence band edge of
the calculated DOS is much sharper than in the experiment, which is likely a
consequence of thermal broadening.

To understand the origin of the observed features close to the valence band
edge, we also present the DFT band structure of the twisted bilayer in Fig.
2(b). This reveals that the highest valence band is extremely flat with a
bandwidth of 2 meV. Approximately 60 meV below the highest valence band we
find another pair of relatively flat bands which contribute to the sharp
feature $V_{1}$. We also identify another sharp feature in the DOS 0.7 eV
below the $V_{1}$, in good agreement with the experiment.

Figure 3: (a),(c),(e) Scanning tunneling spectra acquired at high symmetry
sites in three areas having different twist angles as indicated. Inset:
Topographic maps indicating the points where the data were taken. (b),(d),(f)
dI/dV maps acquired near the valence band edge at bias voltages indicated
(tunneling current $I_{t}=100pA$). The high symmetry sites and the moiré unit
cell are indicated in each map.

Taking advantage of the spatial resolution of a STM/STS experiment, we
measured dI/dV maps corresponding to the feature V1. In Figure 3 we present
scanning tunnelling spectra collected at AA’ and $\mathrm{B^{S/S}}$ stacking
regions taken on three areas with slightly different twist angles:
$57.0^{\circ}\pm 0.3^{\circ}$ [Figure 3(a),(b)], $57.1^{\circ}\pm 0.1^{\circ}$
[Figure 3(c),(d)] and $57.4^{\circ}\pm 0.2^{\circ}$ [Figure 3(e),(f)]. The
insets show the associated topographic image including the high-symmetry
points and unit cell. These maps demonstrate that the feature V1 is associated
with localized states that form a hexagonal pattern with low intensity on the
$\mathrm{B^{S/S}}$ sites and high intensity on the AA’ sites. The spectra
taken at the AA’ and $\mathrm{B^{S/S}}$ sites show that the valence band onset
occurs at the $\mathrm{B^{S/S}}$ site, suggesting that the highest-lying
valence band is localized at this site (indicating that the localized state
associated with V1 is _not_ the highest lying valence band). While this
difference between the two sites is clearly visible for the $57.1^{\circ}\pm
0.1^{\circ}$ system, the separation between the two spectra is much smaller
for the other two systems (for the $57.0^{\circ}\pm 0.3^{\circ}$ system the
onsets at the AA’ and the $\mathrm{B^{S/S}}$ sites occur almost at the same
energy). We also note that more significant differences between the various
twist-angle regions are found further away from the valence band onset
indicating that the electronic structure is highly sensitive to the local
twist angle.

Figure 4: Calculated Local Density Of States (LDOS) from DFT for a free-
standing twisted bilayer (a) including the full atomic relaxations, (b) with
all atomic displacements reduced by a factor of 2 compared to the fully
relaxed result, (c) with all atomic displacements reduced by a factor of 4
compared to the fully relaxed result and (d) without any atomic relaxations.
The valence band maximum is set to 0. The insets show the calculated STM
topographs. The associated color bars indicate variation in height in nm.
(e)-(h): Squared wavefunctions of the highest valence band at the
$\Gamma$-point of the moiré Brillouin zone for the atomic structures shown in
the top panel. (i)-(l): Squared wavefunctions of a deeper lying valence band
at the $\Gamma$-point of the moiré Brillouin zone for the atomic structures
shown in the top panel.

To further understand the experimental findings, we calculated the local
density of states (LDOS) at the AA’ and $\mathrm{B^{S/S}}$ sites using DFT
(see SI for more details). Figure 4(a) shows the result for a fully relaxed
free-standing twisted bilayer (i.e. a bilayer in free space without a
substrate). In the inset of Fig. 4(a), we also show the calculated STM
topography. The LDOS in Figure 4(a) suggests that in this scenario the highest
valence band is highly localized at the $\mathrm{AA}$’ sites, contrary to the
experimental findings of Figure 3. Interestingly, the highest valence band is
localized on the $\mathrm{B^{S/S}}$ site when no atomic relaxations are
included and the layers remain flat, see Fig. 4(d). This suggests that atomic
relaxations in the experimental sample are not as strong as in a free-standing
system. Possible origins of this discrepancy include (i) the influence of the
substrate which is not included in our simulations, (ii) atomic motions due to
thermal fluctuations, and (iii) extrinsic strains. To test this hypothesis, we
also carry out DFT calculations for twisted bilayers in which the
displacements of all atoms from their unrelaxed positions are reduced by a
factor of 2 in Fig. 4(b) and by a factor of 4 in Fig. 4(c), compared to their
values in the fully relaxed free-standing system. In Figure 4(h) we show the
calculation in the absence of relaxation. In these systems, the LDOS in the
AA’ and the $\mathrm{B^{S/S}}$ regions are almost the same near the valence
band onset in good agreement with the experimental findings. These results
clearly demonstrate the high sensitivity of the electronic structure on atomic
relaxations.

Figures 4(e)-(h) demonstrate that the wavefunction of the highest lying
valence band does not correspond to the experimental measurement for the
feature V1 shown in the lower panel of Fig. 3: the wavefunctions do not form a
honeycomb shape, but a triangular pattern. The only exception corresponds to
the wavefunctions obtained when all displacements are reduced by a factor of
4, see Fig. 4(g). However, in this case, the dark regions of the honeycomb
pattern are not associated with $\mathrm{B^{S/S}}$ regions as in the
experiment. Instead, we have found that a deeper-lying valence band has a low
intensity on the Bs/S regions when relaxations are absent or strongly reduced,
see Figs. 4(g) and (h). This state also exhibits a six-fold symmetry around
these BS/S sites as does the honeycomb pattern observed in experiments. We,
therefore, assign the experimental maps of the feature V1 to this deeper-lying
valence band. Considering the thermal broadening in this experiment, it is not
surprising that the contributions from a single narrow flat band would not be
the only ones measured at a particular energy.

In conclusion, we use scanning tunneling microscopy and spectroscopy to
characterize the spatial distribution of the wave functions corresponding to
flat bands in a bilayer WS2 sample twisted $\approx 3^{\circ}$ off the
antiparallel alignment. To understand the localization at the different high-
symmetry stacking regions we performed density-functional theory calculations
which revealed that atomic relaxations play a critical role in determining the
degree to which the wave functions associated with flat bands localize on the
high-symmetry regions. We find that in our samples the atomic relaxations are
much suppressed compared to calculations for the idealized system of a bilayer
in vaccum. These findings provide insight into the spatial distributions of
the wave functions associated with flat bands believed to be responsible for
the occurrence of correlated electron phenomena, demonstrating that atomic
relaxations are a way to control the localization of the wave functions.

## References

  * 1 Shimazaki, Y.; Schwartz, I.; Watanabe, K.; Taniguchi, T.; Kroner, M.; Imamoğlu, A. Strongly correlated electrons and hybrid excitons in a moiré heterostructure. Nature 2020, 580, 472–477.
  * 2 Yankowitz, M.; Chen, S.; Polshyn, H.; Zhang, Y.; Watanabe, K.; Taniguchi, T.; Graf, D.; Young, A. F.; Dean, C. R. Tuning Superconductivity in twisted bilayer graphene. Science 2019, 363, 1059–1064.
  * 3 An, L.; Cai, X.; Pei, D.; Huang, M.; Wu, Z.; Zhou, Z.; Lin, J.; Ying, Z.; Ye, Z.; Feng, X., et al. Interaction effects and superconductivity signatures in twisted double-bilayer WSe2. Nanoscale Horiz. 2020, 5, 1309–1316.
  * 4 Brotons-Gisbert, M.; Baek, H.; Molina-Sánchez, A.; Campbell, A.; Scerri, E.; White, D.; Watanabe, K.; Taniguchi, T.; Bonato, C.; Gerardot, B. D. Spin–layer locking of inter-layer excitons trapped in moiré potentials. Nat. Mater. 2020, 19, 630–636.
  * 5 Choi, J.; Hsu, W.-T.; Lu, L.-S.; Sun, L.; Cheng, H.-Y.; Lee, M.-H.; Quan, J.; Tran, K.; Wang, C.-Y.; Staab, M., et al. Moiré potential impedes interlayer exciton diffusion in van der Waals heterostructures. Sci. Adv. 2020, 6, eaba8866.
  * 6 Tran, K.; Moody, G.; Wu, F.; Lu, X.; Choi, J.; Kim, K.; Rai, A.; Sanchez, D. A.; Quan, J.; Singh, A., et al. Evidence for moiré excitons in van der Waals heterostructures. Nature 2019, 567, 71–75.
  * 7 Das, S.; Dandu, M.; Gupta, G.; Murali, K.; Abraham, N.; Kallatt, S.; Watanabe, K.; Taniguchi, T.; Majumdar, K. Highly tunable layered exciton in bilayer WS2: linear quantum confined stark effect versus electrostatic doping. ACS Photonics 2020, 7, 3386–3393.
  * 8 Naik, M. H.; Kundu, S.; Maity, I.; Jain, M. Origin and evolution of ultraflat bands in twisted bilayer transition metal dichalcogenides: Realization of triangular quantum dots. Phys. Rev. B 2020, 102, 075413.
  * 9 Wang, L.; Shih, E.-M.; Ghiotto, A.; Xian, L.; Rhodes, D. A.; Tan, C.; Claassen, M.; Kennes, D. M.; Bai, Y.; Kim, B., et al. Correlated electronic phases in twisted bilayer transition metal dichalcogenides. Nat. Mater. 2020, 19, 861–866.
  * 10 Li, E.; Hu, J.-X.; Feng, X.; Zhou, Z.; An, L.; Law, K. T.; Wang, N.; Lin, N. Lattice reconstruction induced multiple ultra-flat bands in twisted bilayer WSe2. Nat. Commun. 2021, 12, 1–7.
  * 11 Li, H.; Li, S.; Naik, M. H.; Xie, J.; Li, X.; Wang, J.; Regan, E.; Wang, D.; Zhao, W.; Zhao, S., et al. Imaging moiré flat bands in three-dimensional reconstructed WSe2/WS2 superlattices. Nat. Mater. 2021, 20, 945–950.
  * 12 Vitale, V.; Atalar, K.; Mostofi, A. A.; Lischner, J. Flat band properties of twisted transition metal dichalcogenide homo- and heterobilayers of MoS2, MoSe2, WS2 and WSe2. 2d Mater. 2021, 8, 045010.
  * 13 Waters, D.; Nie, Y.; Lupke, F.; Pan, Y.; Folsch, S.; Lin, Y.-C.; Jariwala, B.; Zhang, K.; Wang, C.; Lv, H., et al. Flat bands and mechanical deformation effects in the moiré superlattice of MoS2-WSe2 heterobilayers. ACS Nano 2020, 14, 7564–7573.
  * 14 Maity, I.; Maiti, P. K.; Krishnamurthy, H. R.; Jain, M. Reconstruction of moiré lattices in twisted transition metal dichalcogenide bilayers. Phys. Rev. B 2021, 103, L121102.
  * 15 Li, H.; Li, S.; Naik, M. H.; Xie, J.; Li, X.; Regan, E.; Wang, D.; Zhao, W.; Yumigeta, K.; Blei, M., et al. Imaging local discharge cascades for correlated electrons in WS2/WSe2 moiré superlattices. Nat. Phys. 2021, 17, 1114–1119.
  * 16 Kobayashi, Y.; Yoshida, S.; Sakurada, R.; Takashima, K.; Yamamoto, T.; Saito, T.; Konabe, S.; Taniguchi, T.; Watanabe, K.; Maniwa, Y., et al. Modulation of electrical potential and conductivity in an atomic-layer semiconductor heterojunction. Sci. Rep. 2016, 6, 31223.
  * 17 Shi, J.; Tong, R.; Zhou, X.; Gong, Y.; Zhang, Z.; Ji, Q.; Zhang, Y.; Fang, Q.; Gu, L.; Wang, X.; Liu, Z.; Zhang, Y. Temperature-Mediated Selective Growth of MoS2/WS2 and WS2/MoS2 Vertical stacks on Au foils for direct photocatalytic applications. Adv. Mater. 2016, 28, 10664–10672.
  * 18 Zhang, F.; Lu, Z.; Choi, Y.; Liu, H.; Zheng, H.; Xie, L.; Park, K.; Jiao, L.; Tao, C. Atomically resolved observation of continuous interfaces between an as-grown MoS2 monolayer and a WS2/MoS2 heterobilayer on SiO2. ACS Appl. Nano Mater. 2018, 1, 2041–2048.
  * 19 Hill, H. M.; Rigosi, A. F.; Rim, K. T.; Flynn, G. W.; Heinz, T. F. Band alignment in MoS2/WS2 transition metal dichalcogenide heterostructures probed by scanning tunneling microscopy and spectroscopy. Nano Lett. 2016, 16, 4831–4837.
  * 20 Pan, Y.; Folsch, S.; Nie, Y.; Waters, D.; Lin, Y.-C.; Jariwala, B.; Zhang, K.; Cho, K.; Robinson, J. A.; Feenstra, R. M. Quantum-confined electronic states arising from the moiré pattern of MoS2–WSe2 heterobilayers. Nano Lett. 2018, 18, 1849–1855.
  * 21 Shabani, S.; Halbertal, D.; Wu, W.; Chen, M.; Liu, S.; Hone, J.; Yao, W.; Basov, D. N.; Zhu, X.; Pasupathy, A. N. Deep moiré potentials in twisted transition metal dichalcogenide bilayers. Nat. Phys. 2021, 17, 720–725.
  * 22 Regan, E. C.; Wang, D.; Jin, C.; Bakti Utama, M. I.; Gao, B.; Wei, X.; Zhao, S.; Zhao, W.; Zhang, Z.; Yumigeta, K., et al. Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices. Nature 2020, 579, 359–363.
  * 23 Tang, Y.; Li, L.; Li, T.; Xu, Y.; Liu, S.; Barmak, K.; Watanabe, K.; Taniguchi, T.; MacDonald, A. H.; Shan, J., et al. Simulation of Hubbard model physics in WSe2/WS2 moiré superlattices. Nature 2020, 579, 353–358.
  * 24 Naik, M. H.; Jain, M. Ultraflatbands and shear solitons in moiré patterns of twisted bilayer transition metal dichalcogenides. Phys. Rev. Lett. 2018, 121, 266401.
  * 25 Van Winkle, M.; Craig, I. M.; Carr, S.; Dandu, M.; Bustillo, K. C.; Ciston, J.; Ophus, C.; Taniguchi, T.; Watanabe, K.; Raja, A., et al. Quantitative imaging of intrinsic and extrinsic strain in transition metal dichalcogenide moiré bilayers. ArXiv:2212.07006 2022.
  * 26 Lin, J.; Han, T.; Piot, B. A.; Wu, Z.; Xu, S.; Long, G.; An, L.; Cheung, P.; Zheng, P.P.; Plochocka, P., et al. Determining interaction enhanced valley susceptibility in spin-valley-locked MoS2. Nano Lett. 2019, 19, 1736–1742.
  * 27 Weston, A.; Zou, Y.; Enaldiev, V.; Summerfield, A.; Clark, N.; Zólyomi, V.; Graham, A.; Yelgel, C.; Magorrian, S.; Zhou, M., et al. Atomic reconstruction in twisted bilayers of transition metal dichalcogenides. Nat. Nanotechnol. 2020, 15, 592–597.
  * 28 Enaldiev, V. V.; Zólyomi, V.; Yelgel, C.; Magorrian, S. J.; Fal’ko, V. I. Stacking domains and dislocation networks in marginally twisted bilayers of transition metal dichalcogenides. Phys. Rev. Lett. 2020, 124, 206101.
  * 29 Molino, L.; Aggarwal, L.; Enaldiev, V.; Plumadore, R.; Fal’ko, V.; Luican-Mayer, A. Ferroelectric switching at symmetry-broken interfaces by local control of dislocation networks. ArXiv:2210.03074 2022.
  * 30 Winkle, M.; Craig, I.; Carr, S.; Dandu, M.; Bustillo, K.; Ciston, J.; Ophus, C.; Taniguchi, T.; Watanabe, K.; Raja, A.; Griffin, S.; Bediako, D. Quantitative imaging of intrinsic and extrinsic strain in transition metal dichalcogenide moiré bilayers. ArXiv:2212.07006 2022.
  * 31 Tersoff, J.; Hamann, D. R. Theory of the scanning tunneling microscope. Phys. Rev. B 1985, 31, 805–813.
  * 32 Boddison-Chouinard, J.; Plumadore, R.; Luican-Mayer, A. Fabricating van der Waals heterostructures with precise rotational alignment. JoVE 2019, 149.
  * 33 Boddison-Chouinard, J.; Scarfe, S.; Watanabe, K.; Taniguchi, T.; Luican-Mayer, A. Flattening van der Waals heterostructure interfaces by local thermal treatment. Appl. Phys. Lett. 2019, 115, 231603.
  * 34 Lindvall, N.; Kalabukhov, A.; Yurgens, A. Cleaning graphene using atomic force microscope. J. Appl. Phys. 2012, 111, 064904.
  * 35 Rosenberger, M. R.; Chuang, H.-J.; McCreary, K. M.; Hanbicki, A. T.; Sivaram, S. V.; Jonker, B. T. Nano-“squeegee” for the creation of clean 2D material interfaces. ACS Appl. Mater. Interfaces. 2018, 10, 10379–10387.
  * 36 Naik, S.; Naik, M. H.; Maity, I.; Jain, M. Twister: Construction and structural relaxation of commensurate moiré superlattices. Comput. Phys. Commun. 2022, 271, 108184.
  * 37 Jiang, J.-W.; Zhou, Y.-P. Handbook of Stillinger-Weber potential parameters for two-dimensional atomic crystals; IntechOpen: Rijeka, 2017.
  * 38 Naik, M. H.; Maity, I.; Maiti, P. K.; Jain, M. Kolmogorov–Crespi potential for multilayer transition-metal dichalcogenides: capturing structural transformations in moiré superlattices. J. Phys. Chem. C 2019, 123, 9770–9778.
  * 39 Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in’t Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comp. Phys. Comm. 2022, 271, 108171.
  * 40 https://lammps.org/. https://lammps.org/.
  * 41 Soler, J. M.; Artacho, E.; Gale, J. D.; García, A.; Junquera, J.; Ordejón, P.; Sánchez-Portal, D. The SIESTA method for ab initio order-N materials simulation. J. Condens. Matter Phys. 2002, 14, 2745–2779.
  * 42 Fernández-Seivane, L.; Oliveira, M. A.; Sanvito, S.; Ferrer, J. On-site approximation for spin–orbit coupling in linear combination of atomic orbitals density functional methods. J. Condens. Matter Phys. 2006, 18, 7999–8013.
  * 43 Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43, 1993–2006.
  * 44 Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B 1981, 23, 5048–5079.

## Methods

The sample was prepared by the dry transfer method. Two monolayers of WS2
(purchased from HQ Graphene) were stacked on each other using the tear-and-
stack technique 32, 33. Graphene (purchased from Graphenea) was used as a
substrate to allow a sufficiently conductive contact for STM measurements.
Both WS2 and graphite crystals were mechanically exfoliated onto Si/SiO2
substrates. Monolayer WS2 was identified optically, and the thickness was
confirmed with atomic force microscopy (AFM). Suitable graphite flakes were
identified and were picked up with a polydimethylsiloxane (PDMS) stamp covered
with Polypropylene carbonate (PPC) film on a glass slide. The graphite on the
polymer stamp was then used to pick up monolayer WS2. The two layers of WS2
correspond to two halves of a large flake, allowing for precise rotational
control. The PPC film was then cut and the 2D material stack was transferred
onto a Si/SiO2 substrate coated with Ti/Au (5 nm/150 nm). The sample was then
annealed under ultra-high vacuum at 300°C for several hours. The surface was
cleaned using an AFM tip 34, 35. STM measurements were performed at pressures
below $10^{-9}$ Torr, at room temperature. For STS we used bias modulation of
5.0 mV, 1.325 kHz and a typical set-point: $V_{B}=-1.3V$, $I_{T}=150pA$.

## Computational Details

##  0.1 Atomic reconstructions

The twisted bilayer of WS2 has been generated using the TWISTER package 36.
The total number of atoms in the unit cell of $56.86^{\circ}$ twisted bilayer
WS2 is 1986. All the atomic relaxations are performed using classical
interatomic potentials fitted to density functional theory calculations. The
intralayer interactions within WS2 are described using a Stillinger- Weber
potential 37, and the interlayer interactions are captured using a Kolmogorov-
Crespi potential 38. All the simulations with classical interatomic potentials
are performed using the LAMMPS package 39, 40. We relax the atoms within a
fixed simulation box with a force tolerance of 10-6 eV/Å for any atom along
any direction.

##  0.2 Electronic structure calculations

All the electronic structure calculations were performed using the SIESTA
package 41. All the calculations included spin-orbit coupling 42. We used
Troullier-Martins pseudopotential 43 and the local density approximation with
the Perdew- Zunger parametrization as the exchange-correlation functional 44.
We used a double-$\zeta$ plus polarization basis for the expansion of
wavefunctions. For the electronic structure calculations we used the $\Gamma$
point in the moiré Brillouin zone to obtain the converged ground state charge
density. A large vacuum spacing of $20$ Å was used in the out-of-plane
direction.

###  0.2.1 Density of states calculations

To obtain well-converged results, the density of states calculations were
performed on a fine $5\times 5\times 1$ k-grid in the moiré BZ (which
corresponds to a $91\times 91\times 1$ k-grid of the monolayer Brillouin
zone). We used a Gaussian with a width of 30 meV to represent the
$\delta$-function in the density of states calculations.

###  0.2.2 Local density of states (LDOS) calculations

To obtain the LDOS for selected stackings, we used a $5\times 5\times 1$
k-grid to sample the moiré Brillouin zone. Furthermore, we average the LDOS
over a 9 Å radius in-plane circular patch around the centre of the selected
high-symmetry stackings.

###  0.2.3 STM topography calculations

The calculations associated with the STM topography were performed only using
the $\Gamma$ point in the moiré Brillouin zone. Furthermore, we used a Voronoi
diagram to coarse-grain atomistic details in the topography.

The codes developed for this work (selected LDOS and STM topography) will be
made publicly available at
https://gitlab.com/_imaity_/Moire_toolkit_electronic.

## Acknowledgement

AL-M, LM, RP, LA acknowledge funding from NSERC Discovery RGPIN-2022-05215,
NRC Quantum Sensing Challenge, Ontario Early Researcher Award ER-16-218. This
project received funding from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie Grant agreement No.
101028468. This work used the ARCHER2 UK National Supercomputing Service
(https://www.archer2.ac.uk) via our membership of the UK’s HEC Materials
Chemistry Consortium, which is funded by EPSRC (EP/R029431).

## Supplementary information for: Moiré flat band in antiparallel twisted
bilayer WS2 relevaled by scanning tunneling microscopy

## Section I: Different moiré patterns observed over the sample

Figure S1: (a,b,c,d) STM topography with four different twist angles as
indicated, corresponding to moiré wavelengths of 14 nm ± 1 nm, 13 nm ± 3 nm,
10.1 nm ± 0.7 nm, 7.9 nm ± 0.4 nm respectively. Data were taken at
$V_{B}=-1.3V$, $I_{t}=100pA$.

## Section II: Spatial variation of the local density of states over moiré
periods

Figure S2: STS line profile crossing 4.5 moiré periods, demonstrating moiré
modulation around $V_{1}=-0.9V$ and $V_{2}=-1.6V$. The above figure indicates
the height profile along the line indicated in the STM topography in Figure
S3(a).

## Section III: Spectroscopic gap calculated from the average spectra

Figure S3: (a) Area in which Figure 2 (a) spectra where collected. Spectra in
Figure 2 a) is the average of the curves collected along the indicated line.
b) bang-gap estimation of scanning tunneling spectroscopy curve corresponding
to Figure 2 (a). In a log plot, linear region of the valence and conduction
band edges and the band gap are identified, and used to fit lines. The
intersection of these lines correspond to the valence and conduction band
edges as appropriate.