# Reliability of fault-tolerant system architectures for automated driving
systems

Tim M. Julitz<EMAIL_ADDRESS>Chair for Product Safety and Quality
Engineering, University of Wuppertal, Germany.  Antoine Tordeux  tordeux@uni-
wuppertal.de Chair for Reliability and Traffic Safety, University of
Wuppertal, Germany.  Manuel Löwer<EMAIL_ADDRESS>Chair for Product
Safety and Quality Engineering, University of Wuppertal, Germany.

###### Abstract

Automated driving functions at high levels of autonomy operate without driver
supervision. The system itself must provide suitable responses in case of
hardware element failures. This requires fault-tolerant approaches using
domain ECUs and multicore processors operating in lockstep mode. The selection
of a suitable architecture for fault-tolerant vehicle systems is currently
challenging. Lockstep CPUs enable the implementation of majority redundancy or
M-out-of-N ($M$oo$N$) architectures. In addition to structural redundancy,
diversity redundancy in the ECU architecture is also relevant to fault
tolerance. Two fault-tolerant ECU architecture groups exist: architectures
with one ECU (system on a chip) and architectures consisting of multiple
communicating ECUs. The single-ECU systems achieve higher reliability, whereas
the multi-ECU systems are more robust against dependent failures, such as
common-cause or cascading failures, due to their increased potential for
diversity redundancy. Yet, it remains not fully understood how different types
of architectures influence the system reliability. The work aims to design
architectures with respect to CPU and sensor number, $M$oo$N$ expression, and
hardware element reliability. The results enable a direct comparison of
different architecture types. We calculate their reliability and quantify the
effort to achieve high safety requirements. Markov processes allow comparing
sensor and CPU architectures by varying the number of components and failure
rates. The objective is to evaluate systems’ survival probability and fault
tolerance and design suitable sensor-CPU architectures. The results show that
the system architecture strongly influences the reliability. However, a
suitable system architecture must have a trade-off between reliability and
self-diagnostics that parallel systems without majority redundancies do not
provide.

###### keywords:

Autonomous driving, Advanced driver-assistance system, Fault tolerance, Fail
operational, Hardware architecture, Markov process

##  1 Introduction

The draft law amending the road traffic code and the compulsory vehicle
insurance in Germany takes the development of automated vehicle systems to the
next level (Ludewig and Grieser, 2021). The law creates the conditions for the
use of highly automated vehicles (SAE level 4, On-Road Automated Driving
(ORAD) Committee (2021)) in public road traffic. Already in 2017, the eighth
amendment to the road traffic code came into force, enabling the operation of
level 3 vehicles. The beginnings of research on vehicle automation go back to
the PROMETHEUS project, which started in 1986 and laid the foundations of
today’s commercial driver assistance systems up to level 2 (Williams, 1988).
Some examples are the distance cruise control Distronic Plus and the emergency
brake assistant Pre-Safe from Daimler. The next milestones in the development
of autonomous vehicles were achieved with the launch of the DARPA Challenges
in 2004. In 2005, driverless vehicles made their way over 212 km through the
Mojave Desert, focusing on autonomous navigation (Crane, 2007). In the follow-
up project in 2007, urban traffic was simulated on an abandoned Air Force base
(Buehler et al., 2009).

A large number of software architectures emerged from the DARPA Challenges,
which have one essential thing in common: The architectures are divided into
modules that fulfil different functions. The modules essentially consist of
localisation, perception and vehicle control (Reke et al., 2020). The 2005
winning vehicle team identified some significant problems. The developed
vehicle was able to successfully drive in a static environment, but navigation
through road traffic is not possible due to the insufficient reliability of
the system (Thrun et al., 2006). The hardware architecture of the vehicle
consisted of six computers that performed various functions. Watchdogs
monitored the states of software and hardware to restart the system in case of
failure.

During the 2007 DARPA Challenge, significant successes were achieved in the
monitored urban environment. The first-place vehicle system relied on
different modes of operation: a normal state and a recovery state (Urmson et
al., 2008). The recovery state is triggered when objects block the planned
path, objects are detected too late, or actions are kinematic infeasible. Four
algorithms are used to return to normal operation, which essentially consist
of replanning paths and increasing the safety margin (Urmson et al., 2008).
The software-based solutions allow increasing robustness. The hardware-based
measures consist of a dual-core CPU and the combination of various sensors.
However the reliability and robustness of the vehicle was not sufficient to
drive in real road traffic, which is considerably more complex than the
monitored environment. Although the vehicle was able to recover from many
fault conditions, the time required to do so was considerable, up to ten
minutes, which is not reasonable in the real environment. Therefore, the
developed system cannot be classified with SAE level 4.

In the last ten years, many projects from industrial and academic organisms
have further advanced the state of the art. On the way to SAE Level 4
driverless operation On-Road Automated Driving (ORAD) Committee (2021),
repeated accidents of autonomous vehicles show that further research is needed
to increase the safety, reliability and robustness of the automated driving
systems Daily et al. (2017). The elimination of the human fallback level
requires fault-tolerant approaches to system modelling that cannot be
implemented by simple redundancy. The principles of fault tolerance are based
on self-diagnosis, reliability and availability evaluation, reconstruction and
error recovery. The reliability is usually increased by structural redundancy.
Availability indicates whether a system is functioning at a certain point in
time and can be influenced by diversity redundancy, operation independence, or
asymmetric system architectures. These principles are already applied in
traditional safety-critical systems, which can be found in aviation, rail
transport, space travel, military or nuclear power plants, among others. They
are also becoming increasingly important in the automotive literature, see,
e.g., Baleani et al. (2003); Kohn et al. (2015); Ishigooka et al. (2018);
Schmid et al. (2019); Lin et al. (2018); Sari (2020). The development of a
fault-tolerant system architecture remains nowadays one of the most important
challenges for the market introduction of autonomous vehicles Daily et al.
(2017).

In this contribution, we aim to analyse and compare different types of
redundant and fault-tolerant architectures including parallel and $M$oo$N$
systems using Markovian processes. The contribution is organised as following.
We present the types of considered architecture models in Sec. 2. The
modelling and analysis of the models using Markovian processes are detailed in
Sec.3. Numerical results for different types of architectures are presented in
Sec. 4 and discussed in Sec. 5.

##  2 Architecture models

Many different redundant system architectures exist in the literature. In the
databases of Springer, the IEEE Xplore Digital Library, the Wiley Online
Library and the ACM Digital Library, the keywords ”fail operational”,
”autonomous driving architecture”, ”fallback strategy”, ”$M$oo$N$ redundancy”,
”system on a chip”, ”reliability” and ”functional safety” are currently
frequently used. Indeed, $M$oo$N$ redundancy configurations allow flexible
operation of systems, suitable with recovery and self-diagnose processes.
These characteristics are well adapted to the automation of complex systems,
and especially automation of the driving. In Sari (2020), the hardware
components of a fault-tolerant system architecture is defined as a combination
of microcontroller units (MCUs) supplied by three independent sensors (camera,
radar, lidar), both configurate in a $M$oo$N$ redundancy. The structure of a
single MCU is shown in fig 1 which will be used as a basis for the following.

Monitoring
---
Sensor
---
Processing
---
MCU
---
Watchdog
---
Domain-ECU
---
Sup
---
Fusion
---
Cam
---
Radar
---
Lidar
---
Figure 1: Structure of a MCU in 1oo1 architecture based on Sari (2020)

The MCUs are powered by independent power supplies and are controlled by
independent watchdogs. They consists of three blocks. A monitoring unit, a
sensor fusion and a processing block, which are supplied with information from
the sensors. This is the 1oo1 MCU architecture which is embedded in a Domain-
Electronic-Control-Unit (Domain-ECU) and provides the basis for various
$M$oo$N$ redundancies. In addition, the following distinguishes between a
sensor architecture and a MCU architecture, which is noted as
$S$oo$N_{S}$/$M$oo$N_{M}$. The sensor architecture specifies how many of the
available sensors must be operational to ensure the functionality of the
system.

The selection of the sensor architecture represents a central research topic
in the development of autonomous vehicles Daily et al. (2017). Two types of
architectures are identified:

  1. 1.

System-on-a-Chip architectures (hereafter referred to as 1-ECU system) and

  2. 2.

Dual-Fail-Safe architectures (hereafter referred to as 2oo2DFS system).

The basic 1-ECU architecture includes a 2oo3 redundancy that requires the
operation of four MCUs on one ECU. There are three independent processing
strings for the sensor information (MCU 1 to 3), each supplied by an
independent power supply (Sup 1 to 3). The information comes from the same
sensors (camera, radar, lidar). These are not additionally redundant. The
redundancy lies in the diversification of different sensors. Furthermore, a
majority decision takes place in MCU 4. For this purpose, the outputs of MCUs
1 to 3 are redundantly fed into processor cores 1 to 3 of MCU 4, which execute
the 2oo3 majority comparison Sari (2020). Subsequently, the 2oo3 process is
repeated a second time by core 4 before the result is passed on as output. The
second type of architecture consists of the 2oo2DFS architecture shown in Sari
(2020). The 2oo2DFS architecture consists of two independent fail-safe
subsystems that are outsourced to different ECUs and communicate with each
other in isolation. In the event of a fault, one ECU takes over the function
of the other. In this operating state, the system no longer runs fail-
operationally, but only fail-safe as long as both ECUs are not restored. The
two fail-safe subsystems are operated via two independent MCUs, which process
the sensor information and are in 2oo2 majority comparison with each other. A
third MCU implements this comparison.

##  3 Markovian analyse

The selection of a suitable architecture for fault-tolerant functions is a
current challenging research topic. Various approaches exist. The differences
lie mainly in the number of MCUs, the number of sensors and the degree of
majority redundancy. For example, the 1-ECU architecture can also be operated
with 2 sensors and 4 MCUs in the 2oo4 majority comparison. Using Markovian
processes, we can determine the reliability of the systems, compare different
sensor and MCU architectures, and determine which is the most suitable. The
phase diagram of a Markov chain for a system with three MCUs and three sensors
is outlined in Fig. 2.

Figure 2: Phase diagram of a Markov chain for a system with three MCUs and
three sensors

We consider in the following a system with $N_{M}$ identical MCUs and $N_{S}$
identical sensors. The MCU and sensor times to failure are respectively
identically exponentially distributed with $\lambda_{M}>0$ the failure rate of
a MCU and $\lambda_{S}>0$ the failure rate of a sensor. We denote $P_{m,s}(t)$
the probability that $m=0,\ldots N_{M}$ MCU(s) and $s=0,\ldots,N_{S}$
sensor(s) are operating at time $t$. We have $\sum_{m,s}P_{m,s}(t)=1$ for all
$t$. The definition of the failure of the entire system depends on the
majority redundancy used. A 2oo3 sensor and 3oo4 MCU architecture is
operational when at least two sensors and three MCUs are functional. We
consider in the following that the MCUs and sensors in, respectively,
$M$oo$N_{M}$ and $S$oo$N_{S}$ configurations. This means that the system fails
as soon as the number of operating MCU(s) $m$ is strictly less than $M$ or the
number of operating sensor(s) $s$ is strictly less than $S$. Therefore the
reliability of the system is

$R(t)=\sum_{m=M}^{N_{M}}\sum_{s=S}^{N_{S}}P_{m,s}(t).$ (1)

There exist $(N_{M}+1)(N_{S}+1)$ possibles states. However, only few
transitions from a state to another are possible. Indeed, the time being
continuous, only one sensor or MCU can fail at a given instant. Furthermore a
failed component can not be restored. Therefore only transitions
$(m,s)\mapsto(m-1,s)$, $m>0$, and $(m,s)\mapsto(m,s-1)$, $s>0$, are possible.
Consequently, the transition matrix of the Markov process associated to the
$M$oo$N$ systems is triangular, and mostly sparse. For instance, the
transition matrix $A$ associated to an architecture consisting of three
sensors and three MCUs as described in Fig. 2 is given by

$A=\text{\footnotesize
$\left[\begin{array}[]{cccccccccccc}-2\lambda_{S}-3\lambda_{M}&3\lambda_{M}&0&0&2\lambda_{S}&0&0&0&0&0&0&0\\\
0&-2\lambda_{S}-2\lambda_{M}&2\lambda_{M}&0&0&2\lambda_{S}&0&0&0&0&0&0\\\
0&0&-2\lambda_{S}-\lambda_{M}&\lambda_{M}&0&0&2\lambda_{S}&0&0&0&0&0\\\
0&0&0&-2\lambda_{S}&0&0&0&2\lambda_{S}&0&0&0&0\\\
0&0&0&0&-\lambda_{S}-3\lambda_{M}&3\lambda_{M}&0&0&\lambda_{S}&0&0&0\\\
0&0&0&0&0&-\lambda_{S}-2\lambda_{M}&2\lambda_{M}&0&0&\lambda_{S}&0&0\\\
0&0&0&0&0&0&-\lambda_{S}-\lambda_{M}&\lambda_{M}&0&0&\lambda_{S}&0\\\
0&0&0&0&0&0&0&-\lambda_{S}&0&0&0&\lambda_{S}\\\
0&0&0&0&0&0&0&0&-3\lambda_{M}&3\lambda_{M}&0&0\\\
0&0&0&0&0&0&0&0&0&-2\lambda_{M}&2\lambda_{M}&0\\\
0&0&0&0&0&0&0&0&0&0&-\lambda_{M}&\lambda_{M}\\\
0&0&0&0&0&0&0&0&0&0&0&0\end{array}\right]$}$ (2)

We denote the vector of probabilities of the different possible states

$\begin{array}[]{l}P=(P_{N_{M},N_{S}},P_{N_{M}-1,N_{S}},\ldots,P_{0,N_{S}},\\\
\hskip
18.49411ptP_{N_{M},N_{S}-1},P_{N_{M}-1,N_{S}-1},\ldots,P_{0,0})^{\text{T}}.\end{array}$
(3)

The state probabilities $P$ of the system and, using Eq. (1), its reliability,
are solution of the Kolmogorov backward linear differential equation

$\dot{P}(t)=A^{\text{T}}P(t),\qquad P(0)=P_{0}.$ (4)

Here $P_{0}$ is the initial state of the system. We obtain directly the
solution

$P(t)=e^{A^{\text{T}}t}P_{0}.$ (5)

The calculation of the matrix exponential can be performed iteratively since
the matrix $A$ is triangular. More generally, the matrix exponential can be
obtained using linear algebra and matrix diagonalisation, Laplace transform,
or numerically on a computing software.

By varying the number of MCUs, sensors and their failure rates $\lambda_{M}$
and $\lambda_{S}$ parameters, the system configuration can be adapted with
regard to reliability-related criteria. The aim is to identify the
configuration of the system architecture with the greatest probability of
survival with sufficient self-diagnosis function. The self-diagnosis consists
of selecting a suitable majority redundancy. Systems with a high proportion of
series ($N$oo$N$) and pure parallel systems (1oo$N$) do not have any self-
diagnosis functions without the use of additional hardware designed for this
purpose and are therefore not considered. For cost reasons, the number of
sensors is generally limited to a maximum of three (cf. camera, radar, lidar).
Additional sensors such as ultrasound do not serve to ensure the safety
objective of the driving task. The number of MCUs is limited to a maximum of
four. This makes it possible to realise up to 3oo4 architectures known from
space travel or nuclear power plants.

##  4 Numerical results

0.00.51.01.52.02.53.00.00.20.40.60.81.0$t~{}\text{in}~{}10^{4}~{}$h$R(t)$1oo1/1oo11oo3/1oo31oo3/2oo31oo3/3oo32oo3/1oo32oo3/2oo32oo3/3oo33oo3/1oo33oo3/2oo33oo3/3oo3
Figure 3: Survival probability $R(t)$ of different majority redundant
$S$oo$N_{S}$/$M$oo$N_{M}$ sensor/MCU architectures with $N_{S}=3$ sensors and
$N_{M}=3$ MCUs. $\lambda_{S}=0.1~{}\cdot~{}10^{-4}~{}$h-1 and
$\lambda_{M}=1~{}\cdot~{}10^{-4}~{}$h-1.
0.00.51.01.52.02.53.00.00.20.40.60.81.0$t~{}\text{in}~{}10^{4}~{}$h$R(t)$1oo1/1oo11oo3/1oo41oo3/2oo41oo3/3oo41oo3/4oo42oo3/1oo42oo3/2oo42oo3/3oo42oo3/4oo43oo3/1oo43oo3/2oo43oo3/3oo43oo3/4oo4
Figure 4: Survival probability $R(t)$ of different majority redundant
$S$oo$N_{S}$/$M$oo$N_{M}$ sensor/MCU architectures with $N_{S}=3$ sensors and
$N_{M}=4$ MCUs. $\lambda_{S}=0.1~{}\cdot~{}10^{-4}~{}$h-1 and
$\lambda_{M}=1~{}\cdot~{}10^{-4}~{}$h-1.

Within the Markovian framework, any majority sensors/MCUs configuration may be
considered. In the following, we analyse numerically systems consisting of up
to three sensors and four MCUs (S3M4). We set
$\lambda_{S}=0.1~{}\cdot~{}10^{-4}~{}$h-1 and
$\lambda_{M}=1~{}\cdot~{}10^{-4}~{}$h-1 and assume that all components alre
initially operating. Figure 3 shows the survival probability $R(t)$
(reliability) of the combination of different majority redundant systems of
sensors and MCUs for an architecture consisting of three sensors and three
MCUs. The leading designation $M$oo$N_{S}$ describes the majority redundancy
of the sensors and the trailing designation $M$oo$N_{M}$ the majority
redundancy of the MCUs. A simple series system without redundancy 1oo1/1oo1 is
listed as a reference (black solid curve). Interestingly, the failure rate of
the sensors must be selected lower than the failure rate of the MCUs in order
to achieve higher reliability values. The 3oo3/3oo3 and 2oo3/3oo3 systems have
the lowest reliability. These are systems with a large proportion of series.
The two curves with the highest reliability are assigned to the 1oo3/1oo3 and
the 2oo3/1oo3 system. However, these are also not suitable for fault-tolerant
operation, as they are systems with an exclusive or large parallel share
(1oo$N$). Parallel systems have no self-diagnosis functions. They are not able
to validate their input independently through majority redundancy. The focus
is therefore on the middle curves in the area of the reference line. Of these,
the grey curve (3oo3/1oo3) has the greatest probability of survival. However,
due to the combination of a pure series and a pure parallel system, it must
also be sorted out. Taking into account the self-diagnosis property, the
2oo3/2oo3 architecture of the S3M3 systems can be said to be the most suitable
for fault-tolerant operation. However, it must be noted that the curve falls
short of the 1oo1/1oo1 reference in the average course of the service life.

Fig. 4 shows the same investigation with four MCUs. Once again, the curves in
the upper and lower reliability range correspond to systems with large series
and large parallel parts with low self-diagnosis capacity, which are to be
sorted out. The curves with suitable majority redundancy are again in the
range of the reference system 1oo1/1oo1. This includes the 3oo3/2oo4 and the
1oo3/2oo4 systems, which are also parallel or series systems. It remains
mostly the 2oo3/2oo4 system. Compared to the 2oo3/2oo3 system of the S3M3
architecture, the 2oo3/2oo4 can be said to be even more suitable, since the
crossing with the reference system 1oo1/1oo1 occurs much later.

##  5 Discussion

The introduction of automated or autonomous vehicle systems of automation SAE
levels four and five requires the development of fault-tolerant system
architectures allowing to dispense with the driver’s fallback level. In this
contribution, different fault-tolerant approaches by majority systems were
analysed and system properties were derived that can be used for further
developments. Research so far has focused on separate aspects of fault
tolerance based on self-diagnosis, reliability, availability, reconstruction
and fault recovery. The present analysis combines the areas of self-diagnosis
and reliability.

The highest reliability is achieved by pure parallel systems (1oo$N$).
However, they are not suitable for fault-tolerant applications due to their
lack of self-diagnostic capability. $M$oo$N$ majority redundancies with M $>1$
validate their input by a majority comparison, Whereby in the event of a
fault, subsystems can be switched off and different operating modes can be
used depending on the system state. The development of diagnostic and
monitoring facilities for fault-tolerant operation of 1oo$N$ systems would
advance the development of fault-tolerant system architecture in terms of
reducing complexity and increasing reliability.

2oo4 architectures are considerably more reliable than 2oo3 architectures.
2oo4 systems are mainly established in the military sector, in space travel
and in nuclear power plants. For the automotive sector, however, mainly only
2oo2 and 2oo3 architectures are discussed. It is recommended that future
research should increasingly consider 2oo4 systems. An important question to
be addressed is whether the reliability gain can be justified in comparison to
the increased costs due to the increased number of components and in
comparison to the greater hardware overhead. If a system based on a single
domain ECU is used, then a 2oo3 sensor/ 2oo4 MCU architecture or a 2oo2
sensor/ 2oo4 MCU architecture provides the best performance in terms of
reliability, taking into account its self-diagnostic capability. Systems with
a large parallel component deliver the highest reliability and systems with a
high serial component deliver the lowest reliability. The reliability changes
only marginally between sensor architectures 2oo2 and 2oo3. From a reliability
point of view, it therefore plays a subordinate role whether two or three
sensors are used.

The components of the system architectures can be divided into three groups
for their design to achieve reliability target values (e.g., ISO 26262 ISO
26262:2018 (2018)), which is essentially based on the order of the minimum
sections associated with the components. In the application case in question,
a factor of approx. 10 between the failure rate of components of the sensor
and the MCU architecture has proven itself. This finding is transferable to
further system designs. It should be noted that the minimum cuts cannot be
blindly used as a basis, since various effects can lead to a component being
more critical with regard to a system failure than the other representatives
of its cut-set order. For this purpose, the use of importance analyses can be
recommended.

For further work, it makes sense to also deal with the other aspects of
availability, reconstruction and fault recovery. The DFS architecture
presented is ideally suited for this purpose. To reduce complexity, the two
ECUs of the DFS system have been actively redundant. However, real fault
tolerance can only be achieved through passive redundancy or stand-by
redundancy. This way, the aspect of reconstruction is taken into account. In
the event of a fault, ECU 1 switches off and ECU 2 takes over operation.
Meanwhile, error correction algorithms (e.g. error correcting codes) restore
the failed components of ECU 1, thus taking the aspect of error correction
into account. The time to repair a component can be implemented as a repair
rate, which allows availability considerations. Through this approach, all
aspects of fault tolerance can be covered.

## References

  * Baleani et al. (2003) Baleani, M., A. Ferrari, L. Mangeruca, A. Sangiovanni-Vincentelli, M. Peri, and S. Pezzini (2003).  Fault-tolerant platforms for automotive safety-critical applications.  In Proceedings of the 2003 international conference on Compilers, architecture and synthesis for embedded systems, pp. 170–177.
  * Buehler et al. (2009) Buehler, M., K. Iagnemma, and S. Singh (2009).  The DARPA urban challenge: autonomous vehicles in city traffic, Volume 56.  springer.
  * Crane (2007) Crane, C. D. (2007).  The 2005 DARPA Grand challenge.  In 2007 International Symposium on Computational Intelligence in Robotics and Automation, pp. nil4–nil4. IEEE.
  * Daily et al. (2017) Daily, M., S. Medasani, R. Behringer, and M. Trivedi (2017).  Self-driving cars.  Computer 50(12), 18–23.
  * Ishigooka et al. (2018) Ishigooka, T., S. Honda, and H. Takada (2018).  Cost-effective redundancy approach for fail-operational autonomous driving system.  In 2018 IEEE 21st International Symposium on Real-Time Distributed Computing (ISORC), pp. 107–115. IEEE.
  * ISO 26262:2018 (2018) ISO 26262:2018 (2018).  Road vehicles — Functional safety.  Standard, International Organization for Standardization, Geneva, CH.
  * Kohn et al. (2015) Kohn, A., M. Käßmeyer, R. Schneider, A. Roger, C. Stellwag, and A. Herkersdorf (2015).  Fail-operational in safety-related automotive multi-core systems.  In 10th IEEE International Symposium on Industrial Embedded Systems (SIES), pp. 1–4. IEEE.
  * Lin et al. (2018) Lin, S.-C., Y. Zhang, C.-H. Hsu, M. Skach, M. E. Haque, L. Tang, and J. Mars (2018).  The architectural implications of autonomous driving: Constraints and acceleration.  In Proceedings of the Twenty-Third International Conference on Architectural Support for Programming Languages and Operating Systems, pp. 751–766.
  * Ludewig and Grieser (2021) Ludewig, J. and G. Grieser (2021).  Entwurf eines Gesetzes zur Änderung des Straßenverkehrsgesetzes und des Pflichtversicherungsgesetzes - Gesetz zum autonomen Fahren.  Technical report, No 19/27439, Bundesministerium für Verkehr und digitale Infrastruktur (German Federal Ministry for Digital and Transport).
  * On-Road Automated Driving (ORAD) Committee (2021) On-Road Automated Driving (ORAD) Committee (2021).  Taxonomy and Definitions for Terms Related to Driving Automation Systems for On-Road Motor Vehicles.
  * Reke et al. (2020) Reke, M., D. Peter, J. Schulte-Tigges, S. Schiffer, A. Ferrein, T. Walter, and D. Matheis (2020).  A self-driving car architecture in ros2.  In 2020 International SAUPEC/RobMech/PRASA Conference, pp. 1–6. IEEE.
  * Sari (2020) Sari, B. (2020).  Fail-operational Safety Architecture for ADAS/AD Systems and a Model-driven Approach for Dependent Failure Analysis.  Springer Nature.
  * Schmid et al. (2019) Schmid, T., S. Schraufstetter, S. Wagner, and D. Hellhake (2019).  A safety argumentation for fail-operational automotive systems in compliance with iso 26262.  In 2019 4th International Conference on System Reliability and Safety (ICSRS), pp. 484–493. IEEE.
  * Thrun et al. (2006) Thrun, S., M. Montemerlo, H. Dahlkamp, D. Stavens, A. Aron, J. Diebel, P. Fong, J. Gale, M. Halpenny, G. Hoffmann, et al. (2006).  Stanley: The robot that won the DARPA Grand Challenge.  Journal of field Robotics 23(9), 661–692.
  * Urmson et al. (2008) Urmson, C., J. Anhalt, D. Bagnell, C. Baker, R. Bittner, M. Clark, J. Dolan, D. Duggins, T. Galatali, C. Geyer, et al. (2008).  Autonomous driving in urban environments: Boss and the urban challenge.  Journal of Field Robotics 25(8), 425–466.
  * Williams (1988) Williams, M. (1988).  PROMETHEUS-The European research programme for optimising the road transport system in Europe.  In IEE Colloquium on Driver Information, pp. 1–1. IET.


# Phonon Chirality Induced by Vibronic-Orbital Coupling

Yun-Yi Pai<EMAIL_ADDRESS>Materials Science and Technology Division, Oak
Ridge National Laboratory, Oak Ridge, TN 37831, USA  Quantum Science Center,
Oak Ridge, Tennessee 37831, USA  Claire E. Marvinney  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Quantum Science Center, Oak Ridge, Tennessee 37831, USA  Liangbo Liang  Center
for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN
37831, USA  Ganesh Pokharel  Materials Department and California Nanosystems
Institute, University of California Santa Barbara, Santa Barbara, CA 93106,
USA  Jie Xing  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  Haoxiang Li  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Xun Li  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  Michael Chilcote  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Quantum Science Center, Oak Ridge, Tennessee 37831, USA  Matthew Brahlek
Materials Science and Technology Division, Oak Ridge National Laboratory, Oak
Ridge, TN 37831, USA  Quantum Science Center, Oak Ridge, Tennessee 37831, USA
Lucas Lindsay  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  Hu Miao  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Athena S. Sefat  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  David Parker  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Stephen D. Wilson  Materials Department and California Nanosystems Institute,
University of California Santa Barbara, Santa Barbara, CA 93106, USA  Benjamin
J. Lawrie<EMAIL_ADDRESS>This manuscript has been authored by UT-
Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of
Energy (DOE). The US government retains and the publisher, by accepting the
article for publication, acknowledges that the US government retains a
nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for US
government purposes. DOE will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
(http://energy.gov/downloads/doe-public-access-plan).  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Quantum Science Center, Oak Ridge, Tennessee 37831, USA

###### Abstract

The notion that phonons can carry pseudo-angular momentum has become popular
in the last decade, with recent research efforts highlighting phonon
chirality, Berry curvature of phonon band structure, and the phonon Hall
effect. When a phonon is resonantly coupled to a crystal electric field
excitation, a so-called vibronic bound state forms. Here, we observe angular
momentum transfer of $\Delta\textbf{J}_{z}=\pm 1\hbar$ between phonons and an
orbital state in a vibronic bound state of a candidate quantum spin liquid.
This observation has profound implications for the engineering of phonon band
structure topology through chiral quasiparticle interactions.

##  I Introduction

Since the seminal work of Bloch in 1929[1], electron-phonon coupling has
emerged as one of the most important topics in modern condensed matter
physics. In conventional superconductors, electron-phonon coupling underlies
Cooper pairing. In semiconductors, it sets an upper bound for electron
mobility. Electron-phonon coupling also governs numerous thermal and spin
relaxation processes in solids. In ionic crystals with f orbitals, spin-orbit
coupling and the crystal electric field (CEF) associated with the ionic or
ligand environment split the electronic wavefunction eigenstates into
manifolds. Atomic displacements can change the ligand environments and hence
the orbital manifolds. This change is typically treated within the adiabatic
Born-Oppenheimer approximation, in which atomic motion is static relative to
the electronic degrees of freedom. However, when phonon-CEF coupling is nearly
resonant, the Born-Oppenheimer approximation is inadequate and a new state – a
vibronic bound state (VBS), with both orbital character and phonon character –
can form.

Phonons themselves can have chirality and carry pseudo angular momentum[2, 3,
4, 5, 6, 7] when inversion symmetry ($\mathcal{P}$) or time reversal symmetry
($\mathcal{T}$) is broken[8, 9]. The ability to control pseudo angular
momentum transfer between phonons and electronic degrees of freedom will
unlock new opportunities in information processing. However, the transfer of
angular momentum between photons, electronic spin, orbital excitations, and
phonons is still not well understood. It was studied, for instance, in the
context of paramagnetic spin relaxation[10, 11], but phonons are not generally
thought to carry macroscopic angular momentum[11, 12]. Phononic pseudo angular
momentum is key to several fundamental effects in physics, including the
microscopic explanation of the phonon Hall effect[13], Einstein-de Haas effect
[14], possible circular heat flow[15], and zero point energies of chiral
phonons [14]. The transfer channel from spin to phonon, on the other hand, is
thought to be relevant in ultra-fast demagnetization processes [16, 17].

Here, we report a change of angular momentum $\Delta\textbf{J}_{z}=\pm 1\hbar$
in an orbital excitation due to phononic coupling in NaYbSe2. NaYbSe2[18, 19,
20, 21, 22] belongs to the family of Yb delafossites with the form
A1+Yb3+X${}^{2-}_{2}$ that have been identified as quantum spin liquid (QSL)
candidates. The QSL is a theoretically proposed phase that is characterized by
massive many-body entanglement[23, 24, 25] with possible application to fault-
tolerant quantum information processing[26]. So far 14 members within this
family have been identified to date[27, 28, 29, 21, 20, 19, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39]. They all have planar triangular lattices with spin
S${}_{\text{eff}}=1/2$ on trigonally distorted YbX6 octahedra.
Antiferromagnetic coupling in the triangular lattice causes geometric
frustration, resulting in a lack of long-range magnetic order down to the
lowest probed temperatures. This family of candidate QSLs is objectively less
defect prone than Yb(Mg, Ga)O4[21], and the breadth of substitutional
composites in the family makes it an ideal platform for the study and control
of QSL excitations.

The effective spin S${}_{\text{eff}}$ = 1/2 of the system comes from the Yb3+
ion, which has the [Xe]4f13 electronic configuration. Yb3+, like all the
elements in the 4f block, has weak exchange coupling and strong spin-orbit
coupling compared to the 3d-block elements. The ground state spin-orbit
manifold of Yb3+ has total spin $J=7/2$. The next manifold $J=5/2$ is more
than an eV above the ground state manifold[40, 41]. The ground-state $J=7/2$
manifold splits into 4 time-reversed Kramers pairs. This degeneracy is a
result of Kramers theorem: trivalent Yb3+ has an odd number of electrons and a
$\mathcal{T}$-symmetric Hamiltonian, so the eigenenergies are doubly
degenerate at zero magnetic field. The CEF splitting within the ground-state
manifold is typically only tens of meV due to the well-shielded nature of the
4f electrons.

VBSs form when strong CEF-phonon coupling is present, but they have only been
reported in a handful of materials[42, 43, 44, 45, 46, 47, 48, 49, 50, 51].
Structurally, NaYbSe2 has less distortion in its YbSe6 octahedra than other Yb
delafossites[27]. Its ground states are reported to be described by spinon
excitations [19] or ferrimagnetic quasistatic and dynamic excitations within a
QSL matrix[52]. While pristine NaYbSe2 is insulating, with increasing pressure
it exhibits increased conductivity and eventually demonstrates a dip in
resistivity attributed to a possible superconducting state[18, 22].

##  II Results and Discussion

Figure 1: (a) Energy diagram of the ground state $J=7/2$ manifold CEF levels,
Raman active phonon modes, and vibronic bound state $\omega$. (b) A
semiclassical representations of the CEF eigenstates $|\psi^{\pm}_{0}\rangle$
and $|\psi^{\pm}_{1,2,3}\rangle$. The eigenstates are linear combinations of
multiplets $m_{J}=-7/2\dots 7/2$. For example,
$|\psi^{+}_{0}\rangle=|\frac{7}{2},\frac{3}{2}\rangle$ is represented by
precessing cones with only $c_{\frac{3}{2}}\neq 0$ for
$|\psi^{+}_{0}\rangle=c_{\frac{3}{2}}|J,\frac{3}{2}\rangle$. Hence all but the
$m_{J}=\frac{3}{2}$ cone are transparent (and all but the $m_{J}=\frac{3}{2}$
circle are empty). The same follows for $|\psi^{-}_{0}\rangle$ and
$|\psi^{\pm}_{1,2,3}\rangle$. (c) The E${}_{\text{2g}}$ eigenspace is spanned
by the basis {$\text{E}_{\text{2g, x}}$, $\text{E}_{\text{2g, y}}$}, which can
also be spanned by {$\text{E}_{\text{2g,+}}$, $\text{E}_{\text{2g,-}}$} =
{$\text{E}_{\text{2g, x}}+i\text{E}_{\text{2g, y}}$, $\text{E}_{\text{2g,
x}}-i\text{E}_{\text{2g, y}}$}. (d) Temperature dependent Raman spectra from
$T=3.3$ K to $T=260$ K. (e) The contour plot of the temperature dependent data
from (d). (f), (g) Peak positions and widths for E${}_{\text{2g}}$ and CEF1
from Bayesian inference. The symbols are medians of the posterior
distribution, and the 68% and 95% highest density intervals (HDIs) are
represented by darker and lighter shading, respectively.

Figure 1a shows the low energy excitations of NaYbSe2, where the ground state
spin-orbit manifold $J=7/2$ splits into 4 time-reversed Kramers pairs
$|\psi_{0}^{\pm}\rangle$, $|\psi_{1}^{\pm}\rangle$, $|\psi_{2}^{\pm}\rangle$,
$|\psi_{3}^{\pm}\rangle$. CEF1, CEF2, and CEF3 describe the transition between
$|\psi_{1}^{\pm}\rangle$, $|\psi_{2}^{\pm}\rangle$, $|\psi_{3}^{\pm}\rangle$
and the ground state $|\psi_{0}^{\pm}\rangle$, respectively. The CEF modes are
observed at 117.2 cm-1 (CEF1), 197.8 cm-1 (CEF2), and 247.0 cm-1 (CEF3) in
Raman spectra acquired at 3.3 K, consistent with previous reports[21, 37].
NaYbSe2 also has two Raman-active phonon modes: E${}_{\text{2g}}$ at 124.4
cm-1 and A${}_{\text{1g}}$ at 172.8 cm-1. Figure 1b shows a semiclassical
spinning top representation of the spin-orbit eigenstates for
$|\psi_{0}^{\pm}\rangle$, $|\psi_{1}^{\pm}\rangle$, $|\psi_{2}^{\pm}\rangle$,
$|\psi_{3}^{\pm}\rangle$ as
$\sum_{m_{j}=-\frac{7}{2},...\frac{7}{2}}c_{m_{j}}|\frac{7}{2},m_{j}\rangle$.
Transparent cones represent $c_{m_{j}}=0$. The blue (pink) cones represent the
$+$ ($-$) branch.

Figure 2: (a) Helicity resolved Raman spectra at $T=4$ K, $B=0$ T. Four
prominent modes are CEF1 (with E${}_{\text{2g}}$ at its shoulder), $\omega$,
CEF2, and CEF3. The Rayleigh scattering and $\omega$ are of ($\sigma^{+}$,
$\sigma^{+}$) and ($\sigma^{-}$, $\sigma^{-}$) co-circular scattering
channels, while CEF1-3 are of ($\sigma^{+}$, $\sigma^{-}$) and ($\sigma^{-}$,
$\sigma^{+}$) cross-circular channels. (b-d) Selection rules for CEF modes and
$\omega$ with virtual states illustrated in light pink.

In this material, the doubly degenerate E${}_{\text{2g}}$ mode is resonant
with the CEF1 state. This resonance gives rise to a VBS, $\omega$ (Figure 1a).
As illustrated in Figure 1c, the E${}_{\text{2g}}$ mode is spanned by the
basis {$\text{E}_{\text{2g, x}}$, $\text{E}_{\text{2g, y}}$}; it can also be
spanned by the chiral basis {$\text{E}_{\text{2g,+}}$,
$\text{E}_{\text{2g,-}}$} = {$\text{E}_{\text{2g, x}}+i\text{E}_{\text{2g,
y}}$, $\text{E}_{\text{2g, x}}-i\text{E}_{\text{2g, y}}$}[53]. This second
basis is the natural one for CrBr3 because of broken
$\mathcal{T}$-symmetry[9]. In NaYbSe2, however, because both $\mathcal{P}$ and
$\mathcal{T}$ are intact, phonons are not intrinsically chiral, and the chiral
basis is not, at first glance, the appropriate choice of basis.

Using temperature dependent Raman spectra, we first verify the primary CEF
excitations and phonon modes. The spectra are shown in Figures 1d (traces) and
1e (contour). The CEF modes are clearly identified at 117.2 cm-1 (CEF1), 197.8
cm-1 (CEF2), and 247.0 cm-1 (CEF3) in the spectra. Additionally, they all
become significantly stronger in intensity, and soften (shift toward lower
energy) as the temperature decreases. On the other hand, the temperature
dependence of the E${}_{\text{2g}}$ and the A${}_{\text{1g}}$ phonon modes is
far less drastic. These CEF and phonon modes are consistent with recently
reported linearly polarized NaYbSe2 Raman spectra[21]. Note that the
E${}_{\text{2g}}$ and CEF1 modes are closely spaced. This is the case for
NaYbSe2[21], CsYbSe2[37], and KYbSe2[39]. Using Bayesian inference with a
Hamiltonian Monte Carlo PyMC3[54], we track the peak parameters from data in
Figure 1d. The extracted peaks and widths for CEF1 and E${}_{\text{2g}}$ are
shown in Figure 1f and 1g, respectively. The symbols represent the median
values from the Bayesian inference. The two shaded bands illustrate the 68%
(darker) (corresponding to 1$\sigma$ in the central limit) and 95% (lighter)
(corresponding to 2$\sigma$) highest density intervals (HDIs) (See
Supplementary Information for other selected modes). The model has large error
bars from $T=175$ K to $225$ K, but the two modes are well resolved elsewhere.
The rapid change in the E${}_{\text{2g}}$ Raman mode energy at 90 K and the
associated overlap of this mode with the CEF1 mode from 100 to 200 K is
noteworthy. In CsYbSe2, there is some evidence for anisotropic negative
thermal expansion in these temperature ranges[55, 33]. In general, a nearly
discontinuous change with temperature of any spectral parameter (such as the
energy of the E${}_{\text{2g}}$ mode) is suggestive of a phase transition,
though we find no direct evidence of any phase transition here.

At the same time, we see the VBS mode $\omega$ at 151.2 cm-1. It has
characteristics in common with both CEFs and phonons: $\omega$ becomes
significantly stronger at low temperatures just as the CEFs do, and it hardens
as the temperature decreases like A${}_{\text{1g}}$ does. The mode $\omega$ is
present in both NaYbSe2 and CsYbSe2, although it is far less prominent in
CsYbSe2. This is likely because the spacing between CEF1 and E${}_{\text{2g}}$
is smaller in NaYbSe2 (7.2 cm-1) than in CsYbSe2 (12.9 cm-1)[37]. Fitting to
the Thalmeier-Fulde description of a magnetoelastic vibronic bound state[46]
yields a coupling strength of 32.0 cm-1 (3.97 meV) for NaYbSe2, which is far
stronger than the 23.6 cm-1 (2.93 meV) reported for CsYbSe2 and comparable to
the roughly 34.0 cm-1 (4.22 meV) coupling strength reported for Ce2O3[43].

Typically, CEF spectroscopy data are fit to a CEF Hamiltonian with Stevens
operators given by the symmetry of the ionic environment of the 4f orbitals
(Methods). This links the CEF peaks to eigenstates $|\psi\rangle$. The
procedure is typically under-constrained, but that is not the case when
optical selection rules are taken into account. Figure 2a illustrates
helicity-resolved Raman spectra acquired at $T=4$ K, with CEF1, CEF2, CEF3 and
$\omega$ highlighted. While Rayleigh scattering is stronger in the co-circular
polarization configuration, CEF1–CEF3 are stronger in the cross-circular
polarization configuration. The mode $\omega$ is also stronger in the co-
circular polarization configuration.

Each right (left) circularly polarized $\sigma^{+}$ ($\sigma^{-}$) photon
carries angular momentum $+\hbar$ ($-\hbar$). The link between helicity and
circular polarization is given by $h$ = $\sigma$ $\cdot$ k, where k is the
momentum of the photon, and $h$ the helicity[56]. Absorption of a $\sigma^{+}$
photon increases the angular momentum of the lattice by $+\hbar$,
corresponding to being acted on by operator $J_{+}$, while scattering a
$\sigma^{+}$ photon corresponds to $J_{-}$ acting on the lattice. Therefore,
if one considers Raman spectra acquired in the ($\sigma_{\text{incident}}$,
$\sigma_{\text{scatter}}$) = ($\sigma^{+}$, $\sigma^{-}$) configuration, as in
Figure 2b, the system obtains $+2\hbar$ of angular momentum, corresponding to
$J_{+}J_{+}$ acting on it (a generalization of the closure condition in Axe et
al. [57], which contracts
$\sum_{\text{virtual}}||\langle\psi_{\text{final}}|M_{i}|\psi_{\text{virtual}}\rangle\langle\psi_{\text{virtual}}|M_{j}|\psi_{\text{initial}}\rangle||$
to $||\langle\psi_{\text{final}}|M_{i}M_{j}|\psi_{\text{initial}}\rangle||$
for any dipolar transition $M_{i}$, $M_{j}$). Therefore, a mode that is active
in the ($\sigma^{+}$, $\sigma^{-}$) configuration must have dominant matrix-
element contributions from
$||\langle\psi_{\text{final}}|J_{+}J_{+}|\psi_{\text{initial}}\rangle||$. So
the transition that accounts for the CEF1 mode, which is active in
($\sigma^{+}$, $\sigma^{-}$), must have dominant matrix elements linked by
$\langle 1|J_{+}J_{+}|0\rangle$ for a ground state wavefunction $|0\rangle$
and first excited state wavefunction $|1\rangle$.

Since all the CEF modes are enhanced in the cross-circular channel and
suppressed in the co-circular Raman spectra, $\langle 1|J_{+}J_{+}|0\rangle$,
$\langle 1|J_{-}J_{-}|0\rangle$, $\langle 2|J_{+}J_{+}|0\rangle$, $\langle
2|J_{-}J_{-}|0\rangle$, $\langle 3|J_{+}J_{+}|0\rangle$, $\langle
3|J_{-}J_{-}|0\rangle$ should be large while $\langle 1|J_{+}J_{-}|0\rangle$,
$\langle 1|J_{-}J_{+}|0\rangle$, $\langle 2|J_{+}J_{-}|0\rangle$, $\langle
2|J_{-}J_{+}|0\rangle$, $\langle 3|J_{+}J_{-}|0\rangle$, $\langle
3|J_{-}J_{+}|0\rangle$ should all be small or zero. Without the constraints
imposed by these selection rules, CEF parameters consistent with the
experimentally observed CEF Raman modes can be obtained regardless of the
assignment of $|\psi_{0}\rangle$, $|\psi_{1}\rangle$, $|\psi_{2}\rangle$, and
$|\psi_{3}\rangle$ to $|0\rangle$, $|1\rangle$, $|2\rangle$, and $|3\rangle$.
However, no parameters satisfy both the energy levels and the helicity
selection rules unless $\psi^{\pm}_{0}$ is the ground state $|0\rangle$ and
the remaining modes are defined as excited states, as shown in Figure 1b (see
the Supplementary Information for further discussion). In contrast with the
CEF modes, the $\omega$ mode requires $J_{+}J_{-}$ and $J_{-}J_{+}$. Because
$\omega$ comes from the parent state CEF1 and E${}_{\text{2g}}$, we consider
the 2$\times$2 dimensional space
$|\psi_{1}^{\pm}\rangle\otimes|\text{E}_{\text{2g}}\rangle$. Because of the
strong extinction ratio between the co-circular and cross-circular channel of
$\omega$, it may seem natural to conclude that $\Delta\textbf{J}_{z}=0$ for
the system before and after the interaction with the photon. But in fact, the
angular momentum change $\Delta\textbf{J}_{z}$ between the ground state and
CEF1 equals $\pm 2\hbar$ , as shown in Figure 2d, so there has to be
additional angular momentum transfer due to the phonon degree of freedom.
Therefore, the chiral basis {$\text{E}_{\text{2g,+}}$,
$\text{E}_{\text{2g,-}}$} is the natural choice for describing the
eigenvibration. For this 2 $\times$ 2 space, however, only 2 states,
$\omega^{+}=$ $|\psi_{1}^{-}\rangle\otimes|\text{E}_{\text{2g,+}}\rangle$ and
$\omega^{-}=$ $|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$, can
be connected from the ground state by $J_{+}J_{-}$ or $J_{-}J_{+}$. For
$|\psi_{1}^{\pm}\rangle$, the extracted $c_{1/2}$ is much smaller than
$c_{7/2}$ and $c_{-{5}/{2}}$ (Methods), so we attribute the observed peaks to
$|\psi^{+}_{0}\rangle\to|\omega^{-}\rangle$ or
$|\psi^{-}_{0}\rangle\to|\omega^{+}\rangle$.

Figure 3: Helicity-resolved magnetic field dependence of CEFs and $\omega$ for
NaYbSe2 at $T=4$ K. (a) CEF1. (b) $\omega$. (c) CEF2. (d) CEF3. The subplots
(i), (ii) are contour plots for helicity-resolved magneto Raman with
polarization configurations that the given mode is active in. (iii) and (iv)
are linecuts at $B$ = 0 T, $B$ = -6 T and $B$ = +6 T. (v) peak positions
extracted from (i) and (ii). The 68% and 95% HDIs are represented by darker
and lighter shading, respectively.

To distinguish the $+$ and $-$ branches, we perform magnetic field-dependent
Raman spectroscopy, which lifts the degeneracy within each of the Kramers
pairs and therefore provides further constraints on the CEF assignments.
Magnetization and specific heat measurements have been previously used to show
that no field-induced ordering is present in NaYbSe2 for
$\textbf{B}\;||\;\textbf{c}$ for $B<9$ T at $T=4$ K [20]. Therefore, with the
6 T magnetic fields accessible in this setup, no field-induced spin ordering
is expected. Figure 3 shows the magnetic field dependence of the CEF levels
for (a) CEF1, (b) $\omega$, (c) CEF2, and (d) CEF3 measurements. Sub-panels
(i, ii) are contour plots and (iii, iv) are linecuts at
$B=0~{}\text{T},\;-6~{}\text{T},$ and $6~{}\text{T}$. Again, the CEF levels
only show up in cross-circular polarization configurations,
($\sigma^{+},\sigma^{-}$) and ($\sigma^{-},\sigma^{+}$), while the $\omega$
mode is stronger in co-circular polarization configurations,
($\sigma^{+},\sigma^{+}$) and ($\sigma^{-},\sigma^{-}$). In the
($\sigma^{+},\sigma^{-}$) configuration and with increasing positive magnetic
field, CEF1 and CEF2 increase in energy while CEF3 decreases. This dependence
is inverted if the magnetic field is reversed ($B\to-B$) or the helicity is
reversed (($\sigma^{+},\sigma^{-}$) $\to$ ($\sigma^{-},\sigma^{+}$)).

Figure 4: From helicity-resolved field dependence, transitions between the $+$
and $-$ Kramers pairs in the eigenstates $|\psi^{\pm}_{0}\rangle$ and
$|\psi^{\pm}_{1,2,3}\rangle$ can be assigned (a) CEF1 (b) CEF2 and (c) CEF3.
(d) and (e) $|\psi_{1}^{-}\rangle\otimes|\text{E}_{\text{2g,+}}\rangle$,
$|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$.

Sub-panel (v) of Figure 3a–d shows the Bayesian-inferred peak positions for
CEF1 (a), $\omega$ (b), CEF2 (c), and CEF3 (d) from the co- and cross-
polarized spectra shown in sub-panels (i) and (ii). The aforementioned trends
for the CEF modes are well-resolved in (v) while $\omega$ shows no field
dependence.

These constraints allow us to make the following assignments: for cross-
circular polarized ($\sigma^{+},\sigma^{-}$), the observed CEF1 comes from the
transition $|\psi^{+}_{0}\rangle\to|\psi^{+}_{1}\rangle$, CEF2 comes from
$|\psi^{+}_{0}\rangle\to|\psi^{+}_{2}\rangle$, and CEF3 comes from
$|\psi^{-}_{0}\rangle\to|\psi^{+}_{3}\rangle$. Similarily, for cross-circular
polarized ($\sigma^{-},\sigma{+}$), the observed CEF1 comes from the
transition $|\psi^{-}_{0}\rangle\to|\psi^{-}_{1}\rangle$, CEF2 comes from
$|\psi^{-}_{0}\rangle\to|\psi^{-}_{2}\rangle$, and CEF3 comes from
$|\psi^{+}_{0}\rangle\to|\psi^{-}_{3}\rangle$. These assignments not only
satisfy the aforementioned selection rule but also capture the correct sign of
the magnetic field dependence when an additional Zeeman dependence
$H_{\textbf{B}}=-g_{J}\mu_{B}\textbf{J}_{z}\cdot\textbf{B}$ is added. Figure
4a–c shows the relative transition probability using these state assignments.
The transition probabilities are consistent with the previously established
set of selection rules and states (e.g., CEF1 is dominated by
$|\psi^{+}_{0}\rangle\to|\psi^{+}_{1}\rangle$ with $J_{+}J_{+}$ and
$|\psi^{-}_{0}\rangle\to|\psi^{-}_{1}\rangle$ with $J_{-}J_{-}$).

Because the $\omega^{\pm}$ show negligible magnetic field dependence, we
assign the observed transitions to
$|\psi_{1}^{+}\rangle\to|\psi_{1}^{-}\rangle\otimes|\text{E}_{\text{2g,+}}\rangle$
(as shown in Figure 2d) and
$|\psi_{1}^{-}\rangle\to|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$.
Because the transition from $|\psi_{1}^{+}\rangle$ to the bare
$|\psi_{1}^{-}\rangle$ requires $\Delta J_{z}=+1~{}\hbar$ but the final change
is $\Delta J_{z}=0$, the phonon subsystem takes $+1\hbar$ and creates one
chiral phonon $|\text{E}_{\text{2g,+}}\rangle$. Similarly, a chiral phonon
$|\text{E}_{\text{2g,-}}\rangle$ is created in the transition
$|\psi_{1}^{-}\rangle\to|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$.
Because the parent space
$|\psi_{1}^{\pm}\rangle\otimes|\text{E}_{\text{2g}}\rangle$ is 2$\times$2
dimensional, two more states
$\overline{\omega}^{+}=|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,+}}\rangle$
and
$\overline{\omega}^{-}=|\psi_{1}^{-}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$
are expected, just as in the original Thalmeier-Fulde hybridization model for
a VBS[46]. These two states may be absent in the measured spectra because they
require non-negligible $c_{1/2}$ for $|\psi_{1}^{\pm}\rangle$ to have a
meaningful matrix element, as shown in the right panel of Figure 2d. This
observation is surprising because NaYbSe2 is $\mathcal{P}$-symmetric and
$\mathcal{T}$-symmetric and because the coupling of two
$\mathcal{T}$-symmetric degenerate subsystems (CEF Kramers pair
$|\psi_{1}^{\pm}\rangle$ and $|\text{E}_{\text{2g}}\rangle$) lifts the
degeneracy and block-diagonalizes the joint system into 2D subspace
$\omega^{\pm}$ and 2D subspace $\overline{\omega}^{\pm}$.

Figures 4d and 4e illustrate the 2 states,
$|\psi_{1}^{-}\rangle\otimes|\text{E}_{\text{2g,+}}\rangle$ and
$|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$ based on the
observations described above. In Figure 4d, when the ground state is
$|\psi_{0}^{+}\rangle$ (shown in light pink), the $J_{+}J_{-}$ or $J_{-}J_{+}$
operator excites the system to the $+$ branch,
$|\psi_{1}^{+}\rangle\otimes|\text{E}_{\text{2g,-}}\rangle$ (where the CEF
excitation corresponds to the $+$ branch shown in blue and the phononic
excitation corresponds to $-$ angular momentum shown in dark pink). Figure 4e
shows the same effect for $|\psi_{0}^{+}\rangle$ to
$|\psi_{1}^{-}\rangle\otimes|\text{E}_{\text{2g,+}}\rangle$. Because the
transition is degenerate for $|\psi^{\mp}_{0}\rangle\to\omega^{\pm}$ with
$J_{+}J_{-}$ or $J_{-}J_{+}$, this eigenspace, in principle, can be transduced
from a QSL ground state. If the QSL ground state is written as
$\sum_{j}\prod_{i}\psi_{0}^{\sigma_{ij}}$ – where $i$ is the site index, $j$
the configuration index, and $\sigma_{ij}=+$, $-$ is the branch – then
$J_{+}J_{-}$ or $J_{-}J_{+}$ should bring the system into
$\sum_{j}\prod_{i}\omega^{\sigma_{ij}}$ within the point spread volume of the
excitation beam, transducing the possible ground state entanglement to the VBS
$\omega^{\pm}$.

##  III Conclusion

We have shown here that the VBS can serve as a flexible probe of phonon
chirality by measuring an angular momentum change of $\delta J_{z}=\pm\hbar$
in a VBS due to coupling to degenerate $E_{\text{2g}}$ phonons and modeling
the interactions with a set of helicity-resolved selection rules between CEF
transitions. It is natural to ask whether this coupling provides alternative
pathways for spin relaxation. If the phonon band structure is modified by this
coupling, it may also be possible to access the resulting Berry curvature and
probe the robustness of potential topological invariants. Additionally, the
observation of photon- and CEF-mediated angular momentum transfer in a
material with $\mathcal{P}$ and $\mathcal{T}$ symmetry suggests that there may
be many more possible routes to the creation of chiral phonons by coupling to
electronic, spin, orbital degrees of freedom in a given system. This finding
thus creates a framework that may ultimately enable the transduction of
possible QSL ground states to experimentally accessible VBSs.

##  IV Methods

###  IV.1 Sample Details

Single crystals of NaYbSe2 were synthesized from NaCl (crystalline powder,
$99+\%$), Yb (ingot, $99.9\%$), and Se (powder, $99.999\%$) via the flux
method. The finely ground flux mixtures of NaCl, Yb, and Se with a molar ratio
of 20:1:2.4 were heated to 850 ∘C slowly in a vacuum. After two weeks, the
furnace was cooled to room temperature at a rate of 40∘C/hr. The laboratory-
grown single crystals of NaYbSe2 were separated from excess alkali halide flux
through washing with deionized water and isopropyl alcohol inside a fume-hood.
CsYbSe2 single crystals were grown using a related flux method that is
described more detail in previous work[36].

###  IV.2 Raman Spectroscopy

Variable temperature Raman spectra were acquired in a Montana Instruments
closed-cycle cryostat using an in-vacuum objective with numerical aperture of
0.85. A 1.5 mW, 532.03 nm continuous wave laser excited the sample in an out-
of-plane back scattering geometry (beam path $\|$ c). Rayleigh scattering was
minimized with either a set of 3 Optigrate volume Bragg gratings or a set of
Semrock RazorEdge ultrasteep dichroic and long-pass edge filters with cutoff
at 90 cm-1, and spectra were acquired with a 30 s exposure time.

Magnetic-field-dependent Raman spectra were acquired at a fixed temperature of
$T=$ 4 K for H $\|$ c in a customized Leiden dilution refrigerator with free
space optical access to the sample at the mixing chamber stage[58]. The
spectra were taken with an Andor Kymera 193 spectrograph (2400 line/mm
grating) and a Newton EMCCD DU970P-BV camera. The same laser and filters were
used as with the variable temperature measurements, with the laser power set
to 1.0 mW and a typical exposure time of 300 s per spectrum. Achromatic half-
wave plates and quarter-wave plates were mounted on rotators for automated
polarization control in both microscopes.

###  IV.3 Crystal Field Hamiltonian

For review of the approaches used here, see, for example, Bartolomé et al.
[59]. The CEF Hamiltonian for NaYbSe2 within the point charge approximation
[60, 21, 27, 34] that describes the ground state manifold $J=7/2$ is:

$\begin{split}H_{\text{CEF}}=&B_{2}^{0}\textbf{O}_{2}^{0}+B_{4}^{0}\textbf{O}_{4}^{0}+B_{4}^{3}\textbf{O}_{4}^{3}\\\
+&B_{6}^{0}\textbf{O}_{6}^{0}+B_{6}^{3}\textbf{O}_{6}^{3}+B_{6}^{6}\textbf{O}_{6}^{6}\
\end{split}$ (1)

The expected helicity dependence of the CEF excitations follows from the
selection rule bridging the two relevant states. The eigenstates for the
energy levels described in equation 1 are:

$\begin{split}|\psi^{+}_{\text{1, 2, 3}}\rangle=&-\alpha
e^{i\phi_{\alpha}}|\frac{7}{2},\frac{7}{2}\rangle+\beta|\frac{7}{2},\frac{1}{2}\rangle+\gamma
e^{-i\phi_{\gamma}}|\frac{7}{2},-\frac{5}{2}\rangle\\\ |\psi^{-}_{\text{1, 2,
3}}\rangle=&\alpha
e^{i\phi_{\alpha}}|\frac{7}{2},-\frac{7}{2}\rangle+\beta|\frac{7}{2},-\frac{1}{2}\rangle-\gamma
e^{i\phi_{\gamma}}|\frac{7}{2},\frac{5}{2}\rangle\\\ \end{split}$ (2)

and

$\begin{split}|\psi^{+}_{0}\rangle=&|\frac{7}{2},\frac{3}{2}\rangle\\\
|\psi^{-}_{0}\rangle=&|\frac{7}{2},-\frac{3}{2}\rangle\\\ \end{split}$ (3)

where $\alpha$, $\beta$, $\gamma$, $\phi_{\alpha}$, $\phi_{\beta}$, and
$\phi_{\gamma}$ $\in\mathbb{R}$ are determined by the CEF parameters
$B_{2}^{0}$, $B_{4}^{0}$, $B_{4}^{3}$, $B_{6}^{0}$, $B_{6}^{3}$ and
$B_{6}^{6}$. The eigenstates are shown in Figure 2a. Due to the 3-fold
symmetry of the Yb3+ environment, a 3-fold periodicity and hence angular
momentum folding, analogous to the the Umklapp process for linear momentum[53]
is expected. Without further constraints, the 6 parameters have enough degrees
of freedom to fit experimentally observed energy levels. More than one set of
CEF parameters that minimize the error may exist and the order of the
eigenstates may not be the same across those sets. For example, the little
group spanned by the special pair with single angular momentum eigenstate
$|\psi^{\pm}_{0}\rangle$ is assigned to CEF1 in Zhang et al. [21] and CEF2 in
Scheie et al. [39]. Schimidt et al. pointed out that they cannot be the ground
state due to observed in-plane field dependence at low temperatures[27]. In
addition to the Zeeman dependence
$H_{\textbf{B}}=-g_{J}\mu_{B}\textbf{J}_{z}\cdot\textbf{B}$, several
additional corrections have been considered. For example, Pocs et al. [34]
considered an $H_{\text{XXZ}}$ term. Zhang et al. considered anisotropic spin-
spin interactions [61].

In the context of magnetostrictive coupling with CEF modes, Callan et al.
considered an additional term
$H_{\text{me}}=-\sum_{\Gamma_{\nu}}\zeta(\Gamma_{\nu})\;u(\Gamma_{\nu})\;Q(\Gamma_{\nu})$
where $\zeta$ is the coupling strength, $u(\Gamma_{\nu})$ are phonon
operators, and $Q(\Gamma_{\nu})$ is the transformed phonon mode octupolar
operator on the CEF manifold. The Callen-Callen[62, 63] magnetoelastic
interaction is quadruplar. The quadruple operator (l = 2) is given by
A${}_{\text{1g}}$ \+ 2E${}_{\text{g}}$. Hence,

$\begin{split}H_{\text{me}}=&-\sum_{\Gamma_{\nu}}\zeta(\Gamma_{\nu})\;u(\Gamma_{\nu})\;Q(\Gamma_{\nu})\\\
Q(A_{\text{1g}})=&3J_{z}^{2}-J^{2}\\\
Q(E_{\text{g}}^{\text{I}})=&\alpha_{1}(J_{z}J_{x}+J_{x}J_{z})+\alpha_{2}(J_{y}J_{z}+J_{z}J_{y})\\\
Q(E_{\text{g}}^{\text{II}})=&\alpha_{3}(J_{x}^{2}-J_{y}^{2})+\alpha_{4}(J_{x}J_{y}+J_{y}J_{x})\\\
\end{split}$ (4)

For $Q(E_{\text{g}}^{\text{II}})$, since
$J_{x}^{2}-J_{y}^{2}=(J_{+}J_{+}+J_{-}J_{-})/2$, it induces
$\Delta\textbf{J}_{z}=\pm 2\hbar$, while for $J_{x}J_{y}+J_{y}J_{x}=iJ_{z}$,
$\Delta\textbf{J}_{z}=0$. For $Q(E_{\text{g}}^{\text{I}})$,
$J_{z}J_{x}+J_{x}J_{z}$ and $J_{z}J_{x}+J_{x}J_{z}$ both induces
$\Delta\textbf{J}_{z}=\pm 1\hbar$.

###  IV.4 Fitting Experimental Data to CEF Hamiltonian

To fit the experimental data, we used a general non-linear model normalization
procedure in Mathematica. To start, a general set of eigenvalues and
eigenvectors as a function of the CEF parameters is obtained by direct
diagonalization of the CEF Hamiltonian. The eigenvectors are sorted based on
their eigenvalues for each test parameter space. Then the cost function is
defined by the sum of the squared errors between transitions and the data,
with both inter-branch and intra-branch transitions considered. The selection
rules are added as a term in the cost function when necessary. The final CEF
parameters are the set that minimized the cost function.

###  IV.5 Bayesian Inference

We employed Bayesian inference to extract spectral parameters such as peak
positions and widths. This approach is based on Bayes’ rule:
$P(\theta|y)=\frac{P(y|\theta)P(\theta)}{P(y)}$. In the context of spectral
data, the priors $P(\theta)$ are the true distribution of the parameters such
as peak position, peak height, width, etc. The likelihood $P(y|\theta)$ is the
experimental data. The posterior $P(\theta|y)$ is the conditional distribution
of the experiment parameters given the experimental data. $P(y)$ is a
normalization factor. The distribution of the priors of the parameters are
assumed to be Gaussian or uniform. Additional noise, offset and slope are
added to capture background unrelated to the peak parameters. To carry out the
inference, the Hamiltonian Monte Carlo python package PyMC3[54] was used. A
hierarchical model was constructed to concurrently extract peak parameters
from a family of spectra such as the temperature dependence or magnetic field
dependence. A no U-Turns (NUTS) sampler was used with 4 chains with 3,000
samples per chain. It takes between 20 min to 3 hours for a peak over a
dataset to converge.

##  V Data availability

The data that support the findings of this study are available from the
corresponding authors upon request.

###### Acknowledgements.

The authors would like to acknowledge insightful discussion with Allen Scheie,
Xinshu Zhang, Yi Luo, Cristian Batista, Alan Tennant and Vyacheslav Bryantsev.
This research was sponsored by the U. S. Department of Energy, Office of
Science, Basic Energy Sciences, Materials Sciences and Engineering Division.
The first-principles phonon calculations and variable-temperature Raman
microscopy were performed at the Center for Nanophase Materials Sciences,
which is a U.S. Department of Energy Office of Science User Facility. SDW and
GP acknowledge support by the US Department of Energy, Office of Basic Energy
Sciences, Division of Materials Sciences and Engineering under award DE-
SC0017752. Postdoctoral research support was provided by the Intelligence
Community Postdoctoral Research Fellowship Program at the Oak Ridge National
Laboratory, administered by Oak Ridge Institute for Science and Education
through an interagency agreement between the U.S. Department of Energy and the
Office of the Director of National Intelligence.

##  VI Author contributions

All authors discussed the results thoroughly. Y.-Y. P., C. E. M., and B. J. L.
performed magneto-Raman measurements. L. L. performed Raman tensor analysis.
G. P. and J. X. grew the samples. Y.-Y. P, L. L. did the data analysis with
inputs from L. L., M. C., B. J. L.. X. L. and L. L. performed DFT calculation.
A. S. S., D. P., S. W. and B. J. L initiated and oversaw the project. Y.-Y.
P., L. L., B. J. L wrote most of the manuscript with contributions from all
authors.

##  VII Competing interests

The authors declare no competing interests.

##  VIII Additional Information

Supplementary information is available for this paper online. Correspondence
and requests for materials should be addressed to Y.-Y. P. or B. J. L.

## References

  * [1] Bloch, F.  Über die quantenmechanik der elektronen in kristallgittern.  _Zeitschrift für Physik_ 52, 555–600 (1929).  URL https://doi.org/10.1007/BF01339455.
  * [2] Chen, H., Zhang, W., Niu, Q. & Zhang, L.  Chiral phonons in two-dimensional materials.  _2D Materials_ 6, 012002 (2018).  URL https://doi.org/10.1088/2053-1583/aaf292.
  * [3] Zhang, L. & Niu, Q.  Chiral phonons at high-symmetry points in monolayer hexagonal lattices.  _Phys. Rev. Lett._ 115, 115502 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevLett.115.115502.
  * [4] Juneja, R. _et al._ Quasiparticle twist dynamics in non-symmorphic materials.  _Materials Today Physics_ 21, 100548 (2021).  URL https://www.sciencedirect.com/science/article/pii/S2542529321002091.
  * [5] Pandey, T., Polanco, C. A., Cooper, V. R., Parker, D. S. & Lindsay, L.  Symmetry-driven phonon chirality and transport in one-dimensional and bulk $\mathrm{B}{\mathrm{a}}_{3}\mathrm{N}$-derived materials.  _Phys. Rev. B_ 98, 241405 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.241405.
  * [6] Gunatilleke, W. D. C. B. _et al._ Intrinsic anharmonicity and thermal properties of ultralow thermal conductivity ${\mathrm{ba}}_{6}{\mathrm{sn}}_{6}{\mathrm{se}}_{13}$.  _Phys. Rev. Materials_ 5, 085002 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevMaterials.5.085002.
  * [7] Coh, S.  Classification of materials with phonon angular momentum and microscopic origin of angular momentum (2021).  eprint 1911.05064.
  * [8] Zhu, H. _et al._ Observation of chiral phonons.  _Science_ 359, 579–582 (2018).  URL https://www.science.org/doi/abs/10.1126/science.aar2711.  eprint https://www.science.org/doi/pdf/10.1126/science.aar2711.
  * [9] Yin, T. _et al._ Chiral phonons and giant magneto-optical effect in crbr3 2d magnet.  _Advanced Materials_ 33, 2101618 (2021).
  * [10] Manenkov, A. A. & Orbach, R.  _Spin-lattice relaxation in ionic solids_ (Harper & Row, 1966).
  * [11] Nakane, J. J. & Kohno, H.  Angular momentum of phonons and its application to single-spin relaxation.  _Phys. Rev. B_ 97, 174403 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.97.174403.
  * [12] Streib, S.  Difference between angular momentum and pseudoangular momentum.  _Phys. Rev. B_ 103, L100409 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.103.L100409.
  * [13] Strohm, C., Rikken, G. L. J. A. & Wyder, P.  Phenomenological evidence for the phonon hall effect.  _Phys. Rev. Lett._ 95, 155901 (2005).  URL https://link.aps.org/doi/10.1103/PhysRevLett.95.155901.
  * [14] Zhang, L. & Niu, Q.  Angular momentum of phonons and the einstein–de haas effect.  _Phys. Rev. Lett._ 112, 085503 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevLett.112.085503.
  * [15] Qin, T., Niu, Q. & Shi, J.  Energy magnetization and the thermal hall effect.  _Phys. Rev. Lett._ 107, 236601 (2011).  URL https://link.aps.org/doi/10.1103/PhysRevLett.107.236601.
  * [16] Chen, Z., Luo, J.-W. & Wang, L.-W.  Revealing angular momentum transfer channels and timescales in the ultrafast demagnetization process of ferromagnetic semiconductors.  _Proceedings of the National Academy of Sciences_ 116, 19258–19263 (2019).  URL https://www.pnas.org/content/116/39/19258.  eprint https://www.pnas.org/content/116/39/19258.full.pdf.
  * [17] Tsatsoulis, T. _et al._ Ultrafast demagnetization after femtosecond laser pulses: Transfer of angular momentum from the electronic system to magnetoelastic spin-phonon modes.  _Phys. Rev. B_ 93, 134411 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.93.134411.
  * [18] Jia, Y.-T. _et al._ Mott transition and superconductivity in quantum spin liquid candidate ${{NaYbSe}}_{2}$.  _Chinese Physics Letters_ 37, 097404 (2020).  URL http://dx.doi.org/10.1088/0256-307X/37/9/097404.
  * [19] Dai, P.-L. _et al._ Spinon fermi surface spin liquid in a triangular lattice antiferromagnet ${{NaYbSe}}_{2}$.  _Phys. Rev. X_ 11, 021044 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevX.11.021044.
  * [20] Ranjith, K. M. _et al._ Anisotropic field-induced ordering in the triangular-lattice quantum spin liquid ${{NaYbSe}}_{2}$.  _Phys. Rev. B_ 100, 224417 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.100.224417.
  * [21] Zhang, Z. _et al._ Crystalline electric field excitations in the quantum spin liquid candidate ${{NaYbSe}}_{2}$.  _Phys. Rev. B_ 103, 035144 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.103.035144.
  * [22] Zhang, Z. _et al._ Pressure induced metallization and possible unconventional superconductivity in spin liquid $naybse_{2}$ (2020).  eprint 2003.11479.
  * [23] Anderson, P.  Resonating valence bonds: A new kind of insulator?  _Materials Research Bulletin_ 8, 153–160 (1973).  URL https://www.sciencedirect.com/science/article/pii/0025540873901670.
  * [24] Knolle, J. & Moessner, R.  A field guide to spin liquids.  _Annual Review of Condensed Matter Physics_ 10, 451–472 (2019).  URL https://doi.org/10.1146/annurev-conmatphys-031218-013401.
  * [25] Savary, L. & Balents, L.  Quantum spin liquids: a review.  _Reports on Progress in Physics_ 80, 016502 (2016).  URL http://dx.doi.org/10.1088/0034-4885/80/1/016502.
  * [26] Barkeshli, M., Berg, E. & Kivelson, S.  Coherent transmutation of electrons into fractionalized anyons.  _Science_ 346, 722–725 (2014).
  * [27] Schmidt, B., Sichelschmidt, J., Ranjith, K. M., Doert, T. & Baenitz, M.  Yb delafossites: Unique exchange frustration of $4f$ spin-$\frac{1}{2}$ moments on a perfect triangular lattice.  _Phys. Rev. B_ 103, 214445 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.103.214445.
  * [28] Baenitz, M. _et al._ ${{NaYbS}}_{2}$: A planar spin-$\frac{1}{2}$ triangular-lattice magnet and putative spin liquid.  _Phys. Rev. B_ 98, 220409 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.220409.
  * [29] Sarkar, R. _et al._ Quantum spin liquid ground state in the disorder free triangular lattice ${{NaYbS}}_{2}$.  _Phys. Rev. B_ 100, 241116 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.100.241116.
  * [30] Ding, L. _et al._ Gapless spin-liquid state in the structurally disorder-free triangular antiferromagnet ${{NaYbO}}_{2}$.  _Phys. Rev. B_ 100, 144432 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.100.144432.
  * [31] Ranjith, K. M. _et al._ Field-induced instability of the quantum spin liquid ground state in the ${J}_{{eff}}=\frac{1}{2}$ triangular-lattice compound ${{NaYbO}}_{2}$.  _Phys. Rev. B_ 99, 180401 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.99.180401.
  * [32] Bordelon, M. M. _et al._ Field-tunable quantum disordered ground state in the triangular-lattice antiferromagnet naybo2.  _Nature Physics_ 15, 1058–1064 (2019).  URL https://doi.org/10.1038/s41567-019-0594-5.
  * [33] Xing, J. _et al._ Field-induced magnetic transition and spin fluctuations in the quantum spin-liquid candidate ${{CsYbSe}}_{2}$.  _Phys. Rev. B_ 100, 220407 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.100.220407.
  * [34] Pocs, C. A. _et al._ Systematic fitting of crystal-field levels and accurate extraction of quantum magnetic models in triangular-lattice delafossites (2021).
  * [35] Xie, T. _et al._ Field-induced spin excitations in the spin-1/2 triangular-lattice antiferromagnet csybse2 (2021).  eprint 2106.12451.
  * [36] Xing, J. _et al._ Crystal synthesis and frustrated magnetism in triangular lattice ${CsRESe_{2}}$ (re = la–lu): Quantum spin liquid candidates ${CsCeSe_{2}}$ and ${CsYbSe_{2}}$.  _ACS Materials Letters_ 2, 71–75 (2020).  URL https://doi.org/10.1021/acsmaterialslett.9b00464.
  * [37] Pai, Y.-Y. _et al._ Mesoscale interplay between phonons and crystal electric field excitations in quantum spin liquid candidate csybse2.  _J. Mater. Chem. C_ – (2022).  URL http://dx.doi.org/10.1039/D1TC05934C.
  * [38] Xing, J., Sanjeewa, L. D., May, A. F. & Sefat, A. S.  Synthesis and anisotropic magnetism in quantum spin liquid candidates $A$YbSe2 ($A$ = K and Rb) (2021).  eprint 2109.00384.
  * [39] Scheie, A. O. _et al._ Witnessing quantum criticality and entanglement in the triangular antiferromagnet kybse2 (2021).  eprint 2109.11527.
  * [40] Koningstein, J. A.  Electronic raman spectra. i. raman transitions of trivalent ytterbium, europium, and neodymium in yttrium gallium garnet.  _The Journal of Chemical Physics_ 46, 2811–2816 (1967).  URL https://doi.org/10.1063/1.1841119.  eprint https://doi.org/10.1063/1.1841119.
  * [41] Petit, V., Camy, P., Doualan, J.-L., Portier, X. & Moncorgé, R.  Spectroscopy of ${\text{yb}}^{3+}:{\text{caf}}_{2}$: From isolated centers to clusters.  _Phys. Rev. B_ 78, 085131 (2008).  URL https://link.aps.org/doi/10.1103/PhysRevB.78.085131.
  * [42] Čermák, P. _et al._ Magnetoelastic hybrid excitations in ${CeAuAl_{3}}$.  _Proceedings of the National Academy of Sciences_ 116, 6695–6700 (2019).
  * [43] Sethi, A., Slimak, J. E., Kolodiazhnyi, T. & Cooper, S. L.  Emergent vibronic excitations in the magnetodielectric regime of ${{Ce}}_{2}{{O}}_{3}$.  _Phys. Rev. Lett._ 122, 177601 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.177601.
  * [44] Loewenhaupt, M., Rainford, B. D. & Steglich, F.  Dynamic jahn-teller effect in a rare-earth compound: Ce${{Al}}_{2}$.  _Phys. Rev. Lett._ 42, 1709–1712 (1979).  URL https://link.aps.org/doi/10.1103/PhysRevLett.42.1709.
  * [45] Güntherodt, G., Jayaraman, A., Batlogg, G., Croft, M. & Melczer, E.  Raman scattering from coupled phonon and electronic crystal-field excitations in Ce${{Al}}_{2}$.  _Phys. Rev. Lett._ 51, 2330–2332 (1983).  URL https://link.aps.org/doi/10.1103/PhysRevLett.51.2330.
  * [46] Thalmeier, P. & Fulde, P.  Bound state between a crystal-field excitation and a phonon in Ce${{Al}}_{2}$.  _Phys. Rev. Lett._ 49, 1588–1591 (1982).  URL https://link.aps.org/doi/10.1103/PhysRevLett.49.1588.
  * [47] Adroja, D. _et al._ Vibron quasibound state in the noncentrosymmetric tetragonal heavy-fermion compound $CeCuAl_{3}$.  _Physical review letters_ 108, 216402 (2012).
  * [48] Aksenov, V., Goremychkin, E., Mühle, E., Frauenheim, T. & Bührer, W.  Coupled quadrupole-phonon excitations: Inelastic neutron scattering on van vleck paramagnet ${PrNi_{5}}$.  _Physica B+C_ 120, 310–313 (1983).  URL https://www.sciencedirect.com/science/article/pii/0378436383903984.
  * [49] Heyen, E. T., Wegerer, R. & Cardona, M.  Coupling of phonons to crystal-field excitations in ${{NdBa}}_{2}$${{Cu}}_{3}$${{O}}_{7{-}{\delta}}$.  _Phys. Rev. Lett._ 67, 144–147 (1991).  URL https://link.aps.org/doi/10.1103/PhysRevLett.67.144.
  * [50] Gaudet, J. _et al._ Magnetoelastically induced vibronic bound state in the spin-ice pyrochlore ${{Ho}}_{2}{{Ti}}_{2}{{O}}_{7}$.  _Phys. Rev. B_ 98, 014419 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.014419.
  * [51] Güntherodt, G., Jayaraman, A., Anastassakis, E., Bucher, E. & Bach, H.  Effect of configuration crossover on the electronic Raman scattering by $4f$ multiplets.  _Phys. Rev. Lett._ 46, 855–858 (1981).  URL https://link.aps.org/doi/10.1103/PhysRevLett.46.855.
  * [52] Zhu, Z. H. _et al._ Fluctuating magnetic droplets immersed in a sea of quantum spin liquid (2021).  eprint 2112.06523.
  * [53] Tatsumi, Y., Kaneko, T. & Saito, R.  Conservation law of angular momentum in helicity-dependent raman and rayleigh scattering.  _Phys. Rev. B_ 97, 195444 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.97.195444.
  * [54] Salvatier, J., Wiecki, T. V. & Fonnesbeck, C.  Probabilistic programming in python using pymc3.  _PeerJ Computer Science_ 2, e55 (2016).  URL https://doi.org/10.7717/peerj-cs.55.
  * [55] Deng, B. & Ibers, J. A.  CsYbSe2.  _Acta Crystallographica Section E_ 61, i15–i17 (2005).  URL https://doi.org/10.1107/S1600536805001157.
  * [56] Alpeggiani, F., Bliokh, K. Y., Nori, F. & Kuipers, L.  Electromagnetic helicity in complex media.  _Phys. Rev. Lett._ 120, 243605 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.243605.
  * [57] Axe, J. D.  Two-photon processes in complex atoms.  _Phys. Rev._ 136, A42–A45 (1964).  URL https://link.aps.org/doi/10.1103/PhysRev.136.A42.
  * [58] Lawrie, B. J., Feldman, M., Marvinney, C. E. & Pai, Y. Y.  Free-space confocal magneto-optical spectroscopies at millikelvin temperatures (2021).  eprint 2103.06851.
  * [59] Bartolomé, E., Arauzo, A., Luzón, J., Bartolomé, J. & Bartolomé, F.  Chapter 1 - magnetic relaxation of lanthanide-based molecular magnets.  vol. 26 of _Handbook of Magnetic Materials_ , 1–289 (Elsevier, 2017).  URL https://www.sciencedirect.com/science/article/pii/S1567271917300033.
  * [60] Jensen, J., Mackintosh, A. & Mackintosh, B.  _Rare Earth Magnetism: Structures and Excitations_.  International Series of Monographs on Physics (Clarendon Press, 1991).  URL https://books.google.com/books?id=LbTvAAAAMAAJ.
  * [61] Zhang, Z. _et al._ Effective magnetic hamiltonian at finite temperatures for rare-earth chalcogenides.  _Phys. Rev. B_ 103, 184419 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.103.184419.
  * [62] Callen, E. & Callen, H. B.  Magnetostriction, forced magnetostriction, and anomalous thermal expansion in ferromagnets.  _Phys. Rev._ 139, A455–A471 (1965).  URL https://link.aps.org/doi/10.1103/PhysRev.139.A455.
  * [63] Lovesey, S. W. & Staub, U.  Magnetoelastic model for the relaxation of lanthanide ions in ${\mathrm{yba}}_{2}{\mathrm{cu}}_{3}{\mathrm{o}}_{7-\delta}$ observed by neutron scattering.  _Phys. Rev. B_ 61, 9130–9139 (2000).  URL https://link.aps.org/doi/10.1103/PhysRevB.61.9130.


and ${\mathbf{r}}_{2}$, together with
$||({\mathbf{a}}\times{\mathbf{s}}_{1})||=\sin(\eta_{{\mathbf{a}}{\mathbf{s}}_{1}})$
and
$||({\mathbf{s}}_{2}\times{\mathbf{b}})||=\sin(\eta_{{\mathbf{s}}_{2}{\mathbf{b}}})$
from the normalizations of cross products, (A.64) simplifies to

$\cos(\eta_{{\mathbf{u}}{\mathbf{v}}})=\\{({\mathbf{a}}\cdot{\mathbf{s}}_{1})({\mathbf{s}}_{2}\cdot{\mathbf{b}})-({\mathbf{a}}\times{\mathbf{s}}_{1})\cdot({\mathbf{s}}_{2}\times{\mathbf{b}})\\}.$
(A.66)

This can be further simplified by using the vector identity

$({\mathbf{a}}\times{\mathbf{s}}_{1})\cdot({\mathbf{s}}_{2}\times{\mathbf{b}})=({\mathbf{a}}\cdot{\mathbf{s}}_{2})({\mathbf{s}}_{1}\cdot{\mathbf{b}})-({\mathbf{a}}\cdot{\mathbf{b}})({\mathbf{s}}_{1}\cdot{\mathbf{s}}_{2}).$
(A.67)

Consequently, we arrive at the angle $\eta_{{\mathbf{u}}{\mathbf{v}}}$
expressed in Eq. (99) of Section 5:

$\eta_{{\mathbf{u}}{\mathbf{v}}}({\mathbf{s}}_{1},\,{\mathbf{s}}_{2})=\cos^{-1}\big{\\{}({\mathbf{a}}\cdot{\mathbf{s}}_{1})({\mathbf{s}}_{2}\cdot{\mathbf{b}})-({\mathbf{a}}\cdot{\mathbf{s}}_{2})({\mathbf{s}}_{1}\cdot{\mathbf{b}})+({\mathbf{a}}\cdot{\mathbf{b}})({\mathbf{s}}_{1}\cdot{\mathbf{s}}_{2})\big{\\}}.$
(A.68)

On the other hand, solving equation (A.65) for
${\mathbf{r}}_{0}({\mathbf{s}}_{1},\,{\mathbf{s}}_{2})$ gives

${\mathbf{r}}_{0}({\mathbf{s}}_{1},\,{\mathbf{s}}_{2})=\frac{{\mathbf{r}}_{1}\cos\left(\,\eta_{{\mathbf{s}}_{2}{\mathbf{b}}}\right)\sin\left(\,\eta_{{\mathbf{a}}{\mathbf{s}}_{1}}\right)+{\mathbf{r}}_{2}\cos\left(\,\eta_{{\mathbf{a}}{\mathbf{s}}_{1}}\right)\sin\left(\,\eta_{{\mathbf{s}}_{2}{\mathbf{b}}}\right)-({\mathbf{r}}_{1}\times{\mathbf{r}}_{2})\sin\left(\,\eta_{{\mathbf{a}}{\mathbf{s}}_{1}}\right)\sin\left(\,\eta_{{\mathbf{s}}_{2}{\mathbf{b}}}\right)}{\sin\left(\,\eta_{{\mathbf{a}}{\mathbf{b}}}\right)}.$
(A.69)

Again using the definitions (9) and (13) of the axis vectors
${\mathbf{r}}_{1}$ and ${\mathbf{r}}_{2}$ and
$||({\mathbf{a}}\times{\mathbf{s}}_{1})||=\sin(\eta_{{\mathbf{a}}{\mathbf{s}}_{1}})$
and
$||({\mathbf{s}}_{2}\times{\mathbf{b}})||=\sin(\eta_{{\mathbf{s}}_{2}{\mathbf{b}}})$
from the normalizations of the cross products, (A.69) simplifies to

${\mathbf{r}}_{0}({\mathbf{s}}_{1},\,{\mathbf{s}}_{2})=\frac{({\mathbf{a}}\cdot{\mathbf{s}}_{1})({\mathbf{s}}_{2}\times{\mathbf{b}})+({\mathbf{s}}_{2}\cdot{\mathbf{b}})({\mathbf{a}}\times{\mathbf{s}}_{1})-({\mathbf{a}}\times{\mathbf{s}}_{1})\times({\mathbf{s}}_{2}\times{\mathbf{b}})}{\sin\left(\,\eta_{{\mathbf{a}}{\mathbf{b}}}\right)},$
(A.70)

which is Eq. (100) of Section 5. With the above explicit expressions of
${\mathbf{r}}_{0}({\mathbf{s}}_{1},\,{\mathbf{s}}_{2})$, for
${\mathbf{s}}_{1}=\,{\mathbf{s}}_{2}$ the quaternion (A.63) reduces to

${\mathbf{q}}(\eta_{{\mathbf{a}}{\mathbf{b}}},\,{\mathbf{r}}_{0})=\cos(\eta_{{\mathbf{a}}{\mathbf{b}}})+(I_{3}{\mathbf{r}}_{0})\,\sin(\eta_{{\mathbf{a}}{\mathbf{b}}}),$
(A.71)

because, as evident from (A.68), the angle $\eta_{{\mathbf{u}}{\mathbf{v}}}$
reduces to $\eta_{{\mathbf{a}}{\mathbf{b}}}$ if
${\mathbf{s}}_{1}=\,{\mathbf{s}}_{2}$ holds, as required by conservation of
angular momentum. We thus see that the product
${\mathbf{q}}(\eta_{{\mathbf{a}}{\mathbf{s}}_{1}},\,{\mathbf{r}}_{1})\,{\mathbf{q}}(\eta_{{\mathbf{s}}_{2}{\mathbf{b}}},\,{\mathbf{r}}_{2})$
in Eq. (117) reduces to the quaternion
${\mathbf{q}}(\eta_{{\mathbf{a}}{\mathbf{b}}},\,{\mathbf{r}}_{0})$ in Eq.
(118) for ${\mathbf{s}}_{1}=\,{\mathbf{s}}_{2}$.

## References

  * 1 Einstein A, Podolsky B, Rosen N. Can quantum-mechanical description of physical reality be considered complete? Phys Rev. 1935; 47:777-780. DOI: 10.1103/PhysRev.47.777.
  * 2 Bell JS. On the Einstein-Podolsky-Rosen paradox. Physics. 1964; 1(3):195-200.
  * 3 Aspect A, Grangier P, Roger G. Experimental realization of Einstein-Podolsky-Rosen-Bohm Gedankenexperiment: a new violation of Bell’s inequalities. Phys Rev Lett. 1982; 49:91-94. DOI: 10.1103/PhysRevLett.49.91.
  * 4 Christian J. Macroscopic observability of spinorial sign changes under ${2\pi}$ rotations. Int J Theor Phys. 2015; 54:20-46. DOI: 10.1007/s10773-014-2412-2.
  * 5 Christian J. Quantum correlations are weaved by the spinors of the Euclidean primitives. R Soc Open Sci. 2018; 5(1): 180526. DOI: 10.1098/rsos.180526.
  * 6 Christian J. Bell’s theorem versus local realism in a quaternionic model of physical space. IEEE Access. 2019; 7:133388-133409. DOI: 10.1109/ACCESS.2019.2941275.
  * 7 Christian J. Dr. Bertlmann’s socks in the quaternionic world of ambidextral reality. IEEE Access. 2020; 8:191028-191048. DOI: 10.1109/ACCESS.2020.3031734.
  * 8 Pais A. ‘Subtle is the Lord’: the science and the life of Albert Einstein. Oxford, UK: Clarendon Press; 1982, 460.
  * 9 Doran C, Lasenby AN. Geometric Algebra for Physicists. Cambridge, UK: Cambridge University Press; 2003.
  * 10 Christian J. Exactly soluble sector of quantum gravity. Phys Rev D. 1997; 56:4844-4877. DOI: 10.1103/PhysRevD.56.4844.
  * 11 d’Inverno R. Introducing Einstein’s Relativity. Oxford, UK: Clarendon Press; 1992.
  * 12 Di Valentino E, Melchiorri A, Silk J. Planck evidence for a closed Universe and a possible crisis for cosmology. Nature Astronomy. 2019; DOI: 10.1038/s41550-019-0906-9.
  * 13 Handley W. Curvature tension: evidence for a closed universe. Phys. Rev. D. 2019: 103:L041301. DOI: 10.1103/PhysRevD.103.L041301.
  * 14 Bell JS. Speakable and Unspeakable in Quantum Mechanics. Cambridge, UK: Cambridge University Press; 1987.
  * 15 Christian J. Disproof of Bell’s Theorem: Illuminating the Illusion of Entanglement. Second Edition. Boca Raton, Florida: Brownwalker Press; 2014, Chapter 6.
  * 16 Lasenby AN, A 1d up approach to conformal geometric algebra: application in line fitting and quantum mechanics. Adv Appl Clifford Algebras. 2020; 30:22. DOI: 10.1007/s00006-020-1046-0.
  * 17 Nakahara M. Geometry, Topology and Physics. New York, New York: Adam Hilger; 1990.
  * 18 Bloch WG. The Unimaginable Mathematics of Borges’ Library of Babel. Oxford, UK: Oxford University Press; 2008.
  * 19 Milnor JW. Topology from the Differentiable Viewpoint. Princeton, New Jersey: Princeton University Press; 1997.
  * 20 Dray T, Manogue CA. The Geometry of the Octonions. Singapore, Singapore: World Scientific; 2015, Chapter 5.
  * 21 Sobczyk G. The hyperbolic number plane. The College Mathematics Journal. 1995; 26:268-280. DOI: 10.1080/07468342.1995.11973712.
  * 22 Christian J. Eight-dimensional octonion-like but associative normed division algebra. 2019; https://arxiv.org/abs/1908.06172. Date of Access: 8 Jan 2021.
  * 23 Hurwitz A. Über die Composition der quadratischen Formen von beliebig vielen Variabeln. Nachr. Ges. Wiss. Göttingen. 1898; 1898:309-316.
  * 24 Hurwitz A. Über die Komposition der quadratischen Formen. Math Ann. 1923; 88(1–2):1-25.
  * 25 Baez JC. The octonions. Bull Amer Math Soc. 2002; 39:145-205. DOI: 10.1090/S0273-0979-01-00934-X.
  * 26 Lasenby AN. Problems with handedness and the role of normed division algebras in some recent papers on Bell’s Theorem. 2021; Lecture presented at the AGACSE 2021 conference, Brno, Czech Republic.
  * 27 Christian J. Reply to ‘A 1d up approach to conformal geometric algebra: applications in line fitting and quantum mechanics. ResearchGate. 2020; DOI: 10.13140/RG.2.2.32002.09929/1. Date of Access: 1 May 2020.
  * 28 Christian J. Reply to “comment on ‘Dr. Bertlmann’s socks in the quaternionic world of ambidextral reality”’. IEEE Access. 2021; 9:72161-72171. DOI: 10.1109/ACCESS.2021.3076449.
  * 29 Peres A. Quantum Theory: Concepts and Methods. Dordrecht, The Netherlands: Kluwer; 1993, 160.
  * 30 Christian J. Reply to “Comment on ‘Bell’s theorem versus local realism in a quaternionic model of physical space”’. IEEE Access. 2022; 10:14429-14439. DOI: 10.1109/ACCESS.2022.3146421.
  * 31 Diether III CF. Simulation for "Local origins of quantum correlations rooted in geometric algebra." https://notebookarchive.org/2022-04-2tqt8f0. Date of Access: 12 April 2022.
  * 32 Christian J. Symmetric derivation of the singlet correlations within a quaternionic 3-sphere. 2022; https://doi.org/10.48550/arXiv.2204.10288. Date of Access: 19 Apr 2022.