2022

[1,2]Yueshen Wu These authors contributed equally to this work.

These authors contributed equally to this work.

[1,2,3]Yanpeng Qi

[1,2,6,7]Jun Li

[1]School of Physical Science and Technology, ShanghaiTech University,
Shanghai, 201210, China

[2]ShanghaiTech Laboratory for Topological Physics, ShanghaiTech University,
Shanghai, 201210, China

3]Shanghai Key Laboratory of High-resolution Electron Microscopy, ShanghaiTech
University, Shanghai, 201210, China

4]Research Center for Functional Materials, National Institute for Materials
Science, Tsukuba, 305-0044, Japan

5]International Center for Materials Nanoarchitectonics, National Institute
for Materials Science, Tsukuba, 305-0044, Japan

6]State Key Laboratory of Functional Materials for Informatics, Shanghai
Institute of Microsystem and Information Technology, Chinese Academy of
Sciences, 865 Changning Road, Shanghai, 200050, China

7]Wuhan National High Magnetic Field Center, Huazhong University of Science &
Technology, Wuhan, 430074, China

# Nonreciprocal Charge Transport in Topological Kagome Superconductor CsV3Sb5

<EMAIL_ADDRESS>Qi Wang  Xiang Zhou  Jinghui Wang  Peng Dong
Jiadian He  Yifan Ding  Bolun Teng  Yiwen Zhang  Yifei Li  Chenglong Zhao
Hongti Zhang  Jianpeng Liu<EMAIL_ADDRESS>Kenji Watanabe  Takashi
Taniguchi<EMAIL_ADDRESS>*  *  [  [  [  [  [

###### Abstract

Nonreciprocal charge transport phenomena are widely studied in two-dimensional
superconductors, which demonstrate unidirectional-anisotropy
magnetoresistances as a result of symmetry breaking. Here, we report a strong
nonreciprocal transport phenomenon in superconducting CsV3Sb5 thin flakes. The
second harmonic voltages, mainly originating from the rectification effect of
vortex motion, are unambiguously developed with in-plane and out-of-plane
magnetic fields, and their magnitudes are comparable to those in
noncentrosymmetric superconductors. The second harmonic magnetoresistances
split into several peaks and some of them reverse their signs by ramping the
magnetic field or the current within the superconducting transition. The
nonreciprocity suggests a strong asymmetry in CsV3Sb5. The centrosymmetric
structure and symmetric electronic phases in CsV3Sb5 can hardly induce the
distinct nonreciprocal transport phenomenon, which could be correlated to a
symmetry breaking from an unconventional superconducting order parameter
symmetry.

## Introduction

The nonreciprocal charge transport phenomenon, also called magnetochiral
anisotropy rikken1997observation  or unidirectional magnetoresistance
zhou2021sign , results from inversion and time-reversal symmetry breaking that
can modulate the energy dispersion with unequal electron energies of opposite
momentum k and $-$k or spin $s$ and $-s$ Tokura . Based on this phenomenon, a
wide range of potential applications such as rectifiers, alternating-direct-
current converters, and photodetectors could be engineered Villegas1188 ;
moshchalkov2010nanoscience ; MCAZhang ; nagulu2020non . In superconductors,
the nonreciprocal signal can be considerably enhanced in the superconducting
state by several orders of magnitude due to the significant difference between
Fermi energy and superconducting energy gaps Liueabc7628 ; PhysRevB.98.054510
; daido2022intrinsic , for which the nonreciprocal transport is related to the
chirality of supercurrent or Cooper pairs. Up to now, nonreciprocal transport
measurements have been performed in various superconductors such as the
noncentrosymmetric MoS2 MoS2PRR ; Liueabc7628 , NbSe2 MCAZhang , Nb/V/Ta
superlattice MCAAndo , the polar superconductor SrTiO3 Itahashieaay9120 , and
the artificially engineered MoGe superconducting thin films
lyu2021superconducting ,all of which are basically related to the symmetry
breaking in structures, electronic states, or even anisotropic friction forces
on vortex flowing. On the other hand, the intrinsic vortex motion behavior
correlated with the superconducting gap symmetry can also contribute to the
nonreciprocity Liueabc7628 . In addition, such nonreciprocal transport
phenomenon has been observed in the topological superconducting
heterostructure of Bi2Te3/FeTe MCAYasuda . Thus, these nonreciprocal responses
can be studied as a probe to detect the chirality in an anisotropic
superconducting states.

Recently, a series of kagome metals $A$V3Sb5($A$ = K, Rb, Cs)
PhysRevLett.125.247002  have attracted increasing attentions with the unique
coexistence of topological states, charge density wave (CDW), and multi-band
superconductivity. Among these crystals, CsV3Sb5 reveals the highest of 2.5 K
PhysRevLett.127.177001 . Up to now, fruitful experiments have been actualized
to understand the superconducting gap symmetry of CsV3Sb5, including thermal
conductivity zhao2021nodal , scanning tunneling spectroscopy measurements
xu2021multiband , and penetration length measurements Nodeless . However,
these results provided contradictory conclusions on the gap symmetry
PhysRevLett.127.177001 ; ding2022diagnosis ; kim2022monolayer ;
jiang2022electronic ; Mu_2021_s_wave ; zhao2021nodal ; xu2021multiband ;
Nodeless ; Shunli ; RN985 ; RN986 ; chiral ; PhysRevLett.127.046401 ;
feng2021chiral ; yu2021evidence . On the other hand, since the
superconductivity coexists with CDW, their interplay should be also taken into
account. The investigation of these two competitive orders is one of research
focuses in CsV3Sb5 superconductor kang2022microscopic ; STMCDW ; NMR ;
li2022conjoined ; li2021CDW ; PhysRevX.11.041030 ; PhysRevLett.128.036402 .
For instance, the observation of chiral CDW provides a possible mechanism for
unconventional superconductivity chiral , which induces a pair density wave
phase below the superconducting transition temperature PDW ; ge2022discovery .
Nevertheless, by applying high pressure, an unusual observation of two
superconducting domes was found inside and outside the CDW phase
PhysRevLett.126.247001 ; PhysRevB.103.224513 ; Pressureyu ; Chen_2021 ; wangqi
, indicating the competition between CDW and superconductivity. In addition,
CsV3Sb5 is also a $\bf{Z}$2 topological metal in which the Dirac line locates
between $\bf{K}$ and $\bf{H}$ points, and the topological surface states were
confirmed by the angle-resolved photon electron spectroscopy (ARPES)
PhysRevLett.125.247002 ; HU2021 . Moreover, evidence of anomalous Hall effect
has been observed in the CDW state, which was attributed to the symmetry
breaking of the band structures wang2021electronic ; PhysRevB.104.L041103 ;
fu2021quantum . To address these issues, one of the most promising paths is
investigating the possible chirality of Cooper pairs and vortex dynamics,
which has not yet been reported.

Here, we report the study of the nonreciprocal charge transport of
mechanically exfoliated CsV3Sb5 thin flakes. As reducing the thickness, the
superconducting transition temperature ($T_{\textrm{c}}$) is firstly enhanced
but then suppressed, competing with the emergence of CDW. The nonreciprocal
charge transport is observed when the magnetic field is aligned both within
the $ab$-plane and along the $c$-axis. The second harmonic magnetoresistance
splits into several peaks and reverses their signs by ramping the magnetic
field or the current, suggesting a tunable rectification effect. Our results
reveal detailed information about the chirality and vortex dynamics.

## Results

### Sample characterization

CsV3Sb5 is a layered material with a V-kagome lattice located in the middle of
the unit cell shown in Fig. 1a. Thanks to the weakly bonding strength between
Cs and Sb atoms, atomically thin crystals can be obtained by the mechanical
exfoliation method. Thin flakes were covered by few-layer BN, and then
transferred onto the pre-patterned gold leads as shown in Fig. 1b. Before the
nonreciprocity measurements, we first examined the onset and zero resistance
$T_{\textrm{c}}$s, and their thickness dependence is shown in Fig. 1c. The
onset $T_{\textrm{c}}$ of bulk crystal is 3.5 K, and the $T_{\textrm{c}}$ can
be enhanced up to 4.25 K by shrinking the thickness of the crystal into 80 nm,
followed by a decrease in $T_{\textrm{c}}$ to 2.5 K for the 10-nm-thick one.
The $T_{\textrm{c}}$ of zero-resistance behaves the similarly behavior. The
nonmonotonically evolution of $T_{\textrm{c}}$ with reducing layers may be
correlated to the competition with CDW order. Particularly, in the 80-nm-thick
sample, the transition temperature of CDW is reduced to 77 K, which can be
confirmed by both Hall resistances and d$R$/d$T$ as shown in Fig. 1d. This
phenomenon has been well discussed in recent reports song2021competing ;
song2021competition , in which the surface oxidation that occurred during the
sample fabrication process was atrributed to additional carrier doping and
$T_{\textrm{c}}$ modification song2021competition . In our present work,
however, a piece of BN was covered onto the thin flakes to avoid oxidation,
the $T_{\textrm{c}}$ is generally related to the thickness, we can hardly
ascribe to the sample degradation. Here, we focus on the samples with
$T_{\textrm{c}}$ of above 4 K (a 45-nm-thick device in the manuscript and
others in Supplementary Figure 1 and 2) to study the vortex dynamics and the
chirality of superconductivity from the nonreciprocal charge transport by
probing the second harmonic response.

### Nonreciprocal magnetoresistance

The first and second harmonic voltages ($V^{\omega}$ and $V^{2\omega}$) were
detected using lock-in amplifier. According to the definition of
nonreciprocity, the total voltage $V$ can be considered as
$V=V^{\omega}+V^{2\omega}=V^{\omega}\left(1+\gamma BI\right)$, where, the
coefficient $\gamma$ is equal to $2V^{2\omega}/V^{\omega}BI$. Since the
electronic structure of CsV3Sb5 is quasi-two-dimensional, the nonreciprocal
response could be anisotropic from in-plane to out-of-plane. We then measured
and calculated the temperature dependence of the first and the second harmonic
resistances ($R_{\textrm{xx}}^{\omega}=V^{\omega}/I_{\textrm{0}}$ and
$R_{\textrm{xx}}^{2\omega}=V^{2\omega}/I_{\textrm{0}}$) simultaneously for
both $B/\mkern-4.0mu/ab$-plane and $B/\mkern-4.0mu/c$-axis cases as shown in
Fig. 2a and 2b. At low magnetic field limit for $B/\mkern-4.0mu/ab$-plane
case, the nonreciprocal resistance immediately develops when $R_{\textrm{xx}}$
starts to decrease as marked in the green dashed line. On the other hand, the
$R_{\textrm{xx}}^{2\omega}$ disappears at high magnetic field limit. Similar
behavior is observed for $B/\mkern-4.0mu/c$-axis case as well. Another
important feature is that several times of sign reversals appear at both
$R_{\textrm{xx}}^{2\omega}$-$T$ and $R_{\textrm{xx}}^{2\omega}$-$B$ curves,
which will be discussed in detail below. It should be noted that the
nonreciprocal magnetoresistance is detected just within the superconductivity
transition region. Thus, we focus on the investigation between nonreciprocal
magnetoresistance and superconductivity.

Generally, the superconducting transition of the type-II superconductor
consists of two regions, separated by a mean-field transition temperature
$T_{\textrm{c}}^{\textrm{m}}$ Itahashieaay9120 ; MCAYasuda . The
$T_{\textrm{c}}^{\textrm{m}}$ can be determined as shown in Fig. 2c by using
the Aslamazov-Larkin term
$\sigma_{\textrm{2D}}=A(T-T_{\textrm{c}}^{\textrm{m}})^{-1}$, where
$\sigma_{\textrm{2D}}$ is the two-dimensional conductivity and $A$ is a
fitting parameter. Thus, we can confirm the $T_{\textrm{c}}^{\textrm{m}}$ is
about 4.00 K. Above the $T_{\textrm{c}}^{\textrm{m}}$, the transport
properties are dominated by the amplitude fluctuation of order parameters
which causes a paraconductivity state. Below the
$T_{\textrm{c}}^{\textrm{m}}$, a phase fluctuation generates vortices motion,
which takes charge of the transport properties. It is worth noting that the
resistance is still nonzero until about 3.45 K due to the current-driven
vortex motion. Below 3.45 K, the vortices are frozen as the vortex solid state
(or vortex lattice for the type-II superconductors) and the system is in
superconductivity state completely. Meanwhile, the temperature-dependent
$R_{\textrm{xx}}^{2\omega}$ curve demonstrates four characteristic regions as
well, namely, normal state, paraconductivity, vortex motion, and vortex solid
state according to those of the $R_{\textrm{xx}}^{\omega}-T$ curve. For the
magnetoresistance measurements, the $R_{\textrm{xx}}^{\omega}$ reveals a
typical magnetic-field-dependent behavior as shown in Fig. 2d, while the
$R_{\textrm{xx}}^{2\omega}$ curve demonstrates antisymmetric to the magnetic
field, namely, a typical characteristic of nonreciprocity.

To explore the magnetochiral properties correlated with the superconductivity,
we investigated the field- and temperature-dependent
$R_{\textrm{xx}}^{2\omega}$. Figure 3a and 3b show the
$R_{\textrm{xx}}^{2\omega}-B$ curves in different temperatures under magnetic
fields in the ab-plane and along the c-axis, respectively. The phase diagrams
of $R_{\textrm{xx}}^{2\omega}$ with respect to the magnetic field and the
temperature are shown in Fig. 3c and 3d, respectively. The dashed lines are
the contours of the first harmonic resistance $R_{\textrm{xx}}^{\omega}$ at
zero resistance 0$R_{\textrm{N}}$, half of the normal state resistance
(0.5$R_{\textrm{N}}$), three-quarters of the normal state
resistance(0.75$R_{\textrm{N}}$), and the onset of transition. The
nonreciprocal magnetotransport mainly locates at the region between the
contours of 0$R_{\textrm{N}}$ and 0.5$R_{\textrm{N}}$ when the field is in the
$ab$-plane and along the $c$-axis. Therefore, the nonreciprocal response
basically occurs at the vortex flow regions, indicating that the vortex motion
could be the domination. Above the contour of $R=0.5R_{\textrm{N}}$,
$R_{\textrm{xx}}^{2\omega}$ decays gradually because the magnetic field
diminishes the pinning potentials and the vortices move more freely. Above the
contour of $R=0.75R_{\textrm{N}}$, the free flow regions continuously evolved
into the amplitude fluctuation region and nonreciprocity is almost
indiscernible as demonstrated in Fig. 3a and 3b. However, the
$R_{\textrm{xx}}^{2\omega}$ did not disappear until $R$ reaches
0.95$R_{\textrm{N}}$ (Fig. 2d) at 3.8 K, where the contribution of amplitude
fluctuation of order parameter may not be neglected.

On the other hand, the $R_{\textrm{xx}}^{2\omega}$ under in-plane field is
larger than that along the $c$-axis, suggesting an anisotropic nonreciprocal
magnetotransport. Since the crystal is in a quasi-two-dimensional structure,
anisotropic superconducting properties should also affect the vortex motion.
We then examined the coherence lengths to estimate the size of the vortex from
the anisotropic upper critical field $B_{\textrm{c}}$ (at
$R=0.99R_{\textrm{N}}$). According to the Ginzburg-Landau theory
tinkham2004introduction , the in-plane and out-of-plane coherence lengths are
$\xi_{\text{ab}}^{2}=\frac{\Phi_{\textrm{0}}}{2\pi\mu_{\textrm{0}}B_{\text{c}}^{2}}$
and
$\xi_{\textrm{c}}^{2}=\frac{\Phi_{\textrm{0}}}{2\pi\mu_{\textrm{0}}\sqrt{B_{\textrm{c}}^{\textrm{c}}B_{\textrm{c}}^{\textrm{ab}}}}$,
which are given in Fig.4 a, where $\Phi_{\textrm{0}}$ is flux quanta. Thus,
the anisotropic factor can be estimated as
$\lambda=\xi_{\textrm{c}}/\xi_{\textrm{ab}}$ (=1.72 at 1.8 K), indicating a
weak anisotropy. Meanwhile, the thickness of the sample is comparable to the
diameter of a vortex of $2\xi_{\textrm{ab}}$. Once $B/\mkern-4.0mu/ab$-plane,
$R_{\textrm{xx}}^{2\omega}$ is enhanced due to the size effect and/or possible
anisotropic pinning potential PhysRevB.98.054510 . Consequently, we conclude
that the asymmetric property of pinning potential is three-dimensional.

Moreover, several specific series of peaks are observed other than one main
peak predicted by the theory of vortex dynamics PhysRevB.98.054510 .
Particularly, four series can be identified when $B/\mkern-4.0mu/ab$-plane and
much more series exist in $B/\mkern-4.0mu/c$-axis case. The locations of each
series of peaks are weakly temperature-dependent. In other words, the
nonreciprocity is strong when some number of vortices exist in the device,
indicating that the many-body interaction of vortices may also play a role
de2006controlled . We select three major series of peaks in Fig. 3a to
estimate the $\gamma$ value as $\gamma=\frac{2V^{2\omega}}{V^{\omega}BI}$, and
plot temperature dependence of $\gamma$ in Fig. 4b. While the magnitude of
these three peaks is comparable in the order of several m$\Omega$ at 2.8 K,
the coefficient $\gamma$ varied from about 5 T-1A-1 to near 8000 T-1A-1 for
high-$H$ peaks and low-$H$ peaks, respectively. It is reasonable that the
magnetic field weakens the pinning potential. The temperature plays a similar
role as all these coefficients drastically decrease with increasing
temperature. To compare with other systems, the values of
$\frac{2V^{2\omega}}{V^{\omega}j}$, where $j$ is the current density, are also
estimated as up to about $10^{-11}$m2A-1. We then conclude that the ratchet
effect is remarkable in CsV3Sb5, comparable to the noncentrosymmetric 2D
materials such as MoS2 and NbSe2.

### Current dependent nonreciprocity

We then focus on the current dependence of the nonreciprocity, which should
follow the relation of $R_{\textrm{xx}}^{2\omega}\sim\gamma BI$. The phase
diagrams of the first and the second harmonic signals are plotted in Fig. 5a
and b, respectively. The effect of current on $R_{\textrm{xx}}^{2\omega}$ is
similar to the effect of temperature as it suppresses the superconductivity
and reduces the region of $R_{\textrm{xx}}^{2\omega}$. Meanwhile, several
series of peaks and the corresponding sign reversals are observed as well. The
maximum of the second harmonic resistance is extracted as a function of
current as given in Fig. 5c, where the linear relation at current lower than
0.2 mA is verified for both $B/\mkern-4.0mu/ab$-plane and
$B/\mkern-4.0mu/c$-axis cases. When applying a larger current, the signal
starts to reduce because of the suppression of the rectification effect of
vortex motion by the relative weakening of the pinning potentials MCAZhang ;
Itahashieaay9120 ; PhysRevB.72.064522 ; PhysRevB.42.2639 ;
PhysRevMaterials.4.074003 . The current dependence of coefficient $\gamma$ of
these maximum peaks behaves differently as shown in Fig. 5d, which drops
exponentially as increasing the current when $B/\mkern-4.0mu/$ab-plane, while
no such tendency can be observed clearly for $B/\mkern-4.0mu/$c-axis, possibly
owing to the maximum peaks belonging to different series of peaks. It is worth
noting that the nonreciprocity can be modified by current as well as $B$ and
$T$, providing promising applications such as superconducting diodes MCAAndo .

### Discussion

For the mechanism of the nonreciprocity phenomenon in superconductors, several
explanations have been proposed for superconducting SrTiO3Itahashieaay9120 ,
PbTaSe2ideue2020giant , Bi2Te3/FeTe interfaceMCAYasuda , and NbN thin
filmsnakamura2020 , which can be basically attributed to the intrinsic
symmetry breaking through the vortex dynamics such as viscous vortex flow,
ratchet effect, or specific mechanisms of superconductivity under spontaneous
symmetry-breaking conditions. Under the vortex motion regime, one possibility
is that the vortex liquid state and thermally activated vortex state may
contribute to the response with opposite signs MCAZhang . The viscous vortex
flow is also possible to the nonreciprocal signal, which predicts a monotonic
enhancement of $\gamma$ as decreasing temperature. Nevertheless, one can only
observe an individual nonreciprocity peak for each mechanism instead of multi-
peaks in our present results. Another mechanism is the effect of periodical
pinning potential when the vortices match the lattice constant or artificially
structured pinning sites, resulting in the modulation of condensing
energyshiomi2017oscillatory ; nagao2006periodic ; swiecicki2012strong ;
lyu2021superconducting . After breaking the inversion symmetry of the pinning
potential, a well-controlled many-body interaction can realize sign reversal
in artificial structures de2006controlled . For single crystals, however, the
symmetry of pinning potential should reflect the intrinsic symmetry of the
crystal as well studied on MoS2 MoS2PRR  and PbTaSe2 ideue2020giant , where
the anisotropy of nonreciprocity follows that of crystal. Therefore, the
appearance of series of nonreciprocity peaks in our present results suggests
an intrinsic symmetry breaking as well.

One possible spontaneous symmetry breaking can be from the electronic
structure, which is mainly dominated by the nature of the crystal structure
and quantum phases. Hoshino $et$ $al.$ proposed a possible mechanism
PhysRevB.98.054510  in noncentrosymmetric crystal structures and interfaces
such as Rashba superconductor, and topological surface state with
superconducting proximity effect. For the intrinsic asymmetry like spin-
orbital effects and time-reversal symmetry breaking, an asymmetric energy
dispersion relation can be induced as $E\left(s,\textrm{{k}}\right)\neq
E\left(s,-\textrm{{k}}\right)$ and $E\left(s,\textrm{{k}}\right)\neq
E\left(-s,-\textrm{{k}}\right)$, which is the key to the nonreciprocity. In
bulk CsV3Sb5, the Rashba spin-orbital interaction seems unlikely to exist,
because the crystal structure is centrosymmetric. Alternatively, a recent
study of Raman spectrum in CsV3Sb5 thin films revealed an intrinsic
instability of crystalline structure liu2022anomalous , and STM measurement
also shows an inhomogeneous Cs-surface RN987 . Such degradation may reduce the
symmetry and modify its electronic structure zhang2022emergence , providing a
possibility to induce nonreciprocity. Nevertheless, these effects ought to be
weak since our crystals is relatively thick, and particularly, the samples are
single-crystalline and well protected by BN, the degradation of crystalline
can be basically avoided.

In addition, the charge orders can also modulate the rotational, translation,
and time reversal symmetry. In CsV3Sb5, the nematicity RN987  reveals a
C${}_{\textrm{2}}$ rotational symmetry, and the unidirectional
4$\bf{a}_{\textrm{0}}$ CDW also breaks rotational symmetry li2022rotation .
The amplitudes of three q vectors within the in-plane
2$a_{\textrm{0}}\times$2$a_{\textrm{0}}$ CDW are chiral, which can break the
time-reversal symmetry as well chiral ; wang2021electronic . However, the CDW-
induced lattice distortion CsV3Sb5 was found to be centrosymmetric
PhysRevX.11.041030 , and no any evidence suggests that one of these charge
orders can break the inversion symmetry up to now.

Another possible origin of the spontaneously symmetric breaking is the
unconventional superconducting state. Among a variety of possible pairing
states, superconducting states on topological bands should keep a nontrivial
nature and the pair potential demonstrates anisotropic, for instance a
spinless $p$-wave superconductor mackenzie . Previously, the nonreciprocal
transport has been comprehensively studied on the interfacial systems of
Bi2Te3/FeTe MCAYasuda , which were attributed to inversion symmetry breaking
of the topological superconductivity in the interface. The effect of other
unconventional states, for example, chiral superconductors zinkl2021symmetry ,
still needs further theoretical analysis and experimental evidence. Although
the superconductivity symmetry of $A$V3Sb5 family is still an open question as
discussed in the introduction part, the ARPES and first principles
calculations PhysRevLett.125.247002 ; HU2021 ; RN986 ; PhysRevLett.127.177001
indicate the existence of topological bands and surface states, and a recent
work reported an observation of two-fold symmetry on the $c$-axis resistivity
under in-plane rotation of magnetic field RN985 . Therefore, our nonreciprocal
transport phenomenon should be a promising probe to reflect the anisotropic
nature of the unconventional superconductivity state in CsV3Sb5.

In conclusion, we have carefully investigated the nonreciprocity under the
dissipative state of CsV3Sb5 superconductors. With the reduction of thickness,
we found that the superconductivity is firstly enhanced and then suppressed
due to the competition with CDW. Strong nonreciprocal signals can be detected
when superconductivity occurs, whose strength is comparable to the
artificially structured or noncentrosymmetric superconductors. The second
harmonic resistance can be observed as the magnetic field both in the
$ab$-plane and along the $c$-axis, which can be basically attributed to the
vortex ratchet motion. The nonreciprocal signal with respect to magnetic
fields splits into several series of peaks that are asymmetric with magnetic
fields. We suggest the intrinsic symmetry breaking in superconducting states
may induce the magnetochirality in the CsV3Sb5 superconductor.

## Methods

### Crystal growth

High-quality single crystal were synthesized by Sb flux method and described
in previous workwangqi .

### Device fabrication

The thin crystals of CsV3Sb5 were prepared by mechanical exfoliation from bulk
crystal and transferred on pre-patterned SiO2(300 nm)/Si substrate using
polydimethylsiloxane (PDMS) stamps in a laboratory-made vacuum transfer
system. The pre-patterned leads were fabricated by standard photolithography
method using Photoresist (AZ5214) and subsequent deposition of Ti/Au (5 nm/20
nm) by magnetron sputtering. Some of the samples were then covered by a $h$-BN
capping layer for further degradation. The thicknesses of the samples were
confirmed by AFM after transport measurement.

### Transport measurements

The transport measurement was carried out in PPMS. The four-terminal DC and AC
signal was measured by set of Keithley 2400 and 2182a and set of Keithley 6221
and OE1022 lock-in amplifier, respectively. The first and second harmonic
resistance were defined as $R^{\omega}=V^{\omega}/I_{\textrm{0}}$ and
$R^{2\omega}=V^{2\omega}/I_{\textrm{0}}$, where $I_{\textrm{0}}$ is the
amplitude of the AC current applied and $V^{\omega}$ and $V^{2\omega}$ are the
amplitude of first and second harmonic voltage. The current frequency was set
to be 113 Hz to lower the noise and the phase of second harmonic signal was
set as $\pi/2$.

## Data Availability Statement

The data supporting the findings of this study are available from the
corresponding authors upon reasonable request.

## Acknowledgement

This research was supported in part by the National Natural Science Foundation
of China (Grants No. 12004251, 12104302, 12104303, U1932217, and 11974246),
the Natural Science Foundation of Shanghai (Grant No. 20ZR1436100,
19ZR1477300), the Science and Technology Commission of Shanghai Municipality,
the Shanghai Sailing Program (Grant No. 21YF1429200), the start-up funding
from ShanghaiTech University, Beijing National Laboratory for Condensed Matter
Physics, the Interdisciplinary Program of Wuhan National High Magnetic Field
Center (WHMFC202124), the Elemental Strategy Initiative conducted by the MEXT,
Japan, Grant Number JPMXP0112101001, JSPS KAKENHI Grant Number JP20H00354 and
A3 Foresight by JSPS, and National Key R&D Program of China (Grant No.
2018YFA0704300).

## Competing Interests

We declare no competing interest.

## Author Contributions

Y.W. and Q.W. contributed equally to this work. Y.W. and J.L. designed the
research. Q.W. and Y.Q. grew the single crystals of CsV3Sb5. K.W. and T.T.
grew the single crystals of $h$-BN. Y.W. fabricated the samples and performed
the electrical transport measurements. Y.Z. performed the atomic force
microscopy experiment. Y.Z., P.D., J.H., Y.D., Y.L., B.T., C.Z., and H.Z.
assisted the device fabrication. Y.W., X.Z., J.W., Q.Y., and J.L. contributed
to analyze the data. Y.W. and J.L. wrote the paper.

## References

  *   * (1) Rikken, G. & Raupach, E.  Observation of magneto-chiral dichroism.  _Nature_ 390, 493–494 (1997) .
  * (2) Zhou, X., Zeng, F., Jia, M., Chen, H. & Wu, Y.  Sign reversal of unidirectional magnetoresistance in monocrystalline Fe/Pt bilayers.  _Phys. Rev. B_ 104, 184413 (2021) .
  * (3) Tokura, Y. & Nagaosa, N.  Nonreciprocal responses from non-centrosymmetric quantum materials.  _Nat. Commun._ 9, 3740 (2018) .
  * (4) Villegas, J. E. _et al._ A Superconducting Reversible Rectifier That Controls the Motion of Magnetic Flux Quanta.  _Science_ 302, 1188–1191 (2003) .
  * (5) Moshchalkov, V., Woerdenweber, R. & Lang, W.  _Nanoscience and engineering in superconductivity_ (Springer Science & Business Media, 2010).
  * (6) Zhang, E. _et al._ Nonreciprocal superconducting NbSe2 antenna.  _Nat. Commun._ 11, 5634 (2020) .
  * (7) Nagulu, A., Reiskarimian, N. & Krishnaswamy, H.  Non-reciprocal electronics based on temporal modulation.  _Nat. Electronic._ 3, 241–250 (2020) .
  * (8) Liu, Z. _et al._ Observation of nonreciprocal magnetophonon effect in nonencapsulated few-layered CrI3.  _Sci. Adv._ 6, eabc7628 (2020) .
  * (9) Hoshino, S., Wakatsuki, R., Hamamoto, K. & Nagaosa, N.  Nonreciprocal charge transportin two-dimensional noncentrosymmetric superconductors.  _Phys. Rev. B_ 98, 054510 (2018) .
  * (10) Daido, A., Ikeda, Y. & Yanase, Y.  Intrinsic superconducting diode effect.  _Phys. Rev. Lett._ 128, 037001 (2022) .
  * (11) Itahashi, Y. M., Saito, Y., Ideue, T., Nojima, T. & Iwasa, Y.  Quantum and classical ratchet motions of vortices in a two-dimensional trigonal superconductor.  _Phys. Rev. Res._ 2, 023127 (2020) .
  * (12) Ando, F. _et al._ Observation of superconducting diode effect.  _Nature_ 584, 373–376 (2020) .
  * (13) Itahashi, Y. M. _et al._ Nonreciprocal transport in gate-induced polar superconductor SrTiO3.  _Sci. Adv._ 6, eaay9120 (2020) .
  * (14) Lyu, Y. _et al._ Superconducting diode effect via conformal-mapped nanoholes.  _Nat. Commun._ 12, 2703 (2021) .
  * (15) Yasuda, K. _et al._ Nonreciprocal charge transport at topological insulator/superconductor interface.  _Nat. Commun._ 10, 2734 (2019) .
  * (16) Ortiz, B. R. _et al._ $\mathrm{Cs}{\mathrm{V}}_{3}{\mathrm{Sb}}_{5}$: A ${\mathbb{Z}}_{2}$ Topological Kagome Metal with a Superconducting Ground State.  _Phys. Rev. Lett._ 125, 247002 (2020) .
  * (17) Wu, X. _et al._ Nature of Unconventional Pairing in the Kagome Superconductors A${\mathrm{V}}_{3}{\mathrm{Sb}}_{5}$ (A$=\mathrm{K},\mathrm{Rb},\mathrm{Cs}$).  _Phys. Rev. Lett._ 127, 177001 (2021) .
  * (18) Ding, P., Lee, C. H., Wu, X. & Thomale, R.  Diagnosis of pairing symmetry by vortex and edge spectra in kagome superconductors.  _Phys. Rev. B_ 105, 174518 (2022) .
  * (19) Kim, S.-W., Oh, H., Moon, E.-G. & Kim, Y.  Monolayer Kagome Metals AV3Sb5.  Preprint at 10.48550/arXiv.2202.11521 (2022) .
  * (20) Jiang, H.-M., Yu, S.-L. & Pan, X.-Y.  Electronic structure and spin-lattice relaxation in superconducting vortex states on the kagome lattice near van hove filling.  _Phys. Rev. B_ 106, 014501 (2022) .
  * (21) Mu, C. _et al._ S-Wave Superconductivity in Kagome Metal CsV3Sb5 Revealed by 121/123Sb NQR and 51V NMR Measurements.  _Chin. Phys. Lett._ 38, 077402 (2021) .
  * (22) Li, S. _et al._ Nodal superconductivity and superconducting domes in the topological Kagome metal CsV3Sb5.  Preprint at 10.48550/arXiv.2102.08356 (2021) .
  * (23) Xu, H.-S. _et al._ Multiband Superconductivity with Sign-Preserving Order Parameter in Kagome Superconductor CsV3Sb5.  _Phys. Rev. Lett._ 127, 187004 (2021) .
  * (24) Duan, W. _et al._ Nodeless superconductivity in the kagome metal CsV3Sb5.  _Sci. China Phys. Mech._ 64, 107462 (2021) .
  * (25) Ni, S. _et al._ Anisotropic Superconducting Properties of Kagome Metal CsV3Sb5.  _Chin. Phys. Lett._ 38, 057403 (2021) .
  * (26) Xiang, Y. _et al._ Twofold symmetry of c-axis resistivity in topological kagome superconductor CsV3Sb5 with in-plane rotating magnetic field.  _Nat. Commun._ 12, 6727 (2021) .
  * (27) Kang, M. _et al._ Twofold van Hove singularity and origin of charge order in topological kagome superconductor CsV3Sb5.  _Nat. Phys._ 18, 301–308 (2022) .
  * (28) Jiang, Y.-X. _et al._ Unconventional chiral charge order in kagome superconductor KV3Sb5.  _Nat. Mater._ 20, 1353 (2021) .
  * (29) Tan, H., Liu, Y., Wang, Z. & Yan, B.  Charge Density Waves and Electronic Properties of Superconducting Kagome Metals.  _Phys. Rev. Lett._ 127, 046401 (2021) .
  * (30) Feng, X., Jiang, K., Wang, Z. & Hu, J.  Chiral flux phase in the Kagome superconductor AV3Sb5.  _Sci. Bull._ 66, 1384–1388 (2021) .
  * (31) Yu, L. _et al._ Evidence of a hidden flux phase in the topological kagome metal CsV3Sb5.  Preprint at 10.48550/arXiv.2107.10714 (2021) .
  * (32) Kang, M. _et al._ Microscopic structure of three-dimensional charge order in kagome superconductor AV3Sb5 and its tunability.  Preprint at 10.48550/arXiv.2202.01902 (2022) .
  * (33) Zhao, H. _et al._ Cascade of correlated electron states in the kagome superconductor CsV3Sb5.  _Nature_ 599, 216–221 (2021) .
  * (34) Song, D. _et al._ Orbital ordering and fluctuations in a kagome superconductor CsV3Sb5.  _Sci. China Phys. Mech._ 65, 247462 (2022) .
  * (35) Li, H. _et al._ Conjoined Charge Density Waves in the Kagome Superconductor CsV3Sb5.  Preprint at 10.48550/arXiv.2202.13530 (2022) .
  * (36) Li, H. _et al._ Observation of Unconventional Charge Density Wave without Acoustic Phonon Anomaly in Kagome Superconductors AV3Sb5 (A= Rb, Cs).  _Phys. Rev. X_ 11, 031050 (2021) .
  * (37) Ortiz, B. R. _et al._ Fermi Surface Mapping and the Nature of Charge-Density-Wave Order in the Kagome Superconductor CsV3Sb5.  _Phys. Rev. X_ 11, 041030 (2021) .
  * (38) Lou, R. _et al._ Charge-Density-Wave-Induced Peak-Dip-Hump Structure and the Multiband Superconductivity in a Kagome Superconductor CsV3Sb5.  _Phys. Rev. Lett._ 128, 036402 (2022) .
  * (39) Chen, H. _et al._ Roton pair density wave in a strong-coupling kagome superconductor.  _Nature_ 599, 222–228 (2021) .
  * (40) Ge, J. _et al._ Discovery of charge-4e and charge-6e superconductivity in kagome superconductor CsV3Sb5.  Preprint at 10.48550/arXiv.2201.10352 (2022) .
  * (41) Chen, K. Y. _et al._ Double Superconducting Dome and Triple Enhancement of ${T}_{c}$ in the Kagome Superconductor CsV3Sb5 under High Pressure.  _Phys. Rev. Lett._ 126, 247001 (2021) .
  * (42) Zhang, Z. _et al._ Pressure-induced reemergence of superconductivity in the topological kagome metal CsV3Sb5.  _Phys. Rev. B_ 103, 224513 (2021) .
  * (43) Yu, F. H. _et al._ Unusual competition of superconductivity and charge-density-wave state in a compressed topological kagome metal.  _Nat. Commun._ 12, 3645 (2021) .
  * (44) Chen, X. _et al._ Highly Robust Reentrant Superconductivity in CsV3Sb5 under Pressure.  _Chin. Phys. Lett._ 38, 057402 (2021) .
  * (45) Wang, Q. _et al._ Charge Density Wave Orders and Enhanced Superconductivity under Pressure in the Kagome Metal CsV3Sb5.  _Adv. Mater._ 33, 2102813 (2021) .
  * (46) Hu, Y. _et al._ Topological surface states and flat bands in the kagome superconductor CsV3Sb5.  _Sci. Bull._ 67 (5), 495–500 (2022) .
  * (47) Wang, Z. _et al._ Electronic nature of chiral charge order in kagome superconductor CsV3Sb5.  _Phys. Rev. B_ 104, 075148 (2021) .
  * (48) Yu, F. H. _et al._ Concurrence of anomalous Hall effect and charge density wave in a superconducting topological kagome metal.  _Phys. Rev. B_ 104, L041103 (2021) .
  * (49) Fu, Y. _et al._ Quantum transport evidence of topological band structures of kagome superconductor CsV3Sb5.  _Phys. Rev. Lett._ 127, 207002 (2021) .
  * (50) Song, B. Q. _et al._ Competing superconductivity and charge-density wave in Kagome metal CsV3Sb5: evidence from their evolutions with sample thickness.  Preprint at 10.48550/arXiv.2105.09248 (2021) .
  * (51) Song, Y. _et al._ Competition of superconductivity and charge density wave in selective oxidized CsV3Sb5 thin flakes.  _Phys. Rev. Lett._ 127, 237001 (2021) .
  * (52) Tinkham, M.  _Introduction to superconductivity_ (Courier Corporation, 2004).
  * (53) de Souza Silva, C. C., Van de Vondel, J., Morelle, M. & Moshchalkov, V. V.  Controlled multiple reversals of a ratchet effect.  _Nature_ 440, 651–654 (2006) .
  * (54) Yan, H., Abdelhadi, M. M., Jung, J. A., Willemsen, B. A. & Kihlstrom, K. E.  Exponential dependence of the vortex pinning potential on current density in high-${T}_{c}$ superconductors.  _Phys. Rev. B_ 72, 064522 (2005) .
  * (55) Maley, M. P., Willis, J. O., Lessure, H. & McHenry, M. E.  Dependence of flux-creep activation energy upon current density in grain-aligned ${\mathrm{YBa}}_{2}$${\mathrm{Cu}}_{3}$${\mathrm{O}}_{7\mathrm{-}\mathit{x}}$.  _Phys. Rev. B_ 42, 2639 (1990) .
  * (56) Saito, Y., Itahashi, Y. M., Nojima, T. & Iwasa, Y.  Dynamical vortex phase diagram of two-dimensional superconductivity in gated $\mathrm{MoS}_{2}$.  _Phys. Rev. Mater._ 4, 074003 (2020) .
  * (57) Ideue, T., Koshikawa, S., Namiki, H., Sasagawa, T. & Iwasa, Y.  Giant nonreciprocal magnetotransport in bulk trigonal superconductor PbTaSe2.  _Phys. Rev. Res._ 2, 042046 (2020) .
  * (58) Nakamura, S., Katsumi, K., Terai, H. & Shimano, R.  Nonreciprocal terahertz second-harmonic generation in superconducting nbn under supercurrent injection.  _Phys. Rev. Lett._ 125, 097004 (2020) .
  * (59) Shiomi, Y., Lustikova, J. & Saitoh, E.  Oscillatory Nernst effect in Pt ferrite cuprate-superconductor trilayer films.  _Sci. Rep._ 7, 1–8 (2017) .
  * (60) Nagao, M. _et al._ Periodic oscillations of Josephson-vortex flow resistance in oxygen-deficient YBa2Cu3Ox.  _Phys. Rev. B_ 74, 054502 (2006) .
  * (61) Swiecicki, I. _et al._ Strong field-matching effects in superconducting YBa2Cu3O7-δ films with vortex energy landscapes engineered via masked ion irradiation.  _Phys. Rev. B_ 85, 224502 (2012) .
  * (62) Liu, G. _et al._ Anomalous amplitude modes and thinning-induced loss of crystallinity in the kagome metal CsV3Sb5.  Preprint at 10.48550/arXiv.2201.05330 (2022) .
  * (63) Nie, L. _et al._ Charge-density-wave-driven electronic nematicity in a kagome superconductor.  _Nature_ 604, 59–64 (2022) .
  * (64) Zhang, W. _et al._ Emergence of large quantum oscillation frequencies in thin flakes of a kagome superconductor CsV3Sb5.  Preprint at 10.48550/arXiv.2202.08570 (2022) .
  * (65) Li, H. _et al._ Rotation symmetry breaking in the normal state of a kagome superconductor KV3Sb5.  _Nat. Phys._ 18, 265–270 (2022) .
  * (66) Mackenzie, A. P. & Maeno, Y.  The superconductivity of Sr2RuO4 and the physics of spin-triplet pairing.  _Rev. Mod. Phys._ 75, 657–712 (2003) .
  * (67) Zinkl, B., Hamamoto, K. & Sigrist, M.  Symmetry conditions for the superconducting diode effect in chiral superconductors.  _Phys. Rev. Res._ 4, 033167 (2022) .

## Figure Legends

Figure 1:  Characterization of superconductivity.  a Crystal structure of
CsV3Sb5. b Optical image of a sample covered with $h$-BN. Scale bar: 20
$\mathrm{\mu}$m. c Normalized resistance as a function of temperature with
different thickness. d Temperature dependence of Hall resistance (red) and
differential resistance (blue) of a 80-nm-thick sample under $H=$ 2 T. Here
the resistance is measured with direct current.  Figure 2:  Anisotropic
nonreciprocal magnetoresistance of a 45-nm-thick sample. (a) (b) $R-T$ curves
of first (top) and second (bottom) harmonic signal under $B$ applied along the
$ab$-plane and the $c$-axis, respectively. Scale bar: 2 m$\Omega$. The dashed
line is the guideline of onset temperature of second harmonic signal. (c) The
transition regions for $R_{\textrm{xx}}^{2\omega}-T$ curve under $B=0.1$ T.
Red dashed line is the Aslamazov-Larkin fitting curve of paraconductivity
contribution. (d) The $R_{\textrm{xx}}^{2\omega}-B$ and
$R_{\textrm{xx}}^{\omega}-B$ curves at 3.8 K, where $B$ is applied within the
$ab$-plane.  Figure 3:  The magnetic field and temperature dependence of
$R_{\textrm{xx}}^{2\omega}$ of the 45-nm-thick sample. a and b are the second
harmonic magnetoresistance with the magnetic field applied along the
$ab$-plane and the $c$-axis at different temperatures. Scale bar: 2.5
m$\Omega$. c and d are the color mapping plot of a and b, respectively. The
dashed lines are the contours of first harmonic resistance $R^{\omega}$ of
zero resistance, half of the normal state resistance (0.5$R_{\textrm{N}}$),
0.75$R_{\textrm{N}}$ and onset of transition.  Figure 4:  Coherence lengths
and nonlinear coefficient. a Critical field and coherence length are extracted
for $B/\mkern-4.0mu/ab$-plane and $B/\mkern-4.0mu/c$-axis cases. Black dotted
line shows the thickness of the sample. b The nonlinear coefficient
$\gamma=\frac{2V^{2\omega}}{V^{\omega}BI}$ is calculated from several series
of peaks.  Figure 5:  Current dependence of nonreciprocity of the 45-nm-thick
sample. a and b are color mapping of the $R_{\textrm{xx}}^{2\omega}$ as a
function of magnetic field and current. The dashed lines are the contours of
first harmonic resistance $R_{\textrm{xx}}^{\omega}$ of 0$R_{\textrm{N}}$ ,
0.5$R_{\textrm{N}}$, 0.75$R_{\textrm{N}}$ and onset of transition. c Current
dependence of the maximum second harmonic resistance for
$B/\mkern-4.0mu/ab$-plane and $B/\mkern-4.0mu/c$-axis. Error bars are standard
deviations of the resistances. d The corresponding coefficient $\gamma$ as
estimated by $\gamma=\frac{2V^{2\omega}}{V^{\omega}BI}$.


# Plyades: A Python Library for Space Mission Design

Helge Eichhorn∗†, Reiner Anderl† * Corresponding author<EMAIL_ADDRESS>darmstadt.de† Technische Universität Darmstadt, Department of Computer
Integrated DesignCopyright © 2015 Helge Eichhorn et al. This is an open-access
article distributed under the terms of the Creative Commons Attribution
License, which permits unrestricted use, distribution, and reproduction in any
medium, provided the original author and source are credited.
http://creativecommons.org/licenses/by/3.0/

###### Abstract

Designing a space mission is a computation-heavy task. Software tools that
conduct the necessary numerical simulations and optimizations are therefore
indispensable. The usability of existing software, written in Fortran and
MATLAB, suffers because of high complexity, low levels of abstraction and out-
dated programming practices. We propose Python as a viable alternative for
astrodynamics tools and demonstrate the proof-of-concept library Plyades which
combines powerful features with Pythonic ease of use.

###### Index Terms:

data modeling, object-oriented programming, orbital mechanics, astrodynamics

##  1 Introduction

Designing a space mission trajectory is a computation-heavy task. Software
tools that conduct the necessary numerical simulations and optimizations are
therefore indispensable and high numerical performance is required. No science
mission or spacecraft are exactly the same and during the early mission design
phases the technical capabilities and constraints change frequently. Therefore
high development speed and programmer productivity are required as well. Due
to its supreme numerical performance Fortran has been the top programming
language in many organizations within the astrodynamics community.

At the European Space Operations Center (ESOC) of the European Space Agency
(ESA) a large body of sometimes decades old Fortran77 code remains in use for
mission analysis tasks. While this legacy code mostly fulfills the performance
requirements1 usability and programmer productivity suffer. One reason for
this is that Fortran is a compiled language which hinders exploratory analyses
and rapid prototyping. The low level of abstraction supported by Fortran77 and
programming conventions, like the 6-character variable name limit and fixed-
form source code, that are in conflict with today’s best practices are a more
serious problem, though. The more recent Fortran standards remedy a lot of
these shortcomings, e.g. free-form source in Fortran90 or object-oriented
programming features in Fortran2003, but also introduce new complexity, e.g.
requiring sophisticated build systems for dependency resolution. Compiler
vendors have also been very slow to implement new standards. For example this
year the Intel Fortran compiler achieved full support of the Fortran2003
standard, which was released in 2005
[IFC15].\raisebox{10.00002pt}{\hypertarget{id3}{}}\hyperlink{id1}{1}\raisebox{10.00002pt}{\hypertarget{id3}{}}\hyperlink{id1}{1}footnotetext:
Many routines were not written with thread-safety and re-entrancy in mind and
can therefore not be used in parallel codes.

Due to these reasons Fortran-based tools and libraries have been generally
used together with programming environments with better usability such as
MATLAB. A common approach for developing mission design software at ESOC is
prototyping and implementing downstream processes such as visualization in
MATLAB and then later porting performance-intensive parts or the whole system
to Fortran77. The results are added complexity through the use of the MEX-
interface for integrating Fortran and MATLAB, duplicated effort for porting,
and still a low-level of abstraction because the system design is constrained
by Fortran’s limitations.

Because of the aforementioned problems some organizations explore
possibilities to replace Fortran for future development. The French space
agency CNES (Centre National D’Études Spatiales) for instance uses the Java-
based Orekit library [Ore15] for its flight dynamics systems.

In this paper we show why Python and the scientific Python ecosystem are a
viable choice for the next generation of space mission design software and
present the Plyades library. Plyades is a proof-of-concept implementation of
an object-oriented astrodynamics library in pure Python. It makes use of many
prominent scientific Python libraries such as Numpy, Scipy, Matplotlib, Bokeh,
and Astropy. In the following we discuss the design of the Plyades data model
and conclude the paper with an exemplary analysis.

##  2 Why Python?

Perez, Granger, and Hunter [PGH11] show that the scientific Python ecosystem
has reached a high level of maturity and conclude that "Python has now entered
a phase where it’s clearly a valid choice for high-level scientific code
development, and its use is rapidly growing". This assessment also holds true
for astrodynamics work as most of the required low-level mathematical and
infrastructural building blocks are already available, as shown below:

  * •

Vector algebra: NumPy [WCV11]

  * •

Visualization: Matplotlib [JDH07], Bokeh [BDT15]

  * •

Numerical integration and optimization: SciPy [JOP01]

  * •

High performance numerics: Cython [BBC11], Numba [NDT15]

  * •

Planetary ephemerides: jplephem [BRh15]

Another advantage is Python’s friendliness to beginners. In the United States
Python was the most popular language for teaching introductory computer
science (CS) courses at top-ranked CS-departments in 2014 [PGu14].
Astrodynamicists are rarely computer scientists but mostly aerospace
engineers, physicists and mathematicians. Most graduates of these disciplines
have only little programming experience. Moving to Python could therefore
lower the barrier of entry significantly.

It is also beneficial that the scientific Python ecosystem is open-source
compared to the MATLAB environment and commercial Fortran compilers which
require expensive licenses.

##  3 Design of the Plyades Object Model

The general idea behind the design of the Plyades library is the introduction
of proper abstraction. We developed a domain model based on the following
entities which are part of every analysis in mission design:

  * •

orbits (or trajectories),

  * •

spacecraft state vectors,

  * •

and celestial bodies.

###  3.1 The Body Class

The Body class is a simple helper class that holds physical constants and
other properties of celestial bodies such as planets and moons. These include

  * •

the name of the body,

  * •

the gravitational parameter $\mu$,

  * •

the mean radius $r_{m}$,

  * •

the equatorial radius $r_{e}$,

  * •

the polar radius $r_{p}$,

  * •

the $J_{2}$ coefficient of the body’s gravity potential.

###  3.2 The State Class

To define the state of a spacecraft in space-time we need the following
information

  * •

the position vector ($\vec{r}\in\mathbf{r^{3}}$),

  * •

the velocity vector ($\vec{v}\in\mathbf{r^{3}}$),

  * •

the corresponding moment in time, the so-called epoch,

  * •

the reference frame in which the vectors are defined,

  * •

the central body which also defines the origin of the coordinate system,

  * •

and additional spacecraft status parameters such as mass.

While the information could certainly be stored in a single Numpy-array an
object-oriented programming (OOP) approach offers advantages. Since all
necessary data can be encapsulated in the object most orbital characteristics
can be calculated by calling niladic or monadic instance methods. Keeping the
number of parameters within the application programming interface (API) very
small, as recommended by Robert C. Martin [RCM08], improves usability, e.g.
the user is not required to know the order of the function parameters. OOP
also offers the opportunity to integrate the State class with the Python
object model and the Jupyter notebook to provide rich human-friendly
representations.

State vectors also provide methods for backwards and forwards propagation.
Through propagation trajectories are generated, which are instances of the
Orbit class.

###  3.3 The Orbit Class

In contrast to the State class which represents a single state in space-time
the Orbit class spans a time interval and contains several spacecraft states.
It provides all necessary tools to analyze the evolution of the trajectory
over time including

  * •

quick visualizations in three-dimensional space and two-dimensional
projections,

  * •

evolution of orbital characteristics,

  * •

and determination of intermediate state vectors.

##  4 Exemplary Usage

In this example we use the Plyades library to conduct an analysis of the orbit
of the International Space Station (ISS)2. We obtain the inital state data on
August 28, 2015, 12:00h from NASA realtime trajectory data [NAS15] and use it
to instantiate a Plyades State object as shown
below.\raisebox{10.00002pt}{\hypertarget{id17}{}}\hyperlink{id15}{2}\raisebox{10.00002pt}{\hypertarget{id17}{}}\hyperlink{id15}{2}footnotetext:
A Jupyter notebook with this analysis can be obtained from Github.

iss_r = numpy.array([ -2775.03475, 4524.24941, 4207.43331, ]) *
astropy.units.kmiss_v = numpy.array([ -3.641793088, -5.665088604, 3.679500667,
]) * astropy.units.km/astropy.units.siss_t =
astropy.time.Time(’2015-08-28T12:00:00.000’)frame = ’ECI’body =
plyades.bodies.EARTHiss = plyades.State(iss_r, iss_v, iss_t, frame, body)The
position (iss_r) and velocity (iss_v) vectors use the units functionality from
the Astropy package [ASP13] while the timestamp (iss_t) is an Astropy Time
object. The constant EARTH from the plyades.bodies module is a Body object and
provides Earth’s planetary constants.

The resulting State object contains all data necessary to describe the current
orbit of the spacecraft. Calculations of orbital characteristics are therefore
implemented with the @property decorator, like shown below, and are instantly
available.

@propertydef elements(self): return kepler.elements(self.body.mu, self.r,
self.v)We compute the following orbital elements for the orbit of the ISS:

  * •

Semi-major axis: $a=6777.773$ km

  * •

Eccentricity: $e=0.00109$

  * •

Inclination: $i=51.724$ deg

  * •

Longitude of ascending node: $\Omega=82.803$ deg

  * •

Argument of periapsis: $\omega=101.293$ deg

  * •

True anomaly: $\nu=48.984$ deg

Based on the orbital elements derived quantities like the orbital period can
be determined.

In the idealized two-body problem which assumes a uniform gravity potential
the only orbital element that changes over time is the true anomaly. It is the
angle that defines the position of the spacecraft on the orbital ellipse. By
solving Kepler’s equation we can determine the true anomaly for every point in
time and derive new Cartesian state vectors [DAV13].

kepler_orbit = iss.kepler_orbit()kepler_orbit.plot3d()We now call the
kepler_orbit instance method to solve Kepler’s equation at regular intervals
until one revolution is completed. The trajectory that comprises of the
resulting state vectors is stored in the returned Orbit object. By calling
plot3d we receive a three-dimensional visualization3 of the full orbital
ellipse as shown in figure
1.\raisebox{10.00002pt}{\hypertarget{id21}{}}\hyperlink{id20}{3}\raisebox{10.00002pt}{\hypertarget{id21}{}}\hyperlink{id20}{3}footnotetext:
The visualization suffers from the fact that Matplotlib’s mplot3d toolkit does
not support proper depth buffering. Thus Matplotlib renders the complete orbit
in front of the Earth and the parts of the orbit behind the planet are not
obscured.

Figure 1: A three-dimensional visualization of the orbit based on Matplotlib.

We can achieve a similar result, apart from numerical errors, by numerically
integrating Newton’s equation:

$\vec{\ddot{r}}=-\mu\frac{\vec{r}}{|\vec{r}|^{3}}$ (1)

Plyades uses the DOP853 integrator from the scipy.integrate suite which is an
8th-order Runge-Kutta integrator with Dormand-Prince coefficients. By default
the propagator uses adaptive step-size control and a simple force model that
only considers the uniform gravity potential (see equation 1).

newton_orbit = iss.propagate( iss.period*0.8, max_step=500,
interpolate=100)newton_orbit.plot_plane(plane=’XZ’, show_steps=True)In this
example we propagate for 0.8 revolutions and constrain the step size to 500
seconds to improve accuracy. We also interpolate additional state vectors
between the integrator steps for visualization purposes.

Figure 2: Visualization of a numerically propagated orbit with intermediate
solver steps (+, blue), start point (+, red), and end point (x, red).

The trajectory plot in figure 2 also includes markers for the intermediate
integrator steps.

Since the shape of the Earth is rather an irregular ellipsoid than a sphere
Earth’s gravity potential is also not uniform. We can model the oblateness of
the Earth by including the second dynamic form factor $J_{2}$ in the equations
of motion as shown in equation 2.

$\vec{\ddot{r}}=-\mu\frac{\vec{r}}{|\vec{r}|^{3}}-\frac{3}{2}\frac{\mu
J_{2}R_{e}^{2}}{|\vec{r}|^{5}}\begin{bmatrix}x\left(1-5\frac{z^{2}}{|\vec{r}|^{2}}\right)\\\
y\left(1-5\frac{z^{2}}{|\vec{r}|^{2}}\right)\\\
z\left(3-5\frac{z^{2}}{|\vec{r}|^{2}}\right)\end{bmatrix}$ (2)

When introducing this perturbation we should expect that the properties of the
orbit will change over time. We will now analyze these effects further.

Plyades allows the substitution of force equations with a convenient
decorator-based syntax that is illustrated in the next code listing.

@iss.gravitydef newton_j2(f, t, y, params): r =
np.sqrt(np.square(y[:3]).sum()) mu = params[’body’].mu.value j2 =
params[’body’].j2 r_m = params[’body’].mean_radius.value rx, ry, rz = y[:3]
f[:3] += y[3:] pj = -3/2*mu*j2*r_m**2/r**5 f[3] += -mu*rx/r**3 +
pj*rx*(1-5*rz**2/r**2) f[4] += -mu*ry/r**3 + pj*ry*(1-5*rz**2/r**2) f[5] +=
-mu*rz/r**3 + pj*rz*(3-5*rz**2/r**2)

Figure 3: Visualization of the perturbed orbit.

After propagating over 50 revolutions the perturbation of the orbit is clearly
visible within the visualization in figure 3. A secular (non-periodical)
precession of the orbital plane is visible. Thus a change in the longitude of
the ascending node should be present.

We can plot the longitude of the ascending node by issuing the following
command:

j2_orbit.plot_element(’ascending_node’)The resulting figure 4 shows the
expected secular change of the longitude of the ascending node.

Figure 4: Secular perturbation on the longitude of the ascending node.

##  5 Future Development

As of this writing Plyades has been superseded by the Python Astrodynamics
project [PyA15]. The project aims to merge the three MIT-licensed, Python-
based astrodynamics libraries Plyades, Poliastro [JCR15] and Orbital [FML15]
and provide a comprehensive Python-based astrodynamics toolkit for productive
use.

##  6 Conclusion

In this paper we have discussed the current tools and programming environments
for space mission design. These suffer from high complexity, low levels of
abstraction, low flexibility, and out-dated programming practices. We have
then shown why the maturity and breadth of the scientific Python ecosystem as
well as the usability of the Python programming language make Python a viable
alternative for next generation astrodynamics tools. With the design and
implementation of the proof-of-concept library Plyades we demonstrated that it
is possible to create powerful yet simple to use astrodynamics tools in pure
Python by using scientific Python libraries and following modern best
practices. The Plyades work has lead to the foundation of the Python
Astrodynamics project, an inter-european collaboration, whose goal is the
development of a production-grade Python-based astrodynamics library.

## References

  * [ASP13] The Astropy Collaboration. _Astropy: A community Python package for astronomy_ , Astronomy & Astrophysics, 558(2013):A33.
  * [BBC11] Stefan Behnel, Robert Bradshaw, Craig Citro, Lisandro Dalcin, Dag Sverre Seljebotn and Kurt Smith. _Cython: The Best of Both Worlds_ , Computing in Science and Engineering, 13, 31-39(2011).
  * [BDT15] Bokeh Development Team. _Bokeh: Python library for interactive visualization_ , http://www.bokeh.pydata.org, last visited: November 10, 2015.
  * [BRh15] Brandon Rhodes. _jplephem 2.5_ , https://pypi.python.org/pypi/jplephem, last visited: November 10, 2015.
  * [DAV13] David A. Vallado, Wayne D. McClain. _Fundamentals of Astrodynamics and Applications_ , 4th Edition, Microcosm Press, 2013.
  * [FML15] Frazer McLean. _Orbital_ , https://github.com/RazerM/orbital, last visited: September 17, 2015.
  * [HEi15] Helge Eichhorn. _Plyades: A Python astrodynamics library_ , http://github.com/helgee/plyades, last visited: September 17, 2015.
  * [IFC15] Intel Corporation. _Intel® Fortran Compiler - Support for Fortran language standards_ , https://software.intel.com/en-us/articles/intel-fortran-compiler-support-for-fortran-language-standards, last visited: September 19, 2015.
  * [JCR15] Juan Luis Cano Rodríguez, Jorge Cañardo Alastuey. _Poliastro: Astrodynamics in Python_ , Zenodo, 2015. doi:10.5281/zenodo.17462.
  * [JDH07] John D. Hunter. _Matplotlib: A 2D Graphics Environment_ , Computing in Science & Engineering, 9, 90-95(2007).
  * [JOP01] Eric Jones, Travis Oliphant, Pearu Peterson, et al. _SciPy: Open source scientific tools for Python_ , http://www.scipy.org/, last visited: November 10, 2015.
  * [NAS15] National Aeronautics and Space Association. _ISS Trajectory Data_ , http://spaceflight.nasa.gov/realdata/sightings/SSapplications/Post/JavaSSOP/orbit/ISS/SVPOST.html, last visited: August 28, 2015.
  * [NDT15] Numba Development Team. _Numba_ , http://numba.pydata.org, last visited: November 10, 2015.
  * [Ore15] CS Systèmes d’Information. _Orekit: An accurate and efficient core layer for space flight dynamics applications_ , http://www.orekit.org, last visited: September 17, 2015.
  * [PGH11] Fernando Perez, Brian Granger, John D. Hunter. _Python: An Ecosystem For Scientific Computing_ , Computing in Science & Engineering 13.2(2011):13-21.
  * [PGu14] Philip Guo. _Python is Now the Most Popular Introductory Teaching Language at Top U.S. Universities_ , ACM Communications, July 7, 2014.
  * [PyA15] Juan Luis Cano Rodriguez, Helge Eichhorn, Frazer McLean. _Python Astrodynamics_ , http://www.python-astrodynamics.org, last visited: September 17, 2015.
  * [RCM08] Robert C. Martin. _Clean Code: A Handbook of Agile Software Craftsmanship_ , Prentice Hall, 2008.
  * [WCV11] Stéfan van der Walt, S. Chris Colbert and Gaël Varoquaux. _The NumPy Array: A Structure for Efficient Numerical Computation_ , Computing in Science & Engineering, 13, 22-30(2011).


[a]Nikolas Korzoun

#  PeV Gamma-ray Astronomy With Panoramic Optical SETI Telescopes

Wystan Benbow  Aaron Brown  Gregory Foote  William F. Hanlon  Olivier Hervet
John Hoang  Jamie Holder  Paul Horowitz  Wei Liu  Jérôme Maire  Nicolas Rault-
Wang  Dan Werthimer  James Wiley  David A. Williams  Shelley A. Wright

###### Abstract

The Panoramic Search for Extraterrestrial Intelligence (PANOSETI) experiment
is designed to detect pulsed optical signals on nanosecond timescales.
PANOSETI is therefore sensitive to Cherenkov radiation generated by extensive
air showers, and can be used for gamma-ray astronomy. Each PANOSETI telescope
uses a 0.5 m Fresnel lens to focus light onto a 1024 pixel silicon
photomultiplier camera that images a 9.9∘$\times$9.9∘ square field of view.
Recent detections of PeV gamma-rays from extended sources in the Galactic
Plane motivate constructing an array with effective area and angular
resolution surpassing current observatories. The PANOSETI telescopes are much
smaller and far more affordable than traditional imaging atmospheric Cherenkov
telescopes (IACT), making them ideal instruments to construct such an array.
We present the results of coincident observations between two PANOSETI
telescopes and the gamma-ray observatory VERITAS, along with simulations
characterizing the performance of a PANOSETI IACT array.

##  1 Introduction

In 2021, LHAASO announced the detection of ultrahigh-energy (UHE) photons ($>$
0.1 PeV) from 12 sources along the Galactic Plane [1]. Of these UHE sources,
only the Crab Nebula was firmly established as a counterpart to one of the
emitters. With the release of their first gamma-ray catalog, LHAASO now report
UHE emission from 43 sources [2]. Observatories like LHAASO and HAWC are
uniquely suited to survey a majority of the sky for UHE emission due to their
wide fields of view, comprehensive duty cycles, and large effective areas [2,
3]. These properties together provide sensitivity to the low flux of UHE
photons (5 photons km-2 day-1 at 0.1 PeV) [1]. Many of the new LHAASO sources
are spatially extended, and 32 are unassociated with other objects (5 of which
are UHE) [2]. Resolving these sources will require targeted studies with
instruments that can provide superior angular resolution to LHAASO and HAWC.
Arrays of IACTs benefit from increased angular resolution at high energies as
a result of stereoscopic reconstruction of the imaged air showers, making them
a suitable alternative to particle detectors. The challenge is to design an
experiment that makes use of the IACT technique and also meets the effective
area requirements in a way that is cost effective and easy to deploy.

We propose an excellent way to achieve this goal is to use the technology
developed by PANOSETI, an experiment designed to search for both optical SETI
targets and optical astrophysical transients in underexplored, short time-
domains. The plan is to simultaneously observe the whole sky on nanosecond
timescales, using many small and easily affordable telescopes [4]. This makes
PANOSETI sensitive to Cherenkov photons produced by atmospheric air showers.
The telescopes house four quadrant boards, each with an array of four
8$\times$8 silicon photomultiplier (SiPM) square pixels for a total camera
size of 32$\times$32 pixels. Each SiPM array (Hamamatsu S13361-3050AE-08) can
resolve single photons on the order of a nanosecond. Custom readout boards
shape the signal for triggering an image with precise timing. Light is focused
on the camera through a 0.5 meter, acrylic Fresnel lens (Orafol SC214) [5].
The total field of view is 9.9∘ on a side, giving a pixel width of roughly
0.31∘ [6]. By positioning a cluster of telescopes together and pointing them
in different directions, a "fly’s eye" dome can view the whole sky at once
[5]. Two such domes separated by about a kilometer would be ideal for
observing pulsed, point-like signals. The same telescope components can
instead be used on steerable mounts separated by a few hundred meters to
observe Chernekov showers. Telescopes are currently installed at two locations
at Lick Observatory in California [6], with two more sites being considered
for a pathfinder array that would study other astrophysical phenomena. Here we
consider UHE gamma-ray astronomy as one such application.

To emphasize how PANOSETI’s design contrasts with contemporary IACTs, consider
the specifications of VERITAS, an array of four IACTs located at the Fred
Lawrence Whipple Observatory (FLWO) in southern Arizona. Each telescope is 12
m in diameter, with a camera composed of 499 photomultiplier tubes that
together image a total field of view of 3.5∘. The VERITAS energy sensitivity
ranges from 100 GeV to >30 TeV [7]. To further appreciate how two disparate
telescopes can facilitate the same science, we test the performance of
PANOSETI telescopes in conjunction with VERITAS.

##  2 Coincident Observing

In November 2021, two PANOSETI telescopes were brought to FLWO to observe the
Crab Nebula simultaneously with VERITAS, shown in Figure 1 [8]. VERITAS
triggers at a rate of 400 Hz, and in 6 hours PANOSETI detected over 10,000 of
the showers seen by VERITAS. More than 3,000 of these were seen by both
PANOSETI telescopes. At the time of observations, a readout of the PANOSETI
camera was triggered only in the quadrant of the camera where any pixel passed
the 11.5 photo-electron trigger threshold (updated firmware triggers all four
quadrants). Coincident events were identified from the timing information
between all telescopes, and a moment analysis was performed on the PANOSETI
images to calculate the Hillas parameters for direct comparison of shower
geometry with VERITAS [9, 10]. Measurements of these parameters in the
PANOSETI image give a width of 0.13∘, while VERITAS measures the width to be
0.138∘. VERITAS could determine that three of the showers seen by both
PANOSETI telescopes were very likely from gamma-rays that originated from the
Crab Nebula, all with energies above 10 TeV [8]. See Figure 2 for a comparison
of a 14.8 TeV gamma-ray initiated shower between PANOSETI and VERITAS.

Figure 1: Two PANOSETI telescopes positioned in front of the T4 VERITAS
telescope.  Figure 2: PANOSETI (left) and VERITAS (right) images of a 14.8 TeV
shower originating from the Crab Nebula. The Crab Nebula is positioned in the
center of the PANOSETI camera, but offset from the center of VERITAS camera at
the location of the red star symbol. (A firmware upgrade now triggers all four
quadrants when an event is detected.)

##  3 Simulation

To design an array of PANOSETI telescopes that is optimized for observing
gamma-ray showers, we have started developing simulation tools. Air showers
are simulated using CORSIKA [11] version 7.7410, compiled with options
typically used for Cherenkov telescopes111ATMEXT, CERENKOV, IACT, SLANT, and
VOLUMEDET and a combination of hadronic interaction models222QGSJET-II-04 [12]
and UrQMD 1.3.1 [13, 14]. One telescope is simulated at each of the four
locations at Lick Observatory mentioned earlier (Figure 3, right), with each
telescope placed at the same altitude (1239 m above sea level) for simplicity.
The atmospheric profile is approximated by measurements taken at VERITAS
during winter. The Cherenkov photons are simulated in bunch sizes of 5 to
reduce simulation time without compromising image quality [15].

Figure 3: Simulated images of a 100 TeV gamma-ray shower as seen by four
PANOSETI telescopes at Lick Observatory (left). The telescopes face North and
point to zenith, at the simulated direction of the incident gamma-ray. The
approximate locations of the telescopes, with an arbitrary choice of origin,
are also plotted (right). The location of the core was reconstructed (blue
$+$) to within 5 meters of the true core position (red $\times$), and the
direction was reconstructed to 0.03 degrees away from zenith.

The Cherenkov photon output is then read by a modified version of
corsikaIOreader [16], which is software used to prepare the CORSIKA output for
optics simulation and for filling histograms using ROOT [17]333See also ”ROOT”
[software], Release v6.28/04, 30/05/2023,
https://root.cern/releases/release-62804/. Atmospheric extinction of Cherenkov
photons is calculated at this stage rather than in the CORSIKA simulation. The
modifications to corsikaIOreader then stochastically attenuate photons to
simulate the transmission efficiency of the Fresnel lens and photon detection
efficiency (PDE) of the SiPM. The probability of attenuation is calculated
using the wavelength dependence of both transmission through PMMA (acrylic)
and the SiPM’s specified PDE. Figure 4 plots the expected number of photo-
electrons detected by a PANOSETI telescope as a function of distance to the
shower core.

Figure 4: Lateral distribution of Cherenkov light captured by simulated
PANOSETI telescopes. Atmospheric extinction and photon detection efficiency is
accounted for. Each point is the average of 100 simulated showers. As
configured while observing at VERITAS, the PANOSETI telescopes create useful
images when more than 60 photo-electrons are detected. Showers with energies
below 10 TeV are not bright enough to discern above NSB and cannot be reliably
reconstructed. Higher energy showers will appear brighter at the same
distance, or can be seen from farther away.

The photon arrival directions are scattered according to the measured optical
point spread of the lens [5]. The arrival directions are then used to create
images of the shower rather than computing a full ray-trace through a model of
the telescope. Night sky background (NSB) is simulated by scaling down the NSB
rate measured by VERITAS. Finally, images are cleaned using an aperture
cleaning method [18]. See Figure 5 for an example of this process.

Figure 5: Cherenkov photons from a simulation of a 100 TeV gamma-ray shower in
the field of view of a PANOSETI telescope. Each plot shows sequentially (from
left to right) the process of adding NSB, applying the aperture cleaning, and
parameterizing the cleaned image.

To reconstruct the shower direction and core location, Hillas parameters are
derived from the cleaned images and used to identify the shower axes. The
intersection of shower axes from each unique image pair are averaged together
and weighted by the angle between them, based on the method used by HEGRA
[19]. The Hillas parameterization and shower reconstruction are implemented
with routines developed for the VERITAS custom analysis package, Eventdisplay
[20]. Figure 3 shows an example of reconstruction.

##  4 Discussion

The observations with VERITAS demonstrate that PANOSETI can image Cherenkov
showers above 10 TeV. The real image shown in Figure 2 was reconstructed by
VERITAS to land 125 meters away from the telescope, and has a total integrated
intensity of 83 photo-electrons in the PANOSETI camera. The intensity
predicted by Figure 4 of a shower with the same energy and impact parameter is
roughly in agreement with what was measured. Observations have also
highlighted useful features which are currently under development, such as
full camera readout and different triggering modes. Requiring two pixels to
pass the trigger threshold instead of one will dramatically reduce the chance
of triggering on a fluctuation due to NSB or instrumental artifact, allowing
us to lower the required pixel intensity and therefore decrease the energy
threshold. Simulations also show that UHE Cherenkov showers can be imaged by
PANOSETI without truncation or camera saturation (Figure 3), which allows for
accurate reconstruction of the arrival direction and energy of the primary
photon.

The primary challenge with studying UHE sources is achieving a sufficient
effective area. For an IACT, this is determined by the size of the Cherenkov
light pool. A single PANOSETI telescope can already detect UHE showers from
hundreds of meters away. The four locations at Lick Observatory (Figure 3,
right) were chosen from sites that have preexisting infrastructure to support
deployment. The spacing between the telescopes was chosen to optimize optical
observations of SETI targets, and to maximize the amount of light Figure 4
predicts would be collected from Cherenkov showers. This pathfinder array is
only the first test, and a small array of 24 telescopes spaced a few hundred
meters apart can easily fill a footprint that covers a square kilometer. Such
an array would provide greater effective area and angular resolution than
current experiments for a fraction of the cost. PANOSETI telescopes are much
smaller and far more affordable than traditional IACTs. Two were successfully
deployed at VERITAS with minimal effort in under a day. An array of PANOSETI
telescopes would provide a unique opportunity to resolve the UHE sources, and
especially those that cannot be seen at lower energies.

##  5 Acknowledgments

We thank the Lick Observatory staff and engineers for their help in the
installation of the PANOSETI telescopes at the Astrograph dome and Barnard
sites. We also thank the VERITAS Collaboration and the technical support staff
at the Fred Lawrence Whipple Observatory for their cooperation in obtaining
joint observations and for the use of their data. The PANOSETI research and
instrumentation program is made possible by the enthusiastic support and
interest by Franklin Antonio and the Bloomfield Family Foundation. Harvard
SETI was supported by The Planetary Society and The Bosack/Kruger Charitable
Foundation. UC Berkeley’s SETI efforts involved with PANOSETI are supported by
NSF grant 1407804, the Breakthrough Prize Foundation, and the Marilyn and
Watson Alberts SETI Chair fund.

## References

  * [1] Z. Cao, F.A. Aharonian, Q. An, Axikegu, L.X. Bai, Y.X. Bai et al., _Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 $\gamma$-ray Galactic sources_, _Nature_ 594 (2021) 33.
  * [2] Z. Cao, F. Aharonian, Q. An, Axikegu, Y.X. Bai, Y.W. Bao et al., _The First LHAASO Catalog of Gamma-Ray Sources_ , _arXiv e-prints_ (2023) arXiv:2305.17030 [2305.17030].
  * [3] A. Albert, R. Alfaro, C. Alvarez, J.R.A. Camacho, J.C. Arteaga-Velázquez, K.P. Arunbabu et al., _3HWC: The Third HAWC Catalog of Very-high-energy Gamma-Ray Sources_ , _The Astrophysical Journal_ 905 (2020) 76.
  * [4] S.A. Wright, P. Horowitz, J. Maire, D. Werthimer, F. Antonio, M. Aronson et al., _Panoramic optical and near-infrared SETI instrument: overall specifications and science program_ , in _Ground-based and Airborne Instrumentation for Astronomy VII_ , C.J. Evans, L. Simard and H. Takami, eds., vol. 10702 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 107025I, July, 2018, DOI [1808.05772].
  * [5] J. Maire, S.A. Wright, M. Cosens, F.P. Antonio, M.L. Aronson, S.A. Chaim-Weismann et al., _Panoramic optical and near-infrared SETI instrument: optical and structural design concepts_ , in _Ground-based and Airborne Instrumentation for Astronomy VII_ , C.J. Evans, L. Simard and H. Takami, eds., vol. 10702 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 107025L, July, 2018, DOI [1808.05773].
  * [6] M. Cosens, J. Maire, S.A. Wright, F. Antonio, M. Aronson, S.A. Chaim-Weismann et al., _Panoramic optical and near-infrared SETI instrument: prototype design and testing_ , in _Ground-based and Airborne Instrumentation for Astronomy VII_ , C.J. Evans, L. Simard and H. Takami, eds., vol. 10702 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 107025H, July, 2018, DOI [1808.05774].
  * [7] N. Park and VERITAS Collaboration, _Performance of the VERITAS experiment_ , in _34th International Cosmic Ray Conference (ICRC2015)_ , vol. 34 of _International Cosmic Ray Conference_ , p. 771, July, 2015, DOI [1508.07070].
  * [8] J. Maire, S.A. Wright, J. Holder, D. Anderson, W. Benbow, A. Brown et al., _Panoramic SETI: program update and high-energy astrophysics applications_ , in _Ground-based and Airborne Instrumentation for Astronomy IX_ , C.J. Evans, J.J. Bryant and K. Motohara, eds., vol. 12184 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 121848B, Aug., 2022, DOI [2210.01356].
  * [9] A.M. Hillas, _Cerenkov Light Images of EAS Produced by Primary Gamma Rays and by Nuclei_ , in _19th International Cosmic Ray Conference (ICRC19), Volume 3_ , vol. 3 of _International Cosmic Ray Conference_ , p. 445, Aug., 1985.
  * [10] D.J. Fegan, _TOPICAL REVIEW: $\gamma$/hadron separation at TeV energies_, _Journal of Physics G Nuclear Physics_ 23 (1997) 1013.
  * [11] D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz and T. Thouw, _CORSIKA: a Monte Carlo code to simulate extensive air showers._ (1998).
  * [12] S. Ostapchenko, _Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model_ , _Phys. Rev. D_ 83 (2011) 014018 [1010.1869].
  * [13] S.A. Bass et al., _Microscopic models for ultrarelativistic heavy ion collisions_ , _Prog. Part. Nucl. Phys._ 41 (1998) 255 [nucl-th/9803035].
  * [14] M. Bleicher et al., _Relativistic hadron hadron collisions in the ultrarelativistic quantum molecular dynamics model_ , _J. Phys. G_ 25 (1999) 1859 [hep-ph/9909407].
  * [15] K. Bernlöhr, _Simulation of imaging atmospheric Cherenkov telescopes with CORSIKA and sim_ telarray_ , _Astroparticle Physics_ 30 (2008) 149 [0808.2253].
  * [16] G. Maier, _corsikaIOreader: a tool to read CORSIKA eventio files_ , Nov., 2020.  10.5281/zenodo.4182562.
  * [17] R. Brun and F. Rademakers, _ROOT: An object oriented data analysis framework_ , _Nucl. Instrum. Meth. A_ 389 (1997) 81.
  * [18] M. Wood, T. Jogler, J. Dumm and S. Funk, _Monte Carlo studies of medium-size telescope designs for the Cherenkov Telescope Array_ , _Astroparticle Physics_ 72 (2016) 11 [1506.07476].
  * [19] W. Hofmann, I. Jung, A. Konopelko, H. Krawczynski, H. Lampeitl and G. Pühlhofer, _Comparison of techniques to reconstruct VHE gamma-ray showers from multiple stereoscopic Cherenkov images_ , _Astroparticle Physics_ 12 (1999) 135 [astro-ph/9904234].
  * [20] G. Maier and J. Holder, _Eventdisplay: An Analysis and Reconstruction Package for Ground-based Gamma-ray Astronomy_ , in _35th International Cosmic Ray Conference (ICRC2017)_ , vol. 301 of _International Cosmic Ray Conference_ , p. 747, July, 2017, DOI [1708.04048].